Solved the forcing function boundary problem.

acc/mhd_solver/stencil_process.sps

@@ -275,8 +275,7 @@ forcing(int3 globalVertexIdx, Scalar dt)
     Vector force   = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase);
 
     //Scaling N = magnitude*cs*sqrt(k*cs/dt)  * 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(DCONST_REAL(AC_kaver)*cs); 
-    const Scalar NN = cs*sqrt(kk*cs); // kk should be the average k!!!
     //MV: Like in the Pencil Code. I don't understandf the logic here. 
     force.x = sqrt(dt)*NN*force.x;
     force.y = force.y;

analysis/python/jupyter/notebook_example.ipynb

@@ -24,7 +24,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "meshdir  = \"/scratch/data/mvaisala/forcingtest/\""
+    "meshdir  = \"/scratch/data/mvaisala/forcingtest/\"\n",
+    "#meshdir  = \"/scratch/data/mvaisala/normaltest/\""
    ]
   },
   {
@@ -34,7 +35,7 @@
    "outputs": [],
    "source": [
     "#imesh = 30000\n",
-    "imesh = 30000\n",
+    "imesh = 100\n",
     "mesh = ad.read.Mesh(imesh, fdir=meshdir)"
    ]
   },
@@ -101,15 +102,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "vis.lineplot.plot_ts(ts, lnrho=1, uutot=1, uux=1, uuy=1, ss=1)"
+    "vis.lineplot.plot_ts(ts, lnrho=1, uutot=1, uux=1, uuy=1, uuz=1, ss=1)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

include/astaroth.h

@@ -165,6 +165,7 @@ extern "C" {
         FUNC(AC_k_forcex),\
         FUNC(AC_k_forcey),\
         FUNC(AC_k_forcez),\
+        FUNC(AC_kaver),\
         FUNC(AC_ff_hel_rex),\
         FUNC(AC_ff_hel_rey),\
         FUNC(AC_ff_hel_rez),\
@@ -409,7 +410,7 @@ AcResult acSynchronize(void);
 /** Tool for loading forcing vector information into the device memory
  */
 AcResult acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force, const AcReal3 ff_hel_re,
-                      const AcReal3 ff_hel_im, const AcReal forcing_phase);
+                      const AcReal3 ff_hel_im, const AcReal forcing_phase, const AcReal kaver);
  
 /* End extern "C" */
 #ifdef __cplusplus

src/core/astaroth.cu

@@ -508,7 +508,7 @@ acSynchronize(void)
 //Tool for loading forcing vector information into the device memory
 AcResult
 acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force, const AcReal3 ff_hel_re, 
-             const AcReal3 ff_hel_im, const AcReal forcing_phase)
+             const AcReal3 ff_hel_im, const AcReal forcing_phase, const AcReal kaver)
 {
 
     for (int i = 0; i < num_devices; ++i) {
@@ -526,6 +526,9 @@ acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force, const AcReal
         loadDeviceConstant(devices[i], AC_ff_hel_imx, ff_hel_im.x);
         loadDeviceConstant(devices[i], AC_ff_hel_imy, ff_hel_im.y);
         loadDeviceConstant(devices[i], AC_ff_hel_imz, ff_hel_im.z);
+
+        loadDeviceConstant(devices[i], AC_kaver, kaver);
+
     }
  
     return AC_SUCCESS;

src/standalone/simulation.cc

@@ -132,6 +132,13 @@ helical_forcing_k_generator(const AcReal kmax, const AcReal kmin)
                         kk*sin(theta)*sin(phi),
                         kk*cos(theta)         };
 
+    //printf("k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x, k_force.y, k_force.z);
+
+    //Round the numbers. In that way k(x/y/z) will get complete waves. 
+    k_force.x = round(k_force.x); k_force.y = round(k_force.y); k_force.z = round(k_force.z);
+
+    //printf("After rounding --> k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x, k_force.y, k_force.z);
+
     return k_force;
 }
 
@@ -156,7 +163,8 @@ helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force)
 
 //PC Manual Eq. 223 
 static inline void
-helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force, const AcReal3 e_force, const AcReal relhel)
+helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force, 
+                               const AcReal3 e_force, const AcReal relhel)
 {
     // k dot e 
     AcReal3 kdote;
@@ -191,7 +199,7 @@ helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcR
     //                       k_cross_k_cross_e.y/denominator,
     //                       k_cross_k_cross_e.z/denominator};
 
-    // See PC forcing.f90 rel_hel()
+    // See PC forcing.f90 forcing_hel_both()
     *ff_hel_re = (AcReal3){kabs*k_cross_e.x/denominator, 
                            kabs*k_cross_e.y,
                            kabs*k_cross_e.z};
@@ -377,6 +385,13 @@ run_simulation(void)
     AcReal bin_save_t = mesh_info.real_params[AC_bin_save_t];
     AcReal bin_crit_t = bin_save_t;
 
+    //Placeholders until determined properly
+    AcReal  magnitude = 0.05;
+    AcReal  relhel = 0.5;
+    AcReal  kmin = 0.8;
+    AcReal  kmax = 1.2;
+
+    AcReal kaver = (kmax - kmin)/AcReal(2.0);
 
     /* initialize random seed: */
     srand (312256655);
@@ -387,12 +402,7 @@ run_simulation(void)
         const AcReal dt   = host_timestep(umax, mesh_info);
 
 #if LFORCING
-        //Generate a forcing vectors before canculating an integration step. 
+        //Generate a forcing vector before canculating an integration step. 
-        //Placeholders until determined properly
-        AcReal  magnitude = 0.05;
-        AcReal  relhel = 0.5;
-        AcReal  kmin = 0.9999;
-        AcReal  kmax = 1.0001;
    
         // Generate forcing wave vector k_force
         AcReal3 k_force;
@@ -412,7 +422,7 @@ run_simulation(void)
 
         AcReal3 ff_hel_re, ff_hel_im;
         helical_forcing_special_vector(&ff_hel_re, &ff_hel_im, k_force, e_force, relhel);
-        acForcingVec(magnitude, k_force, ff_hel_re,ff_hel_im, phase);
+        acForcingVec(magnitude, k_force, ff_hel_re, ff_hel_im, phase, kaver);
 #endif
 
         acIntegrate(dt);