From 450526b8e17165cff606b12fc22bdb0508bc1a8b Mon Sep 17 00:00:00 2001
From: Miikka Vaisala <mvaisala@asiaa.sinica.edu.tw>
Date: Fri, 11 Sep 2020 13:23:01 +0800
Subject: [PATCH] B-field calculation to the DSL file.

Still needs corresponding things elsewhere...
---
 acc/mhd_solver/stencil_kernel.ac | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/acc/mhd_solver/stencil_kernel.ac b/acc/mhd_solver/stencil_kernel.ac
index 3aec600..a0b3143 100644
--- a/acc/mhd_solver/stencil_kernel.ac
+++ b/acc/mhd_solver/stencil_kernel.ac
@@ -8,6 +8,7 @@
 #define LFORCING (0)
 #define LUPWD (0)
 #define LSINK (0)
+#define LBFIELD (1)
 
 #define AC_THERMAL_CONDUCTIVITY (0.001) // TODO: make an actual config parameter
 #define H_CONST (0)                     // TODO: make an actual config parameter
@@ -74,6 +75,9 @@ uniform Scalar AC_ampl_uu;
 uniform Scalar AC_angl_uu;
 uniform Scalar AC_lnrho_edge;
 uniform Scalar AC_lnrho_out;
+uniform Scalar AC_ampl_aa;
+uniform Scalar AC_init_k_wave;
+uniform Scalar AC_init_sigma_hel;
 //  Forcing parameters. User configured.
 uniform Scalar AC_forcing_magnitude;
 uniform Scalar AC_relhel;
@@ -133,6 +137,12 @@ uniform ScalarField VTXBUF_LNRHO;
 uniform ScalarField VTXBUF_ACCRETION;
 #endif
 
+#if LBFIELD
+uniform ScalarField BFIELDX;
+uniform ScalarField BFIELDY;
+uniform ScalarField BFIELDZ;
+#endif
+
 #if LUPWD
 
 Preprocessed Scalar
@@ -629,6 +639,15 @@ forcing(int3 globalVertexIdx, Scalar dt)
 }
 #endif // LFORCING
 
+#if LBFIELD
+// Calculate the B-field for VTXBUFF
+Device Vector
+get_bfield(in VectorField aa)
+{
+    return curl(aa);
+}
+#endif
+
 // Declare input and output arrays using locations specified in the
 // array enum in astaroth.h
 in ScalarField lnrho(VTXBUF_LNRHO);
@@ -657,6 +676,11 @@ in ScalarField accretion(VTXBUF_ACCRETION);
 out ScalarField out_accretion(VTXBUF_ACCRETION);
 #endif
 
+#if LBFIELD
+out VectorField b_out(BFIELDX, BFIELDY, BFIELDZ);
+#endif
+
+
 Kernel void
 solve()
 {
@@ -691,4 +715,11 @@ solve()
         out_accretion = out_accretion * AC_dsx * AC_dsy * AC_dsz; // unit is now mass!
     }
 #endif
+
+#if LBFIELD
+    if (step_number == 2) {
+        b_out = get_bfield(aa); 
+    }
+#endif 
+
 }