Cleaning the forcing part.

acc/mhd_solver/stencil_assembly.sas

@@ -1,5 +1,5 @@
 
-#define LUPWD (0)
+#define LUPWD (1)
 
 
 Preprocessed Scalar

acc/mhd_solver/stencil_process.sps

@@ -3,7 +3,7 @@
 #define LTEMPERATURE (0)
 #define LGRAVITY (0)
 #define LFORCING (1)
-#define LUPWD (0)
+#define LUPWD (1)
 
 
 // Declare uniforms (i.e. device constants)
@@ -204,6 +204,19 @@ heat_transfer(in Vector uu, in Scalar lnrho, in Scalar tt)
 #endif
 
 #if LFORCING
+
+Vector
+simple_vortex_forcing(Vector a, Vector b, Scalar magnitude)
+{
+    return magnitude * cross(normalized(b - a), (Vector){0, 0, 1}); // Vortex   
+}
+
+Vector
+simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude)
+{
+    return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow   
+}
+
 Vector
 forcing(int3 globalVertexIdx)
 {
@@ -215,8 +228,11 @@ forcing(int3 globalVertexIdx)
                         (globalVertexIdx.z - nz_min) * dsz}; // sink (current index)
 
     Scalar magnitude = 0.05;
-    // Vector c = magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
-    Vector c = magnitude * cross(normalized(b - a), (Vector){0, 0, 1}); // Vortex
+
+    //Determine that forcing funtion type at this point. 
+    Vector c = simple_vortex_forcing(a, b, magnitude);
+    //Vector c = simple_outward_flow_forcing(a, b, magnitude);
+
     if (is_valid(c)) { return c; }
     else             { return (Vector){0, 0, 0}; }
 }

analysis/python/jupyter/notebook_example.ipynb

@@ -33,7 +33,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#imesh = 16000\n",
+    "#imesh = 30000\n",
     "imesh = 30000\n",
     "mesh = ad.read.Mesh(imesh, fdir=meshdir)"
    ]
@@ -76,7 +76,7 @@
     "vis.slices.plot_3(mesh, mesh.lnrho, title = r'$\\ln_\\rho$')\n",
     "vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$\\rho$')\n",
     "vis.slices.plot_3(mesh, np.exp(mesh.lnrho), title = r'$N_\\mathrm{col}$', slicetype = 'sum')\n",
-    "#vis.slices.plot_3(mesh, mesh.entropy, title = r'$s$')\n"
+    "vis.slices.plot_3(mesh, mesh.ss, title = r'$s$')\n"
    ]
   },
   {