Trying to calculate the forcing scaling.

Causes nans very quickly. Will need to look closer tomorrow again.

acc/mhd_solver/stencil_process.sps

@@ -241,8 +241,13 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto
     return force;
 }
 
+#if LENTROPY
 Vector
-forcing(int3 globalVertexIdx)
+forcing(int3 globalVertexIdx, Scalar dt, in Scalar lnrho, in Scalar ss)
+#else
+Vector
+forcing(int3 globalVertexIdx, Scalar dt)
+#endif
 {
     Vector a = Scalar(.5) * (Vector){globalGrid.n.x * dsx,
                                      globalGrid.n.y * dsy,
@@ -250,7 +255,12 @@ forcing(int3 globalVertexIdx)
     Vector xx = (Vector){(globalVertexIdx.x - nx_min) * dsx,
                         (globalVertexIdx.y - ny_min) * dsy,
                         (globalVertexIdx.z - nz_min) * dsz}; // sink (current index)
-
+#if LENTROPY
+    const Scalar cs2 = cs2_sound * exp(gamma * value(ss) / cp_sound + (gamma - 1) * (value(lnrho) - LNRHO0));
+#else
+    const Scalar cs2 = cs2_sound;
+#endif
+    const Scalar cs = sqrt(cs2);
 
     //Placeholders until determined properly
     Scalar magnitude = DCONST_REAL(AC_forcing_magnitude);
@@ -265,6 +275,13 @@ forcing(int3 globalVertexIdx)
     //Vector force = simple_outward_flow_forcing(a, xx, magnitude);
     Vector force   = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase);
 
+    //Scaling N = magnitude*cs*sqrt(k*cs/dt)
+    const Scalar kk = sqrt(k_force.x*k_force.x + k_force.y*k_force.y + k_force.z*k_force.z);
+    const Scalar NN = cs*sqrt((kk*cs)/dt);
+    force.x = NN*force.x;
+    force.y = NN*force.y;
+    force.z = NN*force.z;
+
     if (is_valid(force)) { return force; }
     else                 { return (Vector){0, 0, 0}; }
 }
@@ -314,7 +331,11 @@ solve(Scalar dt) {
 
     #if LFORCING
     if (step_number == 2) {
-        out_uu = out_uu + dt * forcing(globalVertexIdx);
+        #if LENTROPY
+        out_uu = out_uu + forcing(globalVertexIdx, dt, lnrho, ss);
+        #else
+        out_uu = out_uu + forcing(globalVertexIdx, dt);
+        #endif
     }
     #endif
 }

src/standalone/simulation.cc

@@ -386,7 +386,7 @@ run_simulation(void)
         AcReal  kmax = 1.7;
     
         // Generate forcing wave vector k_force
-        AcReal3 k_force;// = (AcReal3){0.0, 2.0, 0.0};
+        AcReal3 k_force;
         k_force = helical_forcing_k_generator(kmax, kmin);
 
         //Randomize the phase