Merge branch 'master' into node_device_interface_revision_07-23

3rdparty/.gitignore

@@ -0,0 +1,6 @@
+# Ignore everything
+*
+
+# Except:
+!.gitignore
+!setup_dependencies.sh

acc/mhd_solver/stencil_process.sps

@@ -202,9 +202,6 @@ heat_transfer(in Vector uu, in Scalar lnrho, in Scalar tt)
 
 
 #if LFORCING
-
-
-
 Vector
 simple_vortex_forcing(Vector a, Vector b, Scalar magnitude)
 {
@@ -222,7 +219,13 @@ simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude)
 Vector
 helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi)
 {
-
+    // JP: This looks wrong:
+    //    1) Should it be dsx * nx instead of dsx * ny?
+    //    2) Should you also use globalGrid.n instead of the local n?
+    //    3) Also final point: can we do this with vectors/quaternions instead?
+    //       Tringonometric functions are much more expensive and inaccurate/
+    // Fun related article:
+    // https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/
     xx.x = xx.x*(2.0*M_PI/(dsx*(DCONST_INT(AC_ny_max) - DCONST_INT(AC_ny_min))));
     xx.y = xx.y*(2.0*M_PI/(dsy*(DCONST_INT(AC_ny_max) - DCONST_INT(AC_ny_min))));
     xx.z = xx.z*(2.0*M_PI/(dsz*(DCONST_INT(AC_ny_max) - DCONST_INT(AC_ny_min))));

src/core/kernels/integration.cuh

@@ -59,17 +59,18 @@ create_rotz(const AcReal radians)
 }
 
 #if AC_DOUBLE_PRECISION == 0
+/*
+// Fast but inaccurate
 #define sin __sinf
 #define cos __cosf
 #define exp __expf
+*/
+#define sin sinf
+#define cos cosf
+#define exp expf
 #define rsqrt rsqrtf // hardware reciprocal sqrt
 #endif               // AC_DOUBLE_PRECISION == 0
 
-/*
-typedef struct {
-    int i, j, k;
-} int3;*/
-
 /*
  * =============================================================================
  * Level 0 (Input Assembly Stage)

src/standalone/autotest.cc

@@ -422,6 +422,12 @@ check_rk3(const AcMeshInfo& mesh_info)
             // const AcReal dt   = host_timestep(umax, mesh_info);
             const AcReal dt = AcReal(1e-2); // Use a small constant timestep to avoid instabilities
 
+#if LFORCING
+            const ForcingParams forcing_params = generateForcingParams(model_mesh->info);
+            loadForcingParamsToHost(forcing_params, model_mesh);
+            loadForcingParamsToDevice(forcing_params);
+#endif
+
             acIntegrate(dt);
 
             model_rk3(dt, model_mesh);

src/standalone/model/host_forcing.cc

@@ -255,6 +255,30 @@ loadForcingParamsToHost(const ForcingParams& forcing_params, AcMesh* mesh)
     mesh->info.real_params[AC_kaver] = forcing_params.kaver;
 }
 
+void
+loadForcingParamsToHost(const ForcingParams& forcing_params, ModelMesh* mesh)
+{
+    // %JP: Left some regex magic here in case we need to modify the ForcingParams struct
+    // acLoadDeviceConstant\(([A-Za-z_]*), ([a-z_.]*)\);
+    // mesh->info.real_params[$1] = $2;
+    mesh->info.real_params[AC_forcing_magnitude] = forcing_params.magnitude;
+    mesh->info.real_params[AC_forcing_phase]     = forcing_params.phase;
+
+    mesh->info.real_params[AC_k_forcex] = forcing_params.k_force.x;
+    mesh->info.real_params[AC_k_forcey] = forcing_params.k_force.y;
+    mesh->info.real_params[AC_k_forcez] = forcing_params.k_force.z;
+
+    mesh->info.real_params[AC_ff_hel_rex] = forcing_params.ff_hel_re.x;
+    mesh->info.real_params[AC_ff_hel_rey] = forcing_params.ff_hel_re.y;
+    mesh->info.real_params[AC_ff_hel_rez] = forcing_params.ff_hel_re.z;
+
+    mesh->info.real_params[AC_ff_hel_imx] = forcing_params.ff_hel_im.x;
+    mesh->info.real_params[AC_ff_hel_imy] = forcing_params.ff_hel_im.y;
+    mesh->info.real_params[AC_ff_hel_imz] = forcing_params.ff_hel_im.z;
+
+    mesh->info.real_params[AC_kaver] = forcing_params.kaver;
+}
+
 ForcingParams
 generateForcingParams(const AcMeshInfo& mesh_info)
 {

src/standalone/model/host_forcing.h

@@ -28,6 +28,8 @@
 #pragma once
 #include "astaroth.h"
 
+#include "modelmesh.h"
+
 AcReal get_random_number_01();
 
 AcReal3 cross(const AcReal3& a, const AcReal3& b);
@@ -64,5 +66,6 @@ typedef struct {
 void loadForcingParamsToDevice(const ForcingParams& forcing_params);
 
 void loadForcingParamsToHost(const ForcingParams& forcing_params, AcMesh* mesh);
+void loadForcingParamsToHost(const ForcingParams& forcing_params, ModelMesh* mesh);
 
 ForcingParams generateForcingParams(const AcMeshInfo& mesh_info);

src/standalone/model/model_rk3.cc

@@ -560,11 +560,9 @@ momentum(const ModelVectorData& uu, const ModelScalarData& lnrho
 #else
     // !!!!!!!!!!!!!!!!%JP: NOTE TODO IMPORTANT!!!!!!!!!!!!!!!!!!!!!!!!
     // NOT CHECKED FOR CORRECTNESS: USE AT YOUR OWN RISK
-    const ModelMatrix S   = stress_tensor(uu);
+    const ModelMatrix S = stress_tensor(uu);
-    const ModelScalar cs2 = get(AC_cs2_sound) *
-                            expl((get(AC_gamma) - 1) * (value(lnrho) - get(AC_lnrho0)));
 
-    const ModelVector mom = -mul(gradients(uu), value(uu)) - cs2 * gradient(lnrho) +
+    const ModelVector mom = -mul(gradients(uu), value(uu)) - get(AC_cs2_sound) * gradient(lnrho) +
                             get(AC_nu_visc) * (laplace_vec(uu) +
                                                ModelScalar(1. / 3.) * gradient_of_divergence(uu) +
                                                ModelScalar(2.) * mul(S, gradient(lnrho))) +
@@ -662,15 +660,13 @@ is_valid(const ModelVector& a)
 }
 
 #if LFORCING
-// FORCING NOT SUPPORTED FOR AUTOTEST
+ModelVector
-
-static inline ModelVector
 simple_vortex_forcing(ModelVector a, ModelVector b, ModelScalar magnitude)
 {
     return magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex
 }
 
-static inline ModelVector
+ModelVector
 simple_outward_flow_forcing(ModelVector a, ModelVector b, ModelScalar magnitude)
 {
     return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow
@@ -678,22 +674,24 @@ simple_outward_flow_forcing(ModelVector a, ModelVector b, ModelScalar magnitude)
 
 // The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable
 // helicity
-static inline ModelVector
+ModelVector
-helical_forcing(ModelScalar /* magnitude */, ModelVector k_force, ModelVector xx, ModelVector ff_re,
+helical_forcing(ModelScalar magnitude, ModelVector k_force, ModelVector xx, ModelVector ff_re,
                 ModelVector ff_im, ModelScalar phi)
 {
+    (void)magnitude; // WARNING: unused
+    xx.x = xx.x * (2.0l * M_PI / (get(AC_dsx) * (get(AC_ny_max) - get(AC_ny_min))));
+    xx.y = xx.y * (2.0l * M_PI / (get(AC_dsy) * (get(AC_ny_max) - get(AC_ny_min))));
+    xx.z = xx.z * (2.0l * M_PI / (get(AC_dsz) * (get(AC_ny_max) - get(AC_ny_min))));
 
-    xx.x = xx.x * (2.0 * M_PI / (get(AC_dsx) * (get(AC_ny_max) - get(AC_ny_min))));
+    ModelScalar cosl_phi     = cosl(phi);
-    xx.y = xx.y * (2.0 * M_PI / (get(AC_dsy) * (get(AC_ny_max) - get(AC_ny_min))));
+    ModelScalar sinl_phi     = sinl(phi);
-    xx.z = xx.z * (2.0 * M_PI / (get(AC_dsz) * (get(AC_ny_max) - get(AC_ny_min))));
+    ModelScalar cosl_k_dox_x = cosl(dot(k_force, xx));
-
+    ModelScalar sinl_k_dox_x = sinl(dot(k_force, xx));
-    ModelScalar cos_phi     = cosl(phi);
-    ModelScalar sin_phi     = sinl(phi);
-    ModelScalar cos_k_dox_x = cosl(dot(k_force, xx));
-    ModelScalar sin_k_dox_x = sinl(dot(k_force, xx));
     // Phase affect only the x-component
-    ModelScalar real_comp_phase = cos_k_dox_x * cos_phi - sin_k_dox_x * sin_phi;
+    // ModelScalar real_comp       = cosl_k_dox_x;
-    ModelScalar imag_comp_phase = cos_k_dox_x * sin_phi + sin_k_dox_x * cos_phi;
+    // ModelScalar imag_comp       = sinl_k_dox_x;
+    ModelScalar real_comp_phase = cosl_k_dox_x * cosl_phi - sinl_k_dox_x * sinl_phi;
+    ModelScalar imag_comp_phase = cosl_k_dox_x * sinl_phi + sinl_k_dox_x * cosl_phi;
 
     ModelVector force = (ModelVector){ff_re.x * real_comp_phase - ff_im.x * imag_comp_phase,
                                       ff_re.y * real_comp_phase - ff_im.y * imag_comp_phase,
@@ -702,18 +700,17 @@ helical_forcing(ModelScalar /* magnitude */, ModelVector k_force, ModelVector xx
     return force;
 }
 
-static inline ModelVector
+ModelVector
 forcing(int3 globalVertexIdx, ModelScalar dt)
 {
-    /*
+    ModelVector a = ModelScalar(.5) * (ModelVector){get(AC_nx) * get(AC_dsx),
-    ModelVector a  = ModelScalar(.5) * (ModelVector){get(AC_nx) * get(AC_dsx),
                                                     get(AC_ny) * get(AC_dsy),
                                                     get(AC_nz) * get(AC_dsz)}; // source (origin)
-                                                    */
+    (void)a;                                                                   // WARNING: not used
-    ModelVector xx = (ModelVector){
+    ModelVector xx        = (ModelVector){(globalVertexIdx.x - get(AC_nx_min)) * get(AC_dsx),
-        (globalVertexIdx.x - get(AC_nx_min)) * get(AC_dsx),
+                                   (globalVertexIdx.y - get(AC_ny_min)) * get(AC_dsy),
-        (globalVertexIdx.y - get(AC_ny_min) * get(AC_dsy)),
+                                   (globalVertexIdx.z - get(AC_nz_min)) *
-        (globalVertexIdx.z - get(AC_nz_min) * get(AC_dsz))}; // sink (current index)
+                                       get(AC_dsz)}; // sink (current index)
     const ModelScalar cs2 = get(AC_cs2_sound);
     const ModelScalar cs  = sqrtl(cs2);
 
@@ -724,6 +721,11 @@ forcing(int3 globalVertexIdx, ModelScalar dt)
     ModelVector ff_re = (ModelVector){get(AC_ff_hel_rex), get(AC_ff_hel_rey), get(AC_ff_hel_rez)};
     ModelVector ff_im = (ModelVector){get(AC_ff_hel_imx), get(AC_ff_hel_imy), get(AC_ff_hel_imz)};
 
+    (void)phase;   // WARNING: unused with simple forcing. Should be defined in helical_forcing
+    (void)k_force; // WARNING: unused with simple forcing. Should be defined in helical_forcing
+    (void)ff_re;   // WARNING: unused with simple forcing. Should be defined in helical_forcing
+    (void)ff_im;   // WARNING: unused with simple forcing. Should be defined in helical_forcing
+
     // Determine that forcing funtion type at this point.
     // ModelVector force = simple_vortex_forcing(a, xx, magnitude);
     // ModelVector force = simple_outward_flow_forcing(a, xx, magnitude);

src/standalone/renderer.cc

@@ -32,12 +32,13 @@
 #include <string.h> // memcpy
 
 #include "config_loader.h"
+#include "src/core/errchk.h"
+#include "src/core/math_utils.h"
+#include "model/host_forcing.h"
 #include "model/host_memory.h"
 #include "model/host_timestep.h"
 #include "model/model_reduce.h"
 #include "model/model_rk3.h"
-#include "src/core/errchk.h"
-#include "src/core/math_utils.h"
 #include "timer_hires.h"
 
 // Window
@@ -384,6 +385,12 @@ run_renderer(void)
 #if 1
         const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
         const AcReal dt   = host_timestep(umax, mesh_info);
+
+#if LFORCING
+        const ForcingParams forcing_params = generateForcingParams(mesh->info);
+        loadForcingParamsToDevice(forcing_params);
+#endif
+
         acIntegrate(dt);
 #else
         ModelMesh* model_mesh = modelmesh_create(mesh->info);