From 0e0ace397056a60af507de13ac3a468c7a823a27 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Tue, 6 Aug 2019 18:07:29 +0300 Subject: [PATCH 1/5] Pure hydro now works with autotests --- src/standalone/model/model_rk3.cc | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc index 3907f93..e033aa1 100644 --- a/src/standalone/model/model_rk3.cc +++ b/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 ModelScalar cs2 = get(AC_cs2_sound) * - expl((get(AC_gamma) - 1) * (value(lnrho) - get(AC_lnrho0))); + const ModelMatrix S = stress_tensor(uu); - 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))) + From 0b7f43da913c2f4754602f3efc79e6d478b8f582 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Tue, 6 Aug 2019 19:13:06 +0300 Subject: [PATCH 2/5] Updated 3rdparty .gitignore --- 3rdparty/.gitignore | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 3rdparty/.gitignore diff --git a/3rdparty/.gitignore b/3rdparty/.gitignore new file mode 100644 index 0000000..5e50e93 --- /dev/null +++ b/3rdparty/.gitignore @@ -0,0 +1,6 @@ +# Ignore everything +* + +# Except: +!.gitignore +!setup_dependencies.sh From d7e26e8f2141c39653a40499069816b61b50bccc Mon Sep 17 00:00:00 2001 From: jpekkila Date: Tue, 6 Aug 2019 19:15:28 +0300 Subject: [PATCH 3/5] Added forcing from stencil_process.sps to autotests. 3 Tests fail. --- acc/mhd_solver/stencil_defines.h | 2 +- config/astaroth.conf | 2 +- src/standalone/autotest.cc | 6 ++++ src/standalone/model/host_forcing.cc | 24 +++++++++++++ src/standalone/model/host_forcing.h | 3 ++ src/standalone/model/model_rk3.cc | 52 +++++++++++++++------------- src/standalone/renderer.cc | 7 ++++ 7 files changed, 70 insertions(+), 26 deletions(-) diff --git a/acc/mhd_solver/stencil_defines.h b/acc/mhd_solver/stencil_defines.h index 0ab07c2..5124820 100644 --- a/acc/mhd_solver/stencil_defines.h +++ b/acc/mhd_solver/stencil_defines.h @@ -30,7 +30,7 @@ #define LMAGNETIC (1) #define LENTROPY (1) #define LTEMPERATURE (0) -#define LFORCING (0) +#define LFORCING (1) #define LUPWD (0) #define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter diff --git a/config/astaroth.conf b/config/astaroth.conf index 41b7e51..32f50a3 100644 --- a/config/astaroth.conf +++ b/config/astaroth.conf @@ -40,7 +40,7 @@ AC_chi = 0.0001 AC_relhel = 0.0 AC_forcing_magnitude = 1e-5 AC_kmin = 0.8 -AC_kmax = 1.2 +AC_kmax = 1.2 // Entropy diff --git a/src/standalone/autotest.cc b/src/standalone/autotest.cc index 370eef5..56a8b75 100644 --- a/src/standalone/autotest.cc +++ b/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); diff --git a/src/standalone/model/host_forcing.cc b/src/standalone/model/host_forcing.cc index 423bf19..4ed09dd 100644 --- a/src/standalone/model/host_forcing.cc +++ b/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) { diff --git a/src/standalone/model/host_forcing.h b/src/standalone/model/host_forcing.h index da07c4b..383c663 100644 --- a/src/standalone/model/host_forcing.h +++ b/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); diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc index e033aa1..bb6dbc2 100644 --- a/src/standalone/model/model_rk3.cc +++ b/src/standalone/model/model_rk3.cc @@ -660,15 +660,13 @@ is_valid(const ModelVector& a) } #if LFORCING -// FORCING NOT SUPPORTED FOR AUTOTEST - -static inline ModelVector +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 @@ -676,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 -helical_forcing(ModelScalar /* magnitude */, ModelVector k_force, ModelVector xx, ModelVector ff_re, +ModelVector +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)))); - xx.y = xx.y * (2.0 * M_PI / (get(AC_dsy) * (get(AC_ny_max) - get(AC_ny_min)))); - xx.z = xx.z * (2.0 * M_PI / (get(AC_dsz) * (get(AC_ny_max) - get(AC_ny_min)))); - - 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)); + ModelScalar cosl_phi = cosl(phi); + ModelScalar sinl_phi = sinl(phi); + ModelScalar cosl_k_dox_x = cosl(dot(k_force, xx)); + ModelScalar sinl_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 imag_comp_phase = cos_k_dox_x * sin_phi + sin_k_dox_x * cos_phi; + // ModelScalar real_comp = cosl_k_dox_x; + // 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, @@ -700,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) - */ - ModelVector xx = (ModelVector){ - (globalVertexIdx.x - get(AC_nx_min)) * get(AC_dsx), - (globalVertexIdx.y - get(AC_ny_min) * get(AC_dsy)), - (globalVertexIdx.z - get(AC_nz_min) * get(AC_dsz))}; // sink (current index) + (void)a; // WARNING: not used + ModelVector xx = (ModelVector){(globalVertexIdx.x - get(AC_nx_min)) * get(AC_dsx), + (globalVertexIdx.y - get(AC_ny_min)) * get(AC_dsy), + (globalVertexIdx.z - get(AC_nz_min)) * + get(AC_dsz)}; // sink (current index) const ModelScalar cs2 = get(AC_cs2_sound); const ModelScalar cs = sqrtl(cs2); @@ -722,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); diff --git a/src/standalone/renderer.cc b/src/standalone/renderer.cc index 1522cc5..8bda077 100644 --- a/src/standalone/renderer.cc +++ b/src/standalone/renderer.cc @@ -34,6 +34,7 @@ #include "config_loader.h" #include "core/errchk.h" #include "core/math_utils.h" +#include "model/host_forcing.h" #include "model/host_memory.h" #include "model/host_timestep.h" #include "model/model_reduce.h" @@ -383,6 +384,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); From 6b53eb31ef95a8b9612627344263ddc684da3124 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Tue, 6 Aug 2019 19:29:40 +0300 Subject: [PATCH 4/5] Errors with forcing now down from 3 to 1 after switching from fast & inaccurate trig functions to more accurate ones --- src/core/kernels/integration.cuh | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/core/kernels/integration.cuh b/src/core/kernels/integration.cuh index 3503372..fec3b95 100644 --- a/src/core/kernels/integration.cuh +++ b/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) From a6fca069a7650ec8c3d0c0f2e2dde2a794f1812e Mon Sep 17 00:00:00 2001 From: jpekkila Date: Tue, 6 Aug 2019 19:47:03 +0300 Subject: [PATCH 5/5] Added a comment about helical forcing --- acc/mhd_solver/stencil_process.sps | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index eb0fc0c..c3ade38 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/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))));