From 8191c47fa074a96a368e8d386732738b6ee4cc02 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Tue, 25 Jun 2019 19:04:53 +0800 Subject: [PATCH 01/32] Scetching the helical forcing. Not the idead form. Not yet tested. --- acc/mhd_solver/stencil_assembly.sas | 2 +- acc/mhd_solver/stencil_process.sps | 46 ++++++++++++++++++++++++----- 2 files changed, 40 insertions(+), 8 deletions(-) diff --git a/acc/mhd_solver/stencil_assembly.sas b/acc/mhd_solver/stencil_assembly.sas index 06b1063..e2f95f7 100644 --- a/acc/mhd_solver/stencil_assembly.sas +++ b/acc/mhd_solver/stencil_assembly.sas @@ -1,5 +1,5 @@ -#define LUPWD (0) +#define LUPWD (1) Preprocessed Scalar diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 40b28ed..9b6780d 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -2,8 +2,8 @@ #define LENTROPY (1) #define LTEMPERATURE (0) #define LGRAVITY (0) -#define LFORCING (0) -#define LUPWD (0) +#define LFORCING (1) +#define LUPWD (1) // Declare uniforms (i.e. device constants) @@ -203,8 +203,12 @@ heat_transfer(in Vector uu, in Scalar lnrho, in Scalar tt) } #endif + + #if LFORCING + + Vector simple_vortex_forcing(Vector a, Vector b, Scalar magnitude) { @@ -217,24 +221,52 @@ simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude) return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow } + +// The Pencil Code hel_vec(), manual Eq. 222, inspired forcing function with adjustable helicity +Vector +helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi) +{ + Scalar cos_phi = cos(phi); + Scalar sin_phi = sin(phi); + Scalar cos_k_dox_x = cos(dot(k_force, xx)); + Scalar sin_k_dox_x = sin(dot(k_force, xx)); + Scalar real_comp = cos_k_dox_x*cos_phi - sin_k_dox_x*sin_phi; + Scalar imag_comp = cos_k_dox_x*sin_phi + sin_k_dox_x*cos_phi; + + + Vector force = (Vector){ ff_re.x*real_comp - ff_im.x*imag_comp, + ff_re.y*real_comp - ff_im.y*imag_comp, + ff_re.z*real_comp - ff_im.z*imag_comp}; + + return force; +} + Vector forcing(int3 globalVertexIdx) { Vector a = Scalar(.5) * (Vector){globalGrid.n.x * dsx, globalGrid.n.y * dsy, globalGrid.n.z * dsz}; // source (origin) - Vector b = (Vector){(globalVertexIdx.x - nx_min) * dsx, + Vector xx = (Vector){(globalVertexIdx.x - nx_min) * dsx, (globalVertexIdx.y - ny_min) * dsy, (globalVertexIdx.z - nz_min) * dsz}; // sink (current index) + + //Placeholders until determined properly Scalar magnitude = 0.05; + Vector k_force = (Vector){2.0, 0.0, 0.0}; + Vector ff_re = (Vector){0.0, 0.5, 0.0}; + Vector ff_im = (Vector){0.0, 0.8666, 0.0}; + Scalar phase = Scalar(0.79); + //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); + //Vector force = simple_vortex_forcing(a, xx, magnitude); + //Vector force = simple_outward_flow_forcing(a, xx, magnitude); + Vector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phi); - if (is_valid(c)) { return c; } - else { return (Vector){0, 0, 0}; } + if (is_valid(force)) { return force; } + else { return (Vector){0, 0, 0}; } } #endif // LFORCING From 5cacda285076d3e41f59d83ceff42535b87a9eb1 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Wed, 26 Jun 2019 13:15:28 +0800 Subject: [PATCH 02/32] Helical forcing funtion works. But we will need a wavenumber generator to add stochasticity. --- acc/mhd_solver/stencil_process.sps | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 9b6780d..714e742 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -263,7 +263,7 @@ forcing(int3 globalVertexIdx) //Determine that forcing funtion type at this point. //Vector force = simple_vortex_forcing(a, xx, magnitude); //Vector force = simple_outward_flow_forcing(a, xx, magnitude); - Vector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phi); + Vector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase); if (is_valid(force)) { return force; } else { return (Vector){0, 0, 0}; } From 231a8aa06ed5421c7ee66ba37ce9fb4282e0bf83 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Wed, 26 Jun 2019 15:23:46 +0800 Subject: [PATCH 03/32] Trying to figure out how to upload values to GPU. --- acc/mhd_solver/stencil_process.sps | 10 +++++----- src/core/astaroth.cu | 27 +++++++++++++++++++++++++++ src/standalone/simulation.cc | 14 ++++++++++++++ 3 files changed, 46 insertions(+), 5 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 714e742..58cfa96 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -253,11 +253,11 @@ forcing(int3 globalVertexIdx) //Placeholders until determined properly - Scalar magnitude = 0.05; - Vector k_force = (Vector){2.0, 0.0, 0.0}; - Vector ff_re = (Vector){0.0, 0.5, 0.0}; - Vector ff_im = (Vector){0.0, 0.8666, 0.0}; - Scalar phase = Scalar(0.79); + //Scalar magnitude = 0.05; + //Vector k_force = (Vector){2.0, 0.0, 0.0}; + //Vector ff_re = (Vector){0.0, 0.5, 0.0}; + //Vector ff_im = (Vector){0.0, 0.8666, 0.0}; + //Scalar phase = Scalar(0.79); //Determine that forcing funtion type at this point. diff --git a/src/core/astaroth.cu b/src/core/astaroth.cu index 956a24f..369082f 100644 --- a/src/core/astaroth.cu +++ b/src/core/astaroth.cu @@ -493,3 +493,30 @@ acSynchronize(void) return AC_SUCCESS; } + +//Tool for loadin forcing vector information into the device memory +AcResult +acForcingVec(const AcReal magnitude, const AcReal3 k_force, const AcReal3 ff_hel_re, + const AcReal3 ff_hel_im, const AcReal phase) +{ + + loadDeviceConstant(device[i], AC_forcing_magnitude, AC_forcing_magnitude); + loadDeviceConstant(device[i], AC_forcing_phase, AC_forcing_phase ); + + loadDeviceConstant(device[i], AC_k_forcex, k_force.x); + loadDeviceConstant(device[i], AC_k_forcey, k_force.y); + loadDeviceConstant(device[i], AC_k_forcez, k_force.z); + + + loadDeviceConstant(device[i], AC_ff_hel_rex, ff_hel_re.x); + loadDeviceConstant(device[i], AC_ff_hel_rey, ff_hel_re.y); + loadDeviceConstant(device[i], AC_ff_hel_rez, ff_hel_re.z); + + loadDeviceConstant(device[i], AC_ff_hel_imx, ff_hel_im.x); + loadDeviceConstant(device[i], AC_ff_hel_imy, ff_hel_im.y); + loadDeviceConstant(device[i], AC_ff_hel_imz, ff_hel_im.z); + + return AC_SUCCESS; +} + + diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index 57cd0fe..d4cf7b9 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -236,10 +236,24 @@ run_simulation(void) AcReal bin_save_t = mesh_info.real_params[AC_bin_save_t]; AcReal bin_crit_t = bin_save_t; + + /* Step the simulation */ for (int i = 1; i < max_steps; ++i) { const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); const AcReal dt = host_timestep(umax, mesh_info); + +#if LFORCING + //Generate a forcing vectors before canculating an integration step. + //Placeholders until determined properly + AcReal magnitude = 0.05; + AcReal phase = Scalar(0.79); + AcReal3 k_force = (AcReal3){2.0, 0.0, 0.0}; + AcReal3 ff_hel_re = (AcReal3){0.0, 0.5, 0.0}; + AcReal3 ff_hel_im = (AcReal3){0.0, 0.8666, 0.0}; + acForcingVec(magnitude, k_force, ff_hel_re,ff_hel_im, phase); +#endif + acIntegrate(dt); t_step += dt; From be0e46c814a52a144655951664a9af5bdd223fce Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Wed, 26 Jun 2019 17:41:39 +0800 Subject: [PATCH 04/32] Can move forcing vector information now from the host to device. next step in to generate random waves in the CPU with a chosen degree of helicity etc. --- acc/mhd_solver/stencil_process.sps | 10 ++++----- include/astaroth.h | 18 ++++++++++++++++ src/core/astaroth.cu | 33 +++++++++++++++--------------- src/standalone/simulation.cc | 2 +- 4 files changed, 41 insertions(+), 22 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 58cfa96..e4adae7 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -253,11 +253,11 @@ forcing(int3 globalVertexIdx) //Placeholders until determined properly - //Scalar magnitude = 0.05; - //Vector k_force = (Vector){2.0, 0.0, 0.0}; - //Vector ff_re = (Vector){0.0, 0.5, 0.0}; - //Vector ff_im = (Vector){0.0, 0.8666, 0.0}; - //Scalar phase = Scalar(0.79); + Scalar magnitude = DCONST_REAL(AC_forcing_magnitude); + Scalar phase = DCONST_REAL(AC_forcing_phase); + Vector k_force = (Vector){ DCONST_REAL(AC_k_forcex), DCONST_REAL(AC_k_forcey), DCONST_REAL(AC_k_forcez)}; + Vector ff_re = (Vector){DCONST_REAL(AC_ff_hel_rex), DCONST_REAL(AC_ff_hel_rey), DCONST_REAL(AC_ff_hel_rez)}; + Vector ff_im = (Vector){DCONST_REAL(AC_ff_hel_imx), DCONST_REAL(AC_ff_hel_imy), DCONST_REAL(AC_ff_hel_imz)}; //Determine that forcing funtion type at this point. diff --git a/include/astaroth.h b/include/astaroth.h index e6114c7..d9a347f 100644 --- a/include/astaroth.h +++ b/include/astaroth.h @@ -158,6 +158,19 @@ extern "C" { FUNC(AC_angl_uu), \ FUNC(AC_lnrho_edge),\ FUNC(AC_lnrho_out),\ + /* Forcing parameters. User configured. */\ + FUNC(AC_forcing_magnitude),\ + /* Forcing parameters. Set by the generator. */\ + FUNC(AC_forcing_phase),\ + FUNC(AC_k_forcex),\ + FUNC(AC_k_forcey),\ + FUNC(AC_k_forcez),\ + FUNC(AC_ff_hel_rex),\ + FUNC(AC_ff_hel_rey),\ + FUNC(AC_ff_hel_rez),\ + FUNC(AC_ff_hel_imx),\ + FUNC(AC_ff_hel_imy),\ + FUNC(AC_ff_hel_imz),\ /* Additional helper params */\ /* (deduced from other params do not set these directly!) */\ FUNC(AC_G_CONST),\ @@ -373,6 +386,11 @@ AcResult acQuit(void); unless otherwise stated. */ 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); + /* End extern "C" */ #ifdef __cplusplus } diff --git a/src/core/astaroth.cu b/src/core/astaroth.cu index 369082f..4ea78be 100644 --- a/src/core/astaroth.cu +++ b/src/core/astaroth.cu @@ -494,28 +494,29 @@ acSynchronize(void) return AC_SUCCESS; } -//Tool for loadin forcing vector information into the device memory +//Tool for loading forcing vector information into the device memory AcResult -acForcingVec(const AcReal magnitude, const AcReal3 k_force, const AcReal3 ff_hel_re, - const AcReal3 ff_hel_im, const AcReal phase) +acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force, const AcReal3 ff_hel_re, + const AcReal3 ff_hel_im, const AcReal forcing_phase) { - loadDeviceConstant(device[i], AC_forcing_magnitude, AC_forcing_magnitude); - loadDeviceConstant(device[i], AC_forcing_phase, AC_forcing_phase ); + for (int i = 0; i < num_devices; ++i) { + loadDeviceConstant(devices[i], AC_forcing_magnitude, forcing_magnitude); + loadDeviceConstant(devices[i], AC_forcing_phase, forcing_phase ); - loadDeviceConstant(device[i], AC_k_forcex, k_force.x); - loadDeviceConstant(device[i], AC_k_forcey, k_force.y); - loadDeviceConstant(device[i], AC_k_forcez, k_force.z); + loadDeviceConstant(devices[i], AC_k_forcex, k_force.x); + loadDeviceConstant(devices[i], AC_k_forcey, k_force.y); + loadDeviceConstant(devices[i], AC_k_forcez, k_force.z); + loadDeviceConstant(devices[i], AC_ff_hel_rex, ff_hel_re.x); + loadDeviceConstant(devices[i], AC_ff_hel_rey, ff_hel_re.y); + loadDeviceConstant(devices[i], AC_ff_hel_rez, ff_hel_re.z); - loadDeviceConstant(device[i], AC_ff_hel_rex, ff_hel_re.x); - loadDeviceConstant(device[i], AC_ff_hel_rey, ff_hel_re.y); - loadDeviceConstant(device[i], AC_ff_hel_rez, ff_hel_re.z); - - loadDeviceConstant(device[i], AC_ff_hel_imx, ff_hel_im.x); - loadDeviceConstant(device[i], AC_ff_hel_imy, ff_hel_im.y); - loadDeviceConstant(device[i], AC_ff_hel_imz, ff_hel_im.z); - + 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); + } + return AC_SUCCESS; } diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index d4cf7b9..51aff73 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -247,7 +247,7 @@ run_simulation(void) //Generate a forcing vectors before canculating an integration step. //Placeholders until determined properly AcReal magnitude = 0.05; - AcReal phase = Scalar(0.79); + AcReal phase = 0.79; AcReal3 k_force = (AcReal3){2.0, 0.0, 0.0}; AcReal3 ff_hel_re = (AcReal3){0.0, 0.5, 0.0}; AcReal3 ff_hel_im = (AcReal3){0.0, 0.8666, 0.0}; From 76d251cd3eb251147ac0f3a8dce62d8f9bb9d7b6 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Wed, 26 Jun 2019 18:50:42 +0800 Subject: [PATCH 05/32] Makes now special helical forcing vector. --- src/standalone/simulation.cc | 46 ++++++++++++++++++++++++++++++++++-- 1 file changed, 44 insertions(+), 2 deletions(-) diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index 51aff73..fd594eb 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -60,6 +60,44 @@ print_diagnostics(const AcMesh& mesh, const int& step, const AcReal& dt) } */ +//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) +{ + // k dot e + AcReal3 kdote; + kdote.x = k_force.x * e_force.x; + kdote.y = k_force.y * e_force.y; + kdote.z = k_force.z * e_force.z; + + // k cross e + AcReal3 k_cross_e; + k_cross_e.x=k_force.y*e_force.z-k_force.z*e_force.y; + k_cross_e.y=k_force.z*e_force.x-k_force.x*e_force.z; + k_cross_e.z=k_force.x*e_force.y-k_force.y*e_force.x; + + // k cross k cross e + AcReal3 k_cross_k_cross_e; + k_cross_k_cross_e.x=k_force.y*k_cross_e.z-k_force.z*k_cross_e.y; + k_cross_k_cross_e.y=k_force.z*k_cross_e.x-k_force.x*k_cross_e.z; + k_cross_k_cross_e.z=k_force.x*k_cross_e.y-k_force.y*k_cross_e.x; + + // abs(k) + AcReal kabs = sqrt(k_force.x*k_force.x + k_force.y*k_force.y + k_force.z*k_force.z); + + AcReal denominator = sqrt(AcReal(1.0) + relhel*relhel)*(kabs*kabs) + *sqrt(AcReal(1.0) - (kdote.x*kdote.x + kdote.y*kdote.y + kdote.z*kdote.z)/(kabs*kabs)); + + *ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator, + -relhel*kabs*k_cross_e.y/denominator, + -relhel*kabs*k_cross_e.z/denominator}; + + *ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator, + k_cross_k_cross_e.y/denominator, + k_cross_k_cross_e.z/denominator}; + +} + // Write all setting info into a separate ascii file. This is done to guarantee // that we have the data specifi information in the thing, even though in // principle these things are in the astaroth.conf. @@ -248,9 +286,13 @@ run_simulation(void) //Placeholders until determined properly AcReal magnitude = 0.05; AcReal phase = 0.79; + AcReal relhel = 0.5; AcReal3 k_force = (AcReal3){2.0, 0.0, 0.0}; - AcReal3 ff_hel_re = (AcReal3){0.0, 0.5, 0.0}; - AcReal3 ff_hel_im = (AcReal3){0.0, 0.8666, 0.0}; + AcReal3 e_force = (AcReal3){0.0, 2.0, 0.0}; + AcReal3 ff_hel_re, ff_hel_im; + helical_forcing_special_vector(&ff_hel_re, &ff_hel_im, k_force, e_force, relhel); + //AcReal3 ff_hel_re = (AcReal3){0.0, 0.5, 0.0}; + //AcReal3 ff_hel_im = (AcReal3){0.0, 0.8666, 0.0}; acForcingVec(magnitude, k_force, ff_hel_re,ff_hel_im, phase); #endif From 9ae3411ccecb4801dc80f3b92d287965c86d03fb Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Thu, 27 Jun 2019 14:53:36 +0800 Subject: [PATCH 06/32] helical_forcing_e_generator() added Without randomization. Will add next. --- src/standalone/simulation.cc | 69 +++++++++++++++++++++++++++++++++--- 1 file changed, 65 insertions(+), 4 deletions(-) diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index fd594eb..ecaacb7 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -60,6 +60,68 @@ print_diagnostics(const AcMesh& mesh, const int& step, const AcReal& dt) } */ +static inline AcReal3 +cross(const AcReal3& a, const AcReal3& b) +{ + AcReal3 c; + + c.x = a.y * b.z - a.z * b.y; + c.y = a.z * b.x - a.x * b.z; + c.z = a.x * b.y - a.y * b.x; + + return c; +} + +static inline AcReal +dot(const AcReal3& a, const AcReal3& b) +{ + return a.x * b.x + a.y * b.y + a.z * b.z; +} + +static inline AcReal3 +vec_norm(const AcReal3& a) +{ + AcReal3 c; + AcReal norm = dot(a, a); + + c.x = a.x/sqrt(norm); + c.y = a.y/sqrt(norm); + c.z = a.z/sqrt(norm); + + return c; +} + +static inline AcReal3 +vec_multi_scal(const AcReal scal, const AcReal3& a) +{ + AcReal3 c; + + c.x = a.x*scal; + c.y = a.y*scal; + c.z = a.z*scal; + + return c; +} + +//Generate the unit perpendicular unit vector e required for helical forcing +//Addapted from Pencil code forcing.f90 hel_vec() subroutine. +static inline void +helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force) +{ + + AcReal3 k_cross_e = cross(k_force, *e_force); + k_cross_e = vec_norm(k_cross_e); + AcReal3 k_cross_k_cross_e = cross(k_force, k_cross_e); + k_cross_k_cross_e = vec_norm(k_cross_k_cross_e); + AcReal phi = 2.9; //TODO RANDOMIZE [0, 2pi] + AcReal3 ee_tmp1 = vec_multi_scal(cos(phi),k_cross_e); + AcReal3 ee_tmp2 = vec_multi_scal(sin(phi), k_cross_k_cross_e); + + *e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, + ee_tmp1.y + ee_tmp2.y, + ee_tmp1.z + ee_tmp2.z}; +} + //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) @@ -287,12 +349,11 @@ run_simulation(void) AcReal magnitude = 0.05; AcReal phase = 0.79; AcReal relhel = 0.5; - AcReal3 k_force = (AcReal3){2.0, 0.0, 0.0}; - AcReal3 e_force = (AcReal3){0.0, 2.0, 0.0}; + AcReal3 k_force = (AcReal3){0.0, 2.0, 0.0}; + AcReal3 e_force = (AcReal3){0.0, 0.0, 1.0}; AcReal3 ff_hel_re, ff_hel_im; + helical_forcing_e_generator(&e_force, k_force); helical_forcing_special_vector(&ff_hel_re, &ff_hel_im, k_force, e_force, relhel); - //AcReal3 ff_hel_re = (AcReal3){0.0, 0.5, 0.0}; - //AcReal3 ff_hel_im = (AcReal3){0.0, 0.8666, 0.0}; acForcingVec(magnitude, k_force, ff_hel_re,ff_hel_im, phase); #endif From fd6a5df0d62c4c592d82e1713852056ed0111538 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Thu, 27 Jun 2019 15:59:58 +0800 Subject: [PATCH 07/32] helical_forcing_e_generator() randomized. --- src/standalone/simulation.cc | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index ecaacb7..7457d2b 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -60,6 +60,16 @@ print_diagnostics(const AcMesh& mesh, const int& step, const AcReal& dt) } */ +//The is a wrapper for genering random numbers with a chosen system. +static inline AcReal +get_random_number_01() +{ + //TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/ + return AcReal(rand())/AcReal(RAND_MAX); +} + + + static inline AcReal3 cross(const AcReal3& a, const AcReal3& b) { @@ -113,7 +123,7 @@ helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force) k_cross_e = vec_norm(k_cross_e); AcReal3 k_cross_k_cross_e = cross(k_force, k_cross_e); k_cross_k_cross_e = vec_norm(k_cross_k_cross_e); - AcReal phi = 2.9; //TODO RANDOMIZE [0, 2pi] + AcReal phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); AcReal3 ee_tmp1 = vec_multi_scal(cos(phi),k_cross_e); AcReal3 ee_tmp2 = vec_multi_scal(sin(phi), k_cross_k_cross_e); @@ -337,6 +347,9 @@ run_simulation(void) AcReal bin_crit_t = bin_save_t; + /* initialize random seed: */ + srand (312256655); + /* Step the simulation */ for (int i = 1; i < max_steps; ++i) { From 9b2e9d376f6b0a3542f80c66d6b01757f2c92a51 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Thu, 27 Jun 2019 18:12:15 +0800 Subject: [PATCH 08/32] helical_forcing_k_generator() added. Now Helical forcing almost works. I just need scale to force per tiome step correctly. The current formulation is wrong. --- src/standalone/simulation.cc | 46 ++++++++++++++++++++++++++++++++---- 1 file changed, 41 insertions(+), 5 deletions(-) diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index 7457d2b..7c9c209 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -113,6 +113,28 @@ vec_multi_scal(const AcReal scal, const AcReal3& a) return c; } +// Generate forcing wave vector k_force +static inline AcReal3 +helical_forcing_k_generator(const AcReal kmax, const AcReal kmin) +{ + AcReal phi, theta, kk; //Spherical direction coordinates + AcReal3 k_force; //forcing wave vector + + AcReal delta_k = kmax - kmin; + + // Generate vector in spherical coordinates + phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); + theta = AcReal(M_PI)*get_random_number_01(); + kk = delta_k*get_random_number_01() + kmin; + + // Cast into Cartesian form + k_force = (AcReal3){kk*sin(theta)*cos(phi), + kk*sin(theta)*sin(phi), + kk*cos(theta) }; + + return k_force; +} + //Generate the unit perpendicular unit vector e required for helical forcing //Addapted from Pencil code forcing.f90 hel_vec() subroutine. static inline void @@ -350,7 +372,6 @@ run_simulation(void) /* initialize random seed: */ srand (312256655); - /* Step the simulation */ for (int i = 1; i < max_steps; ++i) { const AcReal umax = acReduceVec(RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); @@ -360,12 +381,27 @@ run_simulation(void) //Generate a forcing vectors before canculating an integration step. //Placeholders until determined properly AcReal magnitude = 0.05; - AcReal phase = 0.79; AcReal relhel = 0.5; - AcReal3 k_force = (AcReal3){0.0, 2.0, 0.0}; - AcReal3 e_force = (AcReal3){0.0, 0.0, 1.0}; - AcReal3 ff_hel_re, ff_hel_im; + AcReal kmin = 1.3; + AcReal kmax = 1.7; + + // Generate forcing wave vector k_force + AcReal3 k_force;// = (AcReal3){0.0, 2.0, 0.0}; + k_force = helical_forcing_k_generator(kmax, kmin); + + //Randomize the phase + AcReal phase = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); + + // Generate e for k. Needed for the sake of isotrophy. + AcReal3 e_force; + if ((k_force.y == 0.0) && (k_force.z == 0.0)) { + e_force = (AcReal3){0.0, 1.0, 0.0}; + } else { + e_force = (AcReal3){1.0, 0.0, 0.0}; + } helical_forcing_e_generator(&e_force, k_force); + + 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); #endif From 94a25383a906d73bba175f251d0b1886af85f977 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Thu, 27 Jun 2019 19:20:18 +0800 Subject: [PATCH 09/32] Trying to calculate the forcing scaling. Causes nans very quickly. Will need to look closer tomorrow again. --- acc/mhd_solver/stencil_process.sps | 27 ++++++++++++++++++++++++--- src/standalone/simulation.cc | 2 +- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index e4adae7..c47481b 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/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 } diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index 7c9c209..725a606 100644 --- a/src/standalone/simulation.cc +++ b/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 From f04ef8e64cd0cd24e91c9d2a1c33fdf10ac7d77e Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Fri, 28 Jun 2019 19:23:18 +0800 Subject: [PATCH 10/32] Forcing function issue not yet fully resolved. Now brain hurs. No more today. Break needed. --- acc/mhd_solver/stencil_process.sps | 44 ++++++++++++++---------------- src/standalone/simulation.cc | 29 +++++++++++++------- 2 files changed, 39 insertions(+), 34 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index c47481b..947d532 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -222,32 +222,35 @@ simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude) } -// The Pencil Code hel_vec(), manual Eq. 222, inspired forcing function with adjustable helicity +// The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable helicity Vector helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi) { + + 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)))); + Scalar cos_phi = cos(phi); Scalar sin_phi = sin(phi); - Scalar cos_k_dox_x = cos(dot(k_force, xx)); - Scalar sin_k_dox_x = sin(dot(k_force, xx)); - Scalar real_comp = cos_k_dox_x*cos_phi - sin_k_dox_x*sin_phi; - Scalar imag_comp = cos_k_dox_x*sin_phi + sin_k_dox_x*cos_phi; + Scalar cos_k_dox_x = cos(dot(k_force, xx)); + Scalar sin_k_dox_x = sin(dot(k_force, xx)); + // Phase affect only the x-component + Scalar real_comp = cos_k_dox_x; + Scalar imag_comp = sin_k_dox_x; + Scalar real_comp_phase = cos_k_dox_x*cos_phi - sin_k_dox_x*sin_phi; + Scalar imag_comp_phase = cos_k_dox_x*sin_phi + sin_k_dox_x*cos_phi; - Vector force = (Vector){ ff_re.x*real_comp - ff_im.x*imag_comp, + Vector force = (Vector){ ff_re.x*real_comp_phase - ff_im.x*imag_comp_phase, ff_re.y*real_comp - ff_im.y*imag_comp, ff_re.z*real_comp - ff_im.z*imag_comp}; return force; } -#if LENTROPY -Vector -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, @@ -255,11 +258,7 @@ forcing(int3 globalVertexIdx, Scalar dt) 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 @@ -275,12 +274,13 @@ forcing(int3 globalVertexIdx, Scalar dt) //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) + //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((kk*cs)/dt); - force.x = NN*force.x; - force.y = NN*force.y; - force.z = NN*force.z; + 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; + force.z = force.z; if (is_valid(force)) { return force; } else { return (Vector){0, 0, 0}; } @@ -331,11 +331,7 @@ solve(Scalar dt) { #if LFORCING if (step_number == 2) { - #if LENTROPY - out_uu = out_uu + forcing(globalVertexIdx, dt, lnrho, ss); - #else out_uu = out_uu + forcing(globalVertexIdx, dt); - #endif } #endif } diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index 725a606..9ce5bd5 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -179,17 +179,26 @@ helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcR // abs(k) AcReal kabs = sqrt(k_force.x*k_force.x + k_force.y*k_force.y + k_force.z*k_force.z); - AcReal denominator = sqrt(AcReal(1.0) + relhel*relhel)*(kabs*kabs) - *sqrt(AcReal(1.0) - (kdote.x*kdote.x + kdote.y*kdote.y + kdote.z*kdote.z)/(kabs*kabs)); + AcReal denominator = sqrt(AcReal(1.0) + relhel*relhel)*kabs + *sqrt(kabs*kabs - (kdote.x*kdote.x + kdote.y*kdote.y + kdote.z*kdote.z)); - *ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator, - -relhel*kabs*k_cross_e.y/denominator, - -relhel*kabs*k_cross_e.z/denominator}; + //MV: I suspect there is a typo in the Pencil Code manual! + //*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator, + // -relhel*kabs*k_cross_e.y/denominator, + // -relhel*kabs*k_cross_e.z/denominator}; - *ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator, - k_cross_k_cross_e.y/denominator, - k_cross_k_cross_e.z/denominator}; + //*ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator, + // k_cross_k_cross_e.y/denominator, + // k_cross_k_cross_e.z/denominator}; + // See PC forcing.f90 rel_hel() + *ff_hel_re = (AcReal3){kabs*k_cross_e.x/denominator, + kabs*k_cross_e.y, + kabs*k_cross_e.z}; + + *ff_hel_im = (AcReal3){relhel*k_cross_k_cross_e.x/denominator, + relhel*k_cross_k_cross_e.y, + relhel*k_cross_k_cross_e.z}; } // Write all setting info into a separate ascii file. This is done to guarantee @@ -382,8 +391,8 @@ run_simulation(void) //Placeholders until determined properly AcReal magnitude = 0.05; AcReal relhel = 0.5; - AcReal kmin = 1.3; - AcReal kmax = 1.7; + AcReal kmin = 0.9999; + AcReal kmax = 1.0001; // Generate forcing wave vector k_force AcReal3 k_force; From 9f0be0d9ff333556c1c9c518bc3eab132d396472 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Mon, 1 Jul 2019 11:06:42 +0800 Subject: [PATCH 11/32] Solved the forcing function boundary problem. --- acc/mhd_solver/stencil_process.sps | 3 +- .../python/jupyter/notebook_example.ipynb | 14 +++------ include/astaroth.h | 5 ++-- src/core/astaroth.cu | 5 +++- src/standalone/simulation.cc | 30 ++++++++++++------- 5 files changed, 32 insertions(+), 25 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 947d532..f476921 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/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(kk*cs); // kk should be the average k!!! + const Scalar NN = cs*sqrt(DCONST_REAL(AC_kaver)*cs); //MV: Like in the Pencil Code. I don't understandf the logic here. force.x = sqrt(dt)*NN*force.x; force.y = force.y; diff --git a/analysis/python/jupyter/notebook_example.ipynb b/analysis/python/jupyter/notebook_example.ipynb index e85aef1..f7e4956 100644 --- a/analysis/python/jupyter/notebook_example.ipynb +++ b/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": { diff --git a/include/astaroth.h b/include/astaroth.h index 572d194..3776a60 100644 --- a/include/astaroth.h +++ b/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),\ @@ -408,8 +409,8 @@ 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); +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 AcReal kaver); /* End extern "C" */ #ifdef __cplusplus diff --git a/src/core/astaroth.cu b/src/core/astaroth.cu index d676613..007554e 100644 --- a/src/core/astaroth.cu +++ b/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; diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index 9ce5bd5..dc90dc0 100644 --- a/src/standalone/simulation.cc +++ b/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,13 +402,8 @@ run_simulation(void) const AcReal dt = host_timestep(umax, mesh_info); #if LFORCING - //Generate a forcing vectors 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 a forcing vector before canculating an integration step. + // Generate forcing wave vector k_force AcReal3 k_force; k_force = helical_forcing_k_generator(kmax, kmin); @@ -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); From 0600790f4128ba1db712bbf312576f30fdf4105d Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Mon, 1 Jul 2019 14:19:56 +0800 Subject: [PATCH 12/32] Corrected a bug in the timestep and some scaling problems. Now I can reach a saturated stated in forcing without crashing the code. --- acc/mhd_solver/stencil_process.sps | 12 ++++++------ src/standalone/model/host_timestep.cc | 3 ++- src/standalone/simulation.cc | 4 ++-- 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index f476921..39527e0 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -236,15 +236,15 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto Scalar cos_k_dox_x = cos(dot(k_force, xx)); Scalar sin_k_dox_x = sin(dot(k_force, xx)); // Phase affect only the x-component - Scalar real_comp = cos_k_dox_x; - Scalar imag_comp = sin_k_dox_x; + //Scalar real_comp = cos_k_dox_x; + //Scalar imag_comp = sin_k_dox_x; Scalar real_comp_phase = cos_k_dox_x*cos_phi - sin_k_dox_x*sin_phi; Scalar imag_comp_phase = cos_k_dox_x*sin_phi + sin_k_dox_x*cos_phi; Vector force = (Vector){ ff_re.x*real_comp_phase - ff_im.x*imag_comp_phase, - ff_re.y*real_comp - ff_im.y*imag_comp, - ff_re.z*real_comp - ff_im.z*imag_comp}; + ff_re.y*real_comp_phase - ff_im.y*imag_comp_phase, + ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase}; return force; } @@ -278,8 +278,8 @@ forcing(int3 globalVertexIdx, Scalar dt) const Scalar NN = cs*sqrt(DCONST_REAL(AC_kaver)*cs); //MV: Like in the Pencil Code. I don't understandf the logic here. force.x = sqrt(dt)*NN*force.x; - force.y = force.y; - force.z = force.z; + force.y = sqrt(dt)*NN*force.y; + force.z = sqrt(dt)*NN*force.z; if (is_valid(force)) { return force; } else { return (Vector){0, 0, 0}; } diff --git a/src/standalone/model/host_timestep.cc b/src/standalone/model/host_timestep.cc index 9c51d83..54f47cc 100644 --- a/src/standalone/model/host_timestep.cc +++ b/src/standalone/model/host_timestep.cc @@ -50,7 +50,8 @@ host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info) // New, closer to the actual Courant timestep // See Pencil Code user manual p. 38 (timestep section) const long double uu_dt = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + 0.0l)); - const long double visc_dt = cdtv * dsmin * dsmin / max(max(nu_visc, eta), max(gamma, chi)) + 1; // TODO NOTE: comment the +1 out to get scientifically accurate results + const long double visc_dt = cdtv * dsmin * dsmin / max(max(nu_visc, eta), max(gamma, chi));// + 1; // TODO NOTE: comment the +1 out to get scientifically accurate results + //MV: White the +1? It was messing up my computations! const long double dt = min(uu_dt, visc_dt); return AcReal(timescale) * AcReal(dt); diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index dc90dc0..659e382 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -386,8 +386,8 @@ run_simulation(void) AcReal bin_crit_t = bin_save_t; //Placeholders until determined properly - AcReal magnitude = 0.05; - AcReal relhel = 0.5; + AcReal magnitude = 1e-5; + AcReal relhel = 0.0; AcReal kmin = 0.8; AcReal kmax = 1.2; From d0eb308f17ca2d0ac5c82fa3aa77ca283d348ea6 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Tue, 2 Jul 2019 16:35:14 +0800 Subject: [PATCH 13/32] Better interface to forcing. --- config/astaroth.conf | 4 ++++ include/astaroth.h | 4 +++- src/standalone/simulation.cc | 11 ++++++----- 3 files changed, 13 insertions(+), 6 deletions(-) diff --git a/config/astaroth.conf b/config/astaroth.conf index 6d83764..41b7e51 100644 --- a/config/astaroth.conf +++ b/config/astaroth.conf @@ -38,6 +38,10 @@ AC_chi = 0.0001 // Forcing AC_relhel = 0.0 +AC_forcing_magnitude = 1e-5 +AC_kmin = 0.8 +AC_kmax = 1.2 + // Entropy AC_cp_sound = 1.0 diff --git a/include/astaroth.h b/include/astaroth.h index 3776a60..efc2ad1 100644 --- a/include/astaroth.h +++ b/include/astaroth.h @@ -142,7 +142,6 @@ extern "C" { FUNC(AC_cs_sound), \ FUNC(AC_eta), \ FUNC(AC_mu0), \ - FUNC(AC_relhel), \ FUNC(AC_cp_sound), \ FUNC(AC_gamma), \ FUNC(AC_cv_sound), \ @@ -160,6 +159,9 @@ extern "C" { FUNC(AC_lnrho_out),\ /* Forcing parameters. User configured. */\ FUNC(AC_forcing_magnitude),\ + FUNC(AC_relhel), \ + FUNC(AC_kmin), \ + FUNC(AC_kmax), \ /* Forcing parameters. Set by the generator. */\ FUNC(AC_forcing_phase),\ FUNC(AC_k_forcex),\ diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index 659e382..d40aeca 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -166,6 +166,7 @@ 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) { + // k dot e AcReal3 kdote; kdote.x = k_force.x * e_force.x; @@ -385,11 +386,11 @@ 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 = 1e-5; - AcReal relhel = 0.0; - AcReal kmin = 0.8; - AcReal kmax = 1.2; + //Forcing properties + AcReal relhel = mesh_info.real_params[AC_relhel]; + AcReal magnitude = mesh_info.real_params[AC_forcing_magnitude]; + AcReal kmin = mesh_info.real_params[AC_kmin]; + AcReal kmax = mesh_info.real_params[AC_kmax]; AcReal kaver = (kmax - kmin)/AcReal(2.0); From 4766441ffbbaa1aeb5d7cfb1cc5e7b1187c0f4c2 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Tue, 2 Jul 2019 18:24:41 +0800 Subject: [PATCH 14/32] Tryin to prepare autotest for forcing. --- acc/mhd_solver/stencil_process.sps | 2 +- src/standalone/autotest.cc | 2 + src/standalone/model/model_rk3.cc | 101 ++++++++++++++----- src/standalone/simulation.cc | 150 +---------------------------- 4 files changed, 81 insertions(+), 174 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 3547d5b..6065d24 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -3,7 +3,7 @@ #define LTEMPERATURE (0) #define LGRAVITY (0) #define LFORCING (1) -#define LUPWD (1) +#define LUPWD (0) // Declare uniforms (i.e. device constants) diff --git a/src/standalone/autotest.cc b/src/standalone/autotest.cc index 107b950..4385d5f 100644 --- a/src/standalone/autotest.cc +++ b/src/standalone/autotest.cc @@ -699,6 +699,8 @@ run_autotest(void) VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); const AcReal dt = host_timestep(umax, config); + + // Host integration step model_rk3(dt, model_mesh); boundconds(config, model_mesh); diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc index dd04bcf..afec5d0 100644 --- a/src/standalone/model/model_rk3.cc +++ b/src/standalone/model/model_rk3.cc @@ -709,27 +709,80 @@ is_valid(const ModelVector& a) return is_valid(a.x) && is_valid(a.y) && is_valid(a.z); } -static inline ModelVector -forcing(int3 globalVertexIdx) -{ - // source (origin) - ModelVector a = ModelScalar(.5) * (ModelVector){get(AC_nx) * get(AC_dsx), - get(AC_ny) * get(AC_dsy), - get(AC_nz) * get(AC_dsz)}; - // sink (current index) - ModelVector b = (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)}; - ModelScalar magnitude = 0.05; - // Vector c = magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow - ModelVector c = magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex - if (is_valid(c)) { - return c; - } - else { - return (ModelVector){0, 0, 0}; - } + +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 +simple_outward_flow_forcing(ModelVector a, ModelVector b, ModelScalar magnitude) +{ + return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow +} + + +// 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 ff_im, ModelScalar phi) +{ + + 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 = cos(phi); + ModelScalar sin_phi = sin(phi); + ModelScalar cos_k_dox_x = cos(dot(k_force, xx)); + ModelScalar sin_k_dox_x = sin(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; + + + 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, + ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase}; + + return force; +} + +static inline ModelVector +forcing(int3 globalVertexIdx, ModelScalar dt) +{ + 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) + const ModelScalar cs2 = get(AC_cs2_sound); + const ModelScalar cs = sqrt(cs2); + + //Placeholders until determined properly + ModelScalar magnitude = get(AC_forcing_magnitude); + ModelScalar phase = get(AC_forcing_phase); + ModelVector k_force = (ModelVector){ get(AC_k_forcex), get(AC_k_forcey), get(AC_k_forcez)}; + 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)}; + + + //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); + ModelVector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase); + + //Scaling N = magnitude*cs*sqrt(k*cs/dt) * dt + const ModelScalar NN = cs*sqrt(get(AC_kaver)*cs); + //MV: Like in the Pencil Code. I don't understandf the logic here. + force.x = sqrt(dt)*NN*force.x; + force.y = sqrt(dt)*NN*force.y; + force.z = sqrt(dt)*NN*force.z; + + if (is_valid(force)) { return force; } + else { return (ModelVector){0, 0, 0}; } } static void @@ -795,10 +848,10 @@ solve_beta_step(const int step_number, const ModelScalar dt, const int i, const #if LFORCING if (step_number == 2) { - ModelVector force = forcing((int3){i, j, k}); - out->vertex_buffer[VTXBUF_UUX][idx] += force.x * dt; - out->vertex_buffer[VTXBUF_UUY][idx] += force.y * dt; - out->vertex_buffer[VTXBUF_UUZ][idx] += force.z * dt; + ModelVector force = forcing((int3){i, j, k}, dt); + out->vertex_buffer[VTXBUF_UUX][idx] += force.x ; + out->vertex_buffer[VTXBUF_UUY][idx] += force.y ; + out->vertex_buffer[VTXBUF_UUZ][idx] += force.z ; } #endif } diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index d40aeca..ef13af6 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -33,6 +33,7 @@ #include "model/host_timestep.h" #include "model/model_reduce.h" #include "model/model_rk3.h" +#include "model/host_forcing.h" #include "timer_hires.h" #include @@ -60,155 +61,6 @@ print_diagnostics(const AcMesh& mesh, const int& step, const AcReal& dt) } */ -//The is a wrapper for genering random numbers with a chosen system. -static inline AcReal -get_random_number_01() -{ - //TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/ - return AcReal(rand())/AcReal(RAND_MAX); -} - - - -static inline AcReal3 -cross(const AcReal3& a, const AcReal3& b) -{ - AcReal3 c; - - c.x = a.y * b.z - a.z * b.y; - c.y = a.z * b.x - a.x * b.z; - c.z = a.x * b.y - a.y * b.x; - - return c; -} - -static inline AcReal -dot(const AcReal3& a, const AcReal3& b) -{ - return a.x * b.x + a.y * b.y + a.z * b.z; -} - -static inline AcReal3 -vec_norm(const AcReal3& a) -{ - AcReal3 c; - AcReal norm = dot(a, a); - - c.x = a.x/sqrt(norm); - c.y = a.y/sqrt(norm); - c.z = a.z/sqrt(norm); - - return c; -} - -static inline AcReal3 -vec_multi_scal(const AcReal scal, const AcReal3& a) -{ - AcReal3 c; - - c.x = a.x*scal; - c.y = a.y*scal; - c.z = a.z*scal; - - return c; -} - -// Generate forcing wave vector k_force -static inline AcReal3 -helical_forcing_k_generator(const AcReal kmax, const AcReal kmin) -{ - AcReal phi, theta, kk; //Spherical direction coordinates - AcReal3 k_force; //forcing wave vector - - AcReal delta_k = kmax - kmin; - - // Generate vector in spherical coordinates - phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); - theta = AcReal(M_PI)*get_random_number_01(); - kk = delta_k*get_random_number_01() + kmin; - - // Cast into Cartesian form - k_force = (AcReal3){kk*sin(theta)*cos(phi), - 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; -} - -//Generate the unit perpendicular unit vector e required for helical forcing -//Addapted from Pencil code forcing.f90 hel_vec() subroutine. -static inline void -helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force) -{ - - AcReal3 k_cross_e = cross(k_force, *e_force); - k_cross_e = vec_norm(k_cross_e); - AcReal3 k_cross_k_cross_e = cross(k_force, k_cross_e); - k_cross_k_cross_e = vec_norm(k_cross_k_cross_e); - AcReal phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); - AcReal3 ee_tmp1 = vec_multi_scal(cos(phi),k_cross_e); - AcReal3 ee_tmp2 = vec_multi_scal(sin(phi), k_cross_k_cross_e); - - *e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, - ee_tmp1.y + ee_tmp2.y, - ee_tmp1.z + ee_tmp2.z}; -} - -//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) -{ - - // k dot e - AcReal3 kdote; - kdote.x = k_force.x * e_force.x; - kdote.y = k_force.y * e_force.y; - kdote.z = k_force.z * e_force.z; - - // k cross e - AcReal3 k_cross_e; - k_cross_e.x=k_force.y*e_force.z-k_force.z*e_force.y; - k_cross_e.y=k_force.z*e_force.x-k_force.x*e_force.z; - k_cross_e.z=k_force.x*e_force.y-k_force.y*e_force.x; - - // k cross k cross e - AcReal3 k_cross_k_cross_e; - k_cross_k_cross_e.x=k_force.y*k_cross_e.z-k_force.z*k_cross_e.y; - k_cross_k_cross_e.y=k_force.z*k_cross_e.x-k_force.x*k_cross_e.z; - k_cross_k_cross_e.z=k_force.x*k_cross_e.y-k_force.y*k_cross_e.x; - - // abs(k) - AcReal kabs = sqrt(k_force.x*k_force.x + k_force.y*k_force.y + k_force.z*k_force.z); - - AcReal denominator = sqrt(AcReal(1.0) + relhel*relhel)*kabs - *sqrt(kabs*kabs - (kdote.x*kdote.x + kdote.y*kdote.y + kdote.z*kdote.z)); - - //MV: I suspect there is a typo in the Pencil Code manual! - //*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator, - // -relhel*kabs*k_cross_e.y/denominator, - // -relhel*kabs*k_cross_e.z/denominator}; - - //*ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator, - // k_cross_k_cross_e.y/denominator, - // k_cross_k_cross_e.z/denominator}; - - // 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}; - - *ff_hel_im = (AcReal3){relhel*k_cross_k_cross_e.x/denominator, - relhel*k_cross_k_cross_e.y, - relhel*k_cross_k_cross_e.z}; -} // Write all setting info into a separate ascii file. This is done to guarantee // that we have the data specifi information in the thing, even though in From 334ff868d95c304b433cdc3243a404776766676f Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Tue, 2 Jul 2019 18:46:04 +0800 Subject: [PATCH 15/32] Forcing disabled from autotest and from defaults. It is not suitable function of the autotest tool. If there in really a mandatory need to add it. I will need special help from Johannes. --- acc/mhd_solver/stencil_process.sps | 2 +- include/astaroth.h | 2 +- src/standalone/autotest.cc | 8 +++++++- src/standalone/model/model_rk3.cc | 5 ++++- 4 files changed, 13 insertions(+), 4 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 6065d24..4960554 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -2,7 +2,7 @@ #define LENTROPY (1) #define LTEMPERATURE (0) #define LGRAVITY (0) -#define LFORCING (1) +#define LFORCING (0) #define LUPWD (0) diff --git a/include/astaroth.h b/include/astaroth.h index 45159f5..504f768 100644 --- a/include/astaroth.h +++ b/include/astaroth.h @@ -76,7 +76,7 @@ extern "C" { #define LMAGNETIC (1) #define LENTROPY (1) #define LTEMPERATURE (0) - #define LFORCING (1) + #define LFORCING (0) #endif // clang-format on diff --git a/src/standalone/autotest.cc b/src/standalone/autotest.cc index 4385d5f..0f961fa 100644 --- a/src/standalone/autotest.cc +++ b/src/standalone/autotest.cc @@ -1,5 +1,5 @@ /* - Copyright (C) 2014-2018, Johannes Pekkilae, Miikka Vaeisalae. + Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae. This file is part of Astaroth. @@ -32,6 +32,7 @@ #include "core/math_utils.h" #include "model/host_memory.h" #include "model/host_timestep.h" +#include "model/host_forcing.h" #include "model/model_boundconds.h" #include "model/model_reduce.h" #include "model/model_rk3.h" @@ -699,6 +700,11 @@ run_autotest(void) VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); const AcReal dt = host_timestep(umax, config); +#if LFORCING + + //CURRENTLY AUTOTEST NOT SUPPORTED WITH FORCING!!! + +#endif // Host integration step diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc index afec5d0..6c79908 100644 --- a/src/standalone/model/model_rk3.cc +++ b/src/standalone/model/model_rk3.cc @@ -1,5 +1,5 @@ /* - Copyright (C) 2014-2018, Johannes Pekkilae, Miikka Vaeisalae. + Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae. This file is part of Astaroth. @@ -710,6 +710,8 @@ is_valid(const ModelVector& a) } +#if 0 +//FORCING NOT SUPPORTED FOR AUTOTEST static inline ModelVector simple_vortex_forcing(ModelVector a, ModelVector b, ModelScalar magnitude) @@ -784,6 +786,7 @@ forcing(int3 globalVertexIdx, ModelScalar dt) if (is_valid(force)) { return force; } else { return (ModelVector){0, 0, 0}; } } +#endif static void solve_alpha_step(const int step_number, const ModelScalar dt, const int i, const int j, const int k, From f0d2be831ecb322c341fe47f76de5294d88d1118 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Wed, 3 Jul 2019 09:55:23 +0800 Subject: [PATCH 16/32] host_forcing now committed. Sorry. --- src/standalone/model/host_forcing.cc | 180 +++++++++++++++++++++++++++ src/standalone/model/host_forcing.h | 45 +++++++ 2 files changed, 225 insertions(+) create mode 100644 src/standalone/model/host_forcing.cc create mode 100644 src/standalone/model/host_forcing.h diff --git a/src/standalone/model/host_forcing.cc b/src/standalone/model/host_forcing.cc new file mode 100644 index 0000000..ecac783 --- /dev/null +++ b/src/standalone/model/host_forcing.cc @@ -0,0 +1,180 @@ +/* + Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae. + + This file is part of Astaroth. + + Astaroth is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Astaroth is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Astaroth. If not, see . +*/ + +/** + * @file + * \brief Brief info. + * + * Detailed info. + * + */ +#include "host_forcing.h" + +#include "core/math_utils.h" + + +//The is a wrapper for genering random numbers with a chosen system. +AcReal +get_random_number_01() +{ + //TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/ + return AcReal(rand())/AcReal(RAND_MAX); +} + + + +AcReal3 +cross(const AcReal3& a, const AcReal3& b) +{ + AcReal3 c; + + c.x = a.y * b.z - a.z * b.y; + c.y = a.z * b.x - a.x * b.z; + c.z = a.x * b.y - a.y * b.x; + + return c; +} + +AcReal +dot(const AcReal3& a, const AcReal3& b) +{ + return a.x * b.x + a.y * b.y + a.z * b.z; +} + +AcReal3 +vec_norm(const AcReal3& a) +{ + AcReal3 c; + AcReal norm = dot(a, a); + + c.x = a.x/sqrt(norm); + c.y = a.y/sqrt(norm); + c.z = a.z/sqrt(norm); + + return c; +} + +AcReal3 +vec_multi_scal(const AcReal scal, const AcReal3& a) +{ + AcReal3 c; + + c.x = a.x*scal; + c.y = a.y*scal; + c.z = a.z*scal; + + return c; +} + +// Generate forcing wave vector k_force +AcReal3 +helical_forcing_k_generator(const AcReal kmax, const AcReal kmin) +{ + AcReal phi, theta, kk; //Spherical direction coordinates + AcReal3 k_force; //forcing wave vector + + AcReal delta_k = kmax - kmin; + + // Generate vector in spherical coordinates + phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); + theta = AcReal(M_PI)*get_random_number_01(); + kk = delta_k*get_random_number_01() + kmin; + + // Cast into Cartesian form + k_force = (AcReal3){kk*sin(theta)*cos(phi), + 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; +} + +//Generate the unit perpendicular unit vector e required for helical forcing +//Addapted from Pencil code forcing.f90 hel_vec() subroutine. +void +helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force) +{ + + AcReal3 k_cross_e = cross(k_force, *e_force); + k_cross_e = vec_norm(k_cross_e); + AcReal3 k_cross_k_cross_e = cross(k_force, k_cross_e); + k_cross_k_cross_e = vec_norm(k_cross_k_cross_e); + AcReal phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); + AcReal3 ee_tmp1 = vec_multi_scal(cos(phi),k_cross_e); + AcReal3 ee_tmp2 = vec_multi_scal(sin(phi), k_cross_k_cross_e); + + *e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, + ee_tmp1.y + ee_tmp2.y, + ee_tmp1.z + ee_tmp2.z}; +} + +//PC Manual Eq. 223 +void +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; + kdote.x = k_force.x * e_force.x; + kdote.y = k_force.y * e_force.y; + kdote.z = k_force.z * e_force.z; + + // k cross e + AcReal3 k_cross_e; + k_cross_e.x=k_force.y*e_force.z-k_force.z*e_force.y; + k_cross_e.y=k_force.z*e_force.x-k_force.x*e_force.z; + k_cross_e.z=k_force.x*e_force.y-k_force.y*e_force.x; + + // k cross k cross e + AcReal3 k_cross_k_cross_e; + k_cross_k_cross_e.x=k_force.y*k_cross_e.z-k_force.z*k_cross_e.y; + k_cross_k_cross_e.y=k_force.z*k_cross_e.x-k_force.x*k_cross_e.z; + k_cross_k_cross_e.z=k_force.x*k_cross_e.y-k_force.y*k_cross_e.x; + + // abs(k) + AcReal kabs = sqrt(k_force.x*k_force.x + k_force.y*k_force.y + k_force.z*k_force.z); + + AcReal denominator = sqrt(AcReal(1.0) + relhel*relhel)*kabs + *sqrt(kabs*kabs - (kdote.x*kdote.x + kdote.y*kdote.y + kdote.z*kdote.z)); + + //MV: I suspect there is a typo in the Pencil Code manual! + //*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator, + // -relhel*kabs*k_cross_e.y/denominator, + // -relhel*kabs*k_cross_e.z/denominator}; + + //*ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator, + // k_cross_k_cross_e.y/denominator, + // k_cross_k_cross_e.z/denominator}; + + // 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}; + + *ff_hel_im = (AcReal3){relhel*k_cross_k_cross_e.x/denominator, + relhel*k_cross_k_cross_e.y, + relhel*k_cross_k_cross_e.z}; +} diff --git a/src/standalone/model/host_forcing.h b/src/standalone/model/host_forcing.h new file mode 100644 index 0000000..569485e --- /dev/null +++ b/src/standalone/model/host_forcing.h @@ -0,0 +1,45 @@ + +/* + Copyright (C) 2014-2019, Johannes Pekkilae, Miikka Vaeisalae. + + This file is part of Astaroth. + + Astaroth is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Astaroth is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Astaroth. If not, see . +*/ + +/** + * @file + * \brief Brief info. + * + * Detailed info. + * + */ +#pragma once +#include "astaroth.h" + +AcReal get_random_number_01(); + +AcReal3 cross(const AcReal3& a, const AcReal3& b); + +AcReal dot(const AcReal3& a, const AcReal3& b); + +AcReal3 vec_norm(const AcReal3& a); + +AcReal3 vec_multi_scal(const AcReal scal, const AcReal3& a); + +AcReal3 helical_forcing_k_generator(const AcReal kmax, const AcReal kmin); + +void helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force); + +void helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force, const AcReal3 e_force, const AcReal relhel); From 98713ff9d2ade911472e8a24315eb9334312f1c8 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Wed, 3 Jul 2019 14:49:10 +0800 Subject: [PATCH 17/32] A possible bug note added. Will look into late. --- acc/mhd_solver/stencil_process.sps | 1 + 1 file changed, 1 insertion(+) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 4960554..8fbd8ae 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -69,6 +69,7 @@ momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) { const Scalar cs2 = cs2_sound * exp(gamma * value(ss) / cp_sound + (gamma - 1) * (value(lnrho) - LNRHO0)); const Vector j = (Scalar(1.) / mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density const Vector B = curl(aa); + //TODO: DOES INTHERMAL VERSTION INCLUDE THE MAGNETIC FIELD? const Scalar inv_rho = Scalar(1.) / exp(value(lnrho)); // Regex replace CPU constants with get\(AC_([a-zA-Z_0-9]*)\) From 46a2ef4847df197a2c4cb567239d5cade7d7b15d Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Wed, 3 Jul 2019 15:13:01 +0800 Subject: [PATCH 18/32] Commit demonstration for student. --- doc/testfile.txt | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 doc/testfile.txt diff --git a/doc/testfile.txt b/doc/testfile.txt new file mode 100644 index 0000000..829a508 --- /dev/null +++ b/doc/testfile.txt @@ -0,0 +1,3 @@ +Lorem ipsum + + From 32c85205438564b6e69c827f5b7008d0dddadc64 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Wed, 3 Jul 2019 17:37:27 +0800 Subject: [PATCH 19/32] Removed testfile.txt. Demonstration over. --- doc/testfile.txt | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 doc/testfile.txt diff --git a/doc/testfile.txt b/doc/testfile.txt deleted file mode 100644 index 829a508..0000000 --- a/doc/testfile.txt +++ /dev/null @@ -1,3 +0,0 @@ -Lorem ipsum - - From d4d2680f4058566e8235795262944530f612108b Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 16:19:25 +0300 Subject: [PATCH 20/32] Added a new generic function to the interface (astaroth.h) for loading arbitrary device constants. Also (unintended) autoformatting. --- include/astaroth.h | 10 +++++++--- src/core/astaroth.cu | 20 +++++++++++++------- 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/include/astaroth.h b/include/astaroth.h index 504f768..5d0e404 100644 --- a/include/astaroth.h +++ b/include/astaroth.h @@ -409,11 +409,15 @@ AcResult acQuit(void); unless otherwise stated. */ AcResult acSynchronize(void); +/** */ +AcResult acloadDeviceConstant(const AcRealParam param, const AcReal value); + /** 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 AcReal kaver); - +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 AcReal kaver); + /* End extern "C" */ #ifdef __cplusplus } diff --git a/src/core/astaroth.cu b/src/core/astaroth.cu index a0cf7ae..a912e90 100644 --- a/src/core/astaroth.cu +++ b/src/core/astaroth.cu @@ -506,15 +506,24 @@ acSynchronize(void) return AC_SUCCESS; } -//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, +acloadDeviceConstant(const AcRealParam param, const AcReal value) +{ + for (int i = 0; i < num_devices; ++i) { + loadDeviceConstant(devices[i], param, value); + } + return AC_SUCCESS; +} + +// 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 AcReal kaver) { for (int i = 0; i < num_devices; ++i) { loadDeviceConstant(devices[i], AC_forcing_magnitude, forcing_magnitude); - loadDeviceConstant(devices[i], AC_forcing_phase, forcing_phase ); + loadDeviceConstant(devices[i], AC_forcing_phase, forcing_phase); loadDeviceConstant(devices[i], AC_k_forcex, k_force.x); loadDeviceConstant(devices[i], AC_k_forcey, k_force.y); @@ -529,10 +538,7 @@ acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force, const AcReal loadDeviceConstant(devices[i], AC_ff_hel_imz, ff_hel_im.z); loadDeviceConstant(devices[i], AC_kaver, kaver); - } - + return AC_SUCCESS; } - - From 08e9a32cb14c7a4d98ed888a2df2f40cc9f0ccb1 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 16:37:16 +0300 Subject: [PATCH 21/32] Added a comment about acForcingVec --- include/astaroth.h | 2 +- src/core/astaroth.cu | 6 +++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/include/astaroth.h b/include/astaroth.h index 5d0e404..f5582dc 100644 --- a/include/astaroth.h +++ b/include/astaroth.h @@ -410,7 +410,7 @@ AcResult acQuit(void); AcResult acSynchronize(void); /** */ -AcResult acloadDeviceConstant(const AcRealParam param, const AcReal value); +AcResult acLoadDeviceConstant(const AcRealParam param, const AcReal value); /** Tool for loading forcing vector information into the device memory */ diff --git a/src/core/astaroth.cu b/src/core/astaroth.cu index a912e90..0c49477 100644 --- a/src/core/astaroth.cu +++ b/src/core/astaroth.cu @@ -507,7 +507,7 @@ acSynchronize(void) } AcResult -acloadDeviceConstant(const AcRealParam param, const AcReal value) +acLoadDeviceConstant(const AcRealParam param, const AcReal value) { for (int i = 0; i < num_devices; ++i) { loadDeviceConstant(devices[i], param, value); @@ -516,6 +516,10 @@ acloadDeviceConstant(const AcRealParam param, const AcReal value) } // Tool for loading forcing vector information into the device memory +// %JP: Added a generic function for loading device constants (acLoadDeviceConstant). +// This acForcingVec should go outside the core library since it references user-defined +// parameters such as AC_forcing_magnitude which may not be defined in all projects. +// host_forcing.cc is probably a good place for this. 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 AcReal kaver) From d4968d05832e14ccd231c5b0af4e16a2c02fe021 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 16:38:31 +0300 Subject: [PATCH 22/32] Made the gcc warning flags stricter --- CMakeLists.txt | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 349201e..c1f5824 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -49,7 +49,7 @@ message(STATUS "CMAKE_CXX_COMPILER: " ${CMAKE_CXX_COMPILER_ID}) if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 6.0) # Because of GCC bug 48891 - message(FATAL_ERROR "GCC version 6.0 or higher required") + message(FATAL_ERROR "GCC version 6.0 or higher required") endif() endif() @@ -128,8 +128,7 @@ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}\ set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG}\ -O0 -g") -set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror -Wno-error=unused-parameter\ - -Wno-error=unused-function -Wno-error=unknown-pragmas") +set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror") # Also warn about implicit conversions if the compiler supports it if (${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU") From 25d4b9a0cda84e236a35a58e94ffdacbb00def9a Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 16:54:51 +0300 Subject: [PATCH 23/32] Added compilation warning flags for the Intel compiler. --- CMakeLists.txt | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c1f5824..533f38a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -128,18 +128,14 @@ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}\ set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG}\ -O0 -g") -set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror") - -# Also warn about implicit conversions if the compiler supports it -if (${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU") - set (CXX_FLAGS_WARNING "${CXX_FLAGS_WARNING} -Wdouble-promotion -Wfloat-conversion") +if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") + set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror -Wdouble-promotion -Wfloat-conversion") +elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") + set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror-all -Wsign-conversion") +else() + message(WARNING "Using an unknown compiler. Compilation warning flags were not set.") endif() -# Other flags. -D_FORCE_INLINES is a workaround to some CUDA/C++ "feature" -# which botches the compilation ("memcpy was not declared in this scope") -# (Not required with cc >= 3.0) -#set(CXX_FLAGS_ETC "-D_FORCE_INLINES") - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}\ ${CXX_FLAGS_WARNING}\ ${CXX_FLAGS_ETC}\ From 59ac2647433947fb109b30c510f19f3e180efea0 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 16:57:57 +0300 Subject: [PATCH 24/32] simulation.cc autoformatting --- src/standalone/simulation.cc | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index ef13af6..76cc40f 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -29,11 +29,11 @@ #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" #include "model/model_rk3.h" -#include "model/host_forcing.h" #include "timer_hires.h" #include @@ -61,7 +61,6 @@ print_diagnostics(const AcMesh& mesh, const int& step, const AcReal& dt) } */ - // Write all setting info into a separate ascii file. This is done to guarantee // that we have the data specifi information in the thing, even though in // principle these things are in the astaroth.conf. @@ -204,7 +203,7 @@ run_simulation(void) load_config(&mesh_info); AcMesh* mesh = acmesh_create(mesh_info); - //TODO: This need to be possible to define in astaroth.conf + // TODO: This need to be possible to define in astaroth.conf acmesh_init_to(INIT_TYPE_GAUSSIAN_RADIAL_EXPL, mesh); acInit(mesh_info); @@ -231,23 +230,23 @@ run_simulation(void) acStore(mesh); save_mesh(*mesh, 0, t_step); - const int max_steps = mesh_info.int_params[AC_max_steps]; - const int save_steps = mesh_info.int_params[AC_save_steps]; + const int max_steps = mesh_info.int_params[AC_max_steps]; + const int save_steps = mesh_info.int_params[AC_save_steps]; const int bin_save_steps = mesh_info.int_params[AC_bin_steps]; // TODO Get from mesh_info AcReal bin_save_t = mesh_info.real_params[AC_bin_save_t]; AcReal bin_crit_t = bin_save_t; - //Forcing properties - AcReal relhel = mesh_info.real_params[AC_relhel]; + // Forcing properties + AcReal relhel = mesh_info.real_params[AC_relhel]; AcReal magnitude = mesh_info.real_params[AC_forcing_magnitude]; AcReal kmin = mesh_info.real_params[AC_kmin]; AcReal kmax = mesh_info.real_params[AC_kmax]; - AcReal kaver = (kmax - kmin)/AcReal(2.0); + AcReal kaver = (kmax - kmin) / AcReal(2.0); /* initialize random seed: */ - srand (312256655); + srand(312256655); /* Step the simulation */ for (int i = 1; i < max_steps; ++i) { @@ -255,20 +254,21 @@ run_simulation(void) const AcReal dt = host_timestep(umax, mesh_info); #if LFORCING - //Generate a forcing vector before canculating an integration step. - + // Generate a forcing vector before canculating an integration step. + // Generate forcing wave vector k_force AcReal3 k_force; k_force = helical_forcing_k_generator(kmax, kmin); - //Randomize the phase - AcReal phase = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); + // Randomize the phase + AcReal phase = AcReal(2.0) * AcReal(M_PI) * get_random_number_01(); - // Generate e for k. Needed for the sake of isotrophy. + // Generate e for k. Needed for the sake of isotrophy. AcReal3 e_force; if ((k_force.y == 0.0) && (k_force.z == 0.0)) { e_force = (AcReal3){0.0, 1.0, 0.0}; - } else { + } + else { e_force = (AcReal3){1.0, 0.0, 0.0}; } helical_forcing_e_generator(&e_force, k_force); From af3a1e211e3c5bfdfd0f0f8d0447d99152fc6b4e Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 17:03:26 +0300 Subject: [PATCH 25/32] Suppressed unused variable and function warnings in model_rk3.cc --- src/standalone/model/model_rk3.cc | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc index 6c79908..227a915 100644 --- a/src/standalone/model/model_rk3.cc +++ b/src/standalone/model/model_rk3.cc @@ -697,32 +697,31 @@ entropy(const ModelScalarData& ss, const ModelVectorData& uu, const ModelScalarD */ } -static bool +static inline bool is_valid(const ModelScalar a) { return !isnan(a) && !isinf(a); } -static bool +static inline bool is_valid(const ModelVector& a) { return is_valid(a.x) && is_valid(a.y) && is_valid(a.z); } - #if 0 //FORCING NOT SUPPORTED FOR AUTOTEST static inline ModelVector simple_vortex_forcing(ModelVector a, ModelVector b, ModelScalar magnitude) { - return magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex + return magnitude * cross(normalized(b - a), (ModelVector){0, 0, 1}); // Vortex } static inline ModelVector simple_outward_flow_forcing(ModelVector a, ModelVector b, ModelScalar magnitude) { - return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow + return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow } @@ -746,7 +745,7 @@ helical_forcing(ModelScalar magnitude, ModelVector k_force, ModelVector xx, Mode 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, - ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase}; + ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase}; return force; } @@ -765,20 +764,20 @@ forcing(int3 globalVertexIdx, ModelScalar dt) //Placeholders until determined properly ModelScalar magnitude = get(AC_forcing_magnitude); - ModelScalar phase = get(AC_forcing_phase); + ModelScalar phase = get(AC_forcing_phase); ModelVector k_force = (ModelVector){ get(AC_k_forcex), get(AC_k_forcey), get(AC_k_forcez)}; 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)}; - //Determine that forcing funtion type at this point. + //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); ModelVector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase); //Scaling N = magnitude*cs*sqrt(k*cs/dt) * dt - const ModelScalar NN = cs*sqrt(get(AC_kaver)*cs); - //MV: Like in the Pencil Code. I don't understandf the logic here. + const ModelScalar NN = cs*sqrt(get(AC_kaver)*cs); + //MV: Like in the Pencil Code. I don't understandf the logic here. force.x = sqrt(dt)*NN*force.x; force.y = sqrt(dt)*NN*force.y; force.z = sqrt(dt)*NN*force.z; @@ -849,12 +848,13 @@ solve_beta_step(const int step_number, const ModelScalar dt, const int i, const for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) out->vertex_buffer[w][idx] += beta[step_number] * in.vertex_buffer[w][idx]; + (void)dt; // Suppress unused variable warning if forcing not used #if LFORCING if (step_number == 2) { ModelVector force = forcing((int3){i, j, k}, dt); - out->vertex_buffer[VTXBUF_UUX][idx] += force.x ; - out->vertex_buffer[VTXBUF_UUY][idx] += force.y ; - out->vertex_buffer[VTXBUF_UUZ][idx] += force.z ; + out->vertex_buffer[VTXBUF_UUX][idx] += force.x; + out->vertex_buffer[VTXBUF_UUY][idx] += force.y; + out->vertex_buffer[VTXBUF_UUZ][idx] += force.z; } #endif } From 945751e58554ef733683d19c0431642cb89f4de4 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 17:06:57 +0300 Subject: [PATCH 26/32] Autoformatted host_forcing.cc --- src/standalone/model/host_forcing.cc | 113 +++++++++++++-------------- 1 file changed, 55 insertions(+), 58 deletions(-) diff --git a/src/standalone/model/host_forcing.cc b/src/standalone/model/host_forcing.cc index ecac783..3c4c714 100644 --- a/src/standalone/model/host_forcing.cc +++ b/src/standalone/model/host_forcing.cc @@ -28,17 +28,14 @@ #include "core/math_utils.h" - -//The is a wrapper for genering random numbers with a chosen system. +// The is a wrapper for genering random numbers with a chosen system. AcReal get_random_number_01() { - //TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/ - return AcReal(rand())/AcReal(RAND_MAX); + // TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/ + return AcReal(rand()) / AcReal(RAND_MAX); } - - AcReal3 cross(const AcReal3& a, const AcReal3& b) { @@ -63,9 +60,9 @@ vec_norm(const AcReal3& a) AcReal3 c; AcReal norm = dot(a, a); - c.x = a.x/sqrt(norm); - c.y = a.y/sqrt(norm); - c.z = a.z/sqrt(norm); + c.x = a.x / sqrt(norm); + c.y = a.y / sqrt(norm); + c.z = a.z / sqrt(norm); return c; } @@ -75,9 +72,9 @@ vec_multi_scal(const AcReal scal, const AcReal3& a) { AcReal3 c; - c.x = a.x*scal; - c.y = a.y*scal; - c.z = a.z*scal; + c.x = a.x * scal; + c.y = a.y * scal; + c.z = a.z * scal; return c; } @@ -86,95 +83,95 @@ vec_multi_scal(const AcReal scal, const AcReal3& a) AcReal3 helical_forcing_k_generator(const AcReal kmax, const AcReal kmin) { - AcReal phi, theta, kk; //Spherical direction coordinates - AcReal3 k_force; //forcing wave vector + AcReal phi, theta, kk; // Spherical direction coordinates + AcReal3 k_force; // forcing wave vector AcReal delta_k = kmax - kmin; // Generate vector in spherical coordinates - phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); - theta = AcReal(M_PI)*get_random_number_01(); - kk = delta_k*get_random_number_01() + kmin; + phi = AcReal(2.0) * AcReal(M_PI) * get_random_number_01(); + theta = AcReal(M_PI) * get_random_number_01(); + kk = delta_k * get_random_number_01() + kmin; // Cast into Cartesian form - k_force = (AcReal3){kk*sin(theta)*cos(phi), - kk*sin(theta)*sin(phi), - kk*cos(theta) }; + k_force = (AcReal3){kk * sin(theta) * cos(phi), // + 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); + // 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); + // 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); + // 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; } -//Generate the unit perpendicular unit vector e required for helical forcing -//Addapted from Pencil code forcing.f90 hel_vec() subroutine. +// Generate the unit perpendicular unit vector e required for helical forcing +// Addapted from Pencil code forcing.f90 hel_vec() subroutine. void helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force) { - AcReal3 k_cross_e = cross(k_force, *e_force); - k_cross_e = vec_norm(k_cross_e); + AcReal3 k_cross_e = cross(k_force, *e_force); + k_cross_e = vec_norm(k_cross_e); AcReal3 k_cross_k_cross_e = cross(k_force, k_cross_e); - k_cross_k_cross_e = vec_norm(k_cross_k_cross_e); - AcReal phi = AcReal(2.0)*AcReal(M_PI)*get_random_number_01(); - AcReal3 ee_tmp1 = vec_multi_scal(cos(phi),k_cross_e); - AcReal3 ee_tmp2 = vec_multi_scal(sin(phi), k_cross_k_cross_e); + k_cross_k_cross_e = vec_norm(k_cross_k_cross_e); + AcReal phi = AcReal(2.0) * AcReal(M_PI) * get_random_number_01(); + AcReal3 ee_tmp1 = vec_multi_scal(cos(phi), k_cross_e); + AcReal3 ee_tmp2 = vec_multi_scal(sin(phi), k_cross_k_cross_e); - *e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, - ee_tmp1.y + ee_tmp2.y, - ee_tmp1.z + ee_tmp2.z}; + *e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, ee_tmp1.y + ee_tmp2.y, ee_tmp1.z + ee_tmp2.z}; } -//PC Manual Eq. 223 +// PC Manual Eq. 223 void -helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force, +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 + // k dot e AcReal3 kdote; kdote.x = k_force.x * e_force.x; kdote.y = k_force.y * e_force.y; kdote.z = k_force.z * e_force.z; // k cross e - AcReal3 k_cross_e; - k_cross_e.x=k_force.y*e_force.z-k_force.z*e_force.y; - k_cross_e.y=k_force.z*e_force.x-k_force.x*e_force.z; - k_cross_e.z=k_force.x*e_force.y-k_force.y*e_force.x; + AcReal3 k_cross_e; + k_cross_e.x = k_force.y * e_force.z - k_force.z * e_force.y; + k_cross_e.y = k_force.z * e_force.x - k_force.x * e_force.z; + k_cross_e.z = k_force.x * e_force.y - k_force.y * e_force.x; // k cross k cross e - AcReal3 k_cross_k_cross_e; - k_cross_k_cross_e.x=k_force.y*k_cross_e.z-k_force.z*k_cross_e.y; - k_cross_k_cross_e.y=k_force.z*k_cross_e.x-k_force.x*k_cross_e.z; - k_cross_k_cross_e.z=k_force.x*k_cross_e.y-k_force.y*k_cross_e.x; + AcReal3 k_cross_k_cross_e; + k_cross_k_cross_e.x = k_force.y * k_cross_e.z - k_force.z * k_cross_e.y; + k_cross_k_cross_e.y = k_force.z * k_cross_e.x - k_force.x * k_cross_e.z; + k_cross_k_cross_e.z = k_force.x * k_cross_e.y - k_force.y * k_cross_e.x; // abs(k) - AcReal kabs = sqrt(k_force.x*k_force.x + k_force.y*k_force.y + k_force.z*k_force.z); + AcReal kabs = sqrt(k_force.x * k_force.x + k_force.y * k_force.y + k_force.z * k_force.z); - AcReal denominator = sqrt(AcReal(1.0) + relhel*relhel)*kabs - *sqrt(kabs*kabs - (kdote.x*kdote.x + kdote.y*kdote.y + kdote.z*kdote.z)); + AcReal denominator = sqrt(AcReal(1.0) + relhel * relhel) * kabs * + sqrt(kabs * kabs - + (kdote.x * kdote.x + kdote.y * kdote.y + kdote.z * kdote.z)); - //MV: I suspect there is a typo in the Pencil Code manual! - //*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator, + // MV: I suspect there is a typo in the Pencil Code manual! + //*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator, // -relhel*kabs*k_cross_e.y/denominator, // -relhel*kabs*k_cross_e.z/denominator}; - //*ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator, + //*ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator, // k_cross_k_cross_e.y/denominator, // k_cross_k_cross_e.z/denominator}; // 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}; + *ff_hel_re = (AcReal3){kabs * k_cross_e.x / denominator, kabs * k_cross_e.y, + kabs * k_cross_e.z}; - *ff_hel_im = (AcReal3){relhel*k_cross_k_cross_e.x/denominator, - relhel*k_cross_k_cross_e.y, - relhel*k_cross_k_cross_e.z}; + *ff_hel_im = (AcReal3){relhel * k_cross_k_cross_e.x / denominator, relhel * k_cross_k_cross_e.y, + relhel * k_cross_k_cross_e.z}; } From d54ccc1da8b8227cffad97b950243c85c0f91b18 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 17:10:01 +0300 Subject: [PATCH 27/32] Deprecated a block of old code that was used a long time ago for testing forcing --- src/core/kernels/kernels.cuh | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/core/kernels/kernels.cuh b/src/core/kernels/kernels.cuh index 1ff9051..4c990fb 100644 --- a/src/core/kernels/kernels.cuh +++ b/src/core/kernels/kernels.cuh @@ -739,17 +739,20 @@ read_out(const int idx, AcReal* __restrict__ field[], const int3 handle) * ============================================================================= */ +/* // DEPRECATED static AcReal randf(void) { return AcReal(rand()) / AcReal(RAND_MAX); } +*/ AcResult rk3_step_async(const cudaStream_t stream, const int& step_number, const int3& start, const int3& end, const AcReal dt, VertexBufferArray* buffer) { const dim3 tpb(32, 1, 4); + /* // DEPRECATED /////////////////// Forcing #if LFORCING const AcReal ff_scale = AcReal(.2); @@ -764,6 +767,7 @@ rk3_step_async(const cudaStream_t stream, const int& step_number, const int3& st stream); #endif // LFORCING ////////////////////////// + */ const int nx = end.x - start.x; const int ny = end.y - start.y; From 919d446222b8526fe1d5f4ed7fb50dc9d463c304 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 17:10:40 +0300 Subject: [PATCH 28/32] Commented out unused variables in simulation.cc --- src/standalone/simulation.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index 76cc40f..156e23a 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -238,12 +238,12 @@ run_simulation(void) AcReal bin_crit_t = bin_save_t; // Forcing properties - AcReal relhel = mesh_info.real_params[AC_relhel]; - AcReal magnitude = mesh_info.real_params[AC_forcing_magnitude]; - AcReal kmin = mesh_info.real_params[AC_kmin]; - AcReal kmax = mesh_info.real_params[AC_kmax]; + // AcReal relhel = mesh_info.real_params[AC_relhel]; // JP: gcc warning: unused + // AcReal magnitude = mesh_info.real_params[AC_forcing_magnitude]; // JP: gcc warning: unused + // AcReal kmin = mesh_info.real_params[AC_kmin]; // JP: gcc warning: unused + // AcReal kmax = mesh_info.real_params[AC_kmax]; // JP: gcc warning: unused - AcReal kaver = (kmax - kmin) / AcReal(2.0); + // AcReal kaver = (kmax - kmin) / AcReal(2.0); // JP: gcc warning: unused /* initialize random seed: */ srand(312256655); From d7228f06473effa3c3b7a2ed1fcf8476845113f6 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 17:11:26 +0300 Subject: [PATCH 29/32] Added an explicit cast from double to AcReal to avoid a narrowing conversion error --- src/standalone/model/host_forcing.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/standalone/model/host_forcing.cc b/src/standalone/model/host_forcing.cc index 3c4c714..4ce7892 100644 --- a/src/standalone/model/host_forcing.cc +++ b/src/standalone/model/host_forcing.cc @@ -94,9 +94,9 @@ helical_forcing_k_generator(const AcReal kmax, const AcReal kmin) kk = delta_k * get_random_number_01() + kmin; // Cast into Cartesian form - k_force = (AcReal3){kk * sin(theta) * cos(phi), // - kk * sin(theta) * sin(phi), // - kk * cos(theta)}; + k_force = (AcReal3){AcReal(kk * sin(theta) * cos(phi)), // + AcReal(kk * sin(theta) * sin(phi)), // + AcReal(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); From 8ed947ce98ed04479cbc3dad96f9644a4fe4b124 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 17:13:45 +0300 Subject: [PATCH 30/32] Removed deprecated sinusoidal forcing from kernels.cuh --- src/core/kernels/kernels.cuh | 52 ------------------------------------ 1 file changed, 52 deletions(-) diff --git a/src/core/kernels/kernels.cuh b/src/core/kernels/kernels.cuh index 4c990fb..faf6917 100644 --- a/src/core/kernels/kernels.cuh +++ b/src/core/kernels/kernels.cuh @@ -543,34 +543,6 @@ normalized(const AcReal3& vec) return inv_len * vec; } -// Sinusoidal forcing -// https://arxiv.org/pdf/1704.04676.pdf -// NOTE: This method of forcing is depracated. However, it will remain in here -// until a corresponding scheme exists in the new code. -__constant__ AcReal3 forcing_vec; -__constant__ AcReal forcing_phi; -static __device__ __forceinline__ AcReal3 -DEPRECATED_forcing(const int i, const int j, const int k) -{ -#define DOMAIN_SIZE_X (DCONST_INT(AC_nx) * DCONST_REAL(AC_dsx)) -#define DOMAIN_SIZE_Y (DCONST_INT(AC_ny) * DCONST_REAL(AC_dsy)) -#define DOMAIN_SIZE_Z (DCONST_INT(AC_nz) * DCONST_REAL(AC_dsz)) - const AcReal3 k_vec = (AcReal3){ - (i - DCONST_INT(AC_nx_min)) * DCONST_REAL(AC_dsx) - AcReal(.5) * DOMAIN_SIZE_X, - (j - DCONST_INT(AC_ny_min)) * DCONST_REAL(AC_dsy) - AcReal(.5) * DOMAIN_SIZE_Y, - (k - DCONST_INT(AC_nz_min)) * DCONST_REAL(AC_dsz) - AcReal(.5) * DOMAIN_SIZE_Z}; - AcReal inv_len = reciprocal_len(k_vec); - if (isnan(inv_len) || isinf(inv_len)) - inv_len = 0; - if (inv_len > 2) // hack to make it cool - inv_len = 2; - const AcReal k_dot_x = dot(k_vec, forcing_vec); - - const AcReal waves = cos(k_dot_x) * cos(forcing_phi) - sin(k_dot_x) * sin(forcing_phi); - - return inv_len * inv_len * waves * forcing_vec; -} - // Note: LNT0 and LNRHO0 must be set very carefully: if the magnitude is different that other values // in the mesh, then we will inherently lose precision #define LNT0 (AcReal(0.0)) @@ -739,35 +711,11 @@ read_out(const int idx, AcReal* __restrict__ field[], const int3 handle) * ============================================================================= */ -/* // DEPRECATED -static AcReal -randf(void) -{ - return AcReal(rand()) / AcReal(RAND_MAX); -} -*/ - AcResult rk3_step_async(const cudaStream_t stream, const int& step_number, const int3& start, const int3& end, const AcReal dt, VertexBufferArray* buffer) { const dim3 tpb(32, 1, 4); - /* // DEPRECATED -/////////////////// Forcing -#if LFORCING - const AcReal ff_scale = AcReal(.2); - static AcReal3 ff = ff_scale * (AcReal3){1, 0, 0}; - const AcReal radians = randf() * AcReal(2 * M_PI) / 360 / 8; - const AcMatrix rotz = create_rotz(radians); - ff = mul(rotz, ff); - cudaMemcpyToSymbolAsync(forcing_vec, &ff, sizeof(ff), 0, cudaMemcpyHostToDevice, stream); - - const AcReal ff_phi = AcReal(M_PI); // AcReal(2 * M_PI) * randf(); - cudaMemcpyToSymbolAsync(forcing_phi, &ff_phi, sizeof(ff_phi), 0, cudaMemcpyHostToDevice, - stream); -#endif // LFORCING - ////////////////////////// - */ const int nx = end.x - start.x; const int ny = end.y - start.y; From 81a09501b86e51a647a678158869ac0b88c60e30 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 17:23:37 +0300 Subject: [PATCH 31/32] Removed deprecated LNT0 and LNRHO0 defines, now the actual configuration parameters are used (AC_lnrho0 and AC_lnT0). Also accidental autoformatting again, there seems to be stray spaces before linebreaks in some files which get automatically removed by my text editor --- acc/mhd_solver/stencil_process.sps | 27 +++++++++++++++------------ src/core/kernels/kernels.cuh | 5 ----- src/standalone/model/model_rk3.cc | 19 +++++++------------ 3 files changed, 22 insertions(+), 29 deletions(-) diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 8fbd8ae..e0b3080 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -20,6 +20,9 @@ uniform Scalar dsx; uniform Scalar dsy; uniform Scalar dsz; +uniform Scalar lnT0; +uniform Scalar lnrho0; + uniform int nx_min; uniform int ny_min; uniform int nz_min; @@ -54,9 +57,9 @@ gradients(in Vector uu) Scalar continuity(in Vector uu, in Scalar lnrho) { - return -dot(value(uu), gradient(lnrho)) + return -dot(value(uu), gradient(lnrho)) #if LUPWD - //This is a corrective hyperdiffusion term for upwinding. + //This is a corrective hyperdiffusion term for upwinding. + upwd_der6(uu, lnrho) #endif - divergence(uu); @@ -66,7 +69,7 @@ continuity(in Vector uu, in Scalar lnrho) { Vector momentum(in Vector uu, in Scalar lnrho, in Scalar ss, in Vector aa) { const Matrix S = stress_tensor(uu); - const Scalar cs2 = cs2_sound * exp(gamma * value(ss) / cp_sound + (gamma - 1) * (value(lnrho) - LNRHO0)); + const Scalar cs2 = cs2_sound * exp(gamma * value(ss) / cp_sound + (gamma - 1) * (value(lnrho) - lnrho0)); const Vector j = (Scalar(1.) / mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density const Vector B = curl(aa); //TODO: DOES INTHERMAL VERSTION INCLUDE THE MAGNETIC FIELD? @@ -150,8 +153,8 @@ induction(in Vector uu, in Vector aa) { #if LENTROPY Scalar lnT( in Scalar ss, in Scalar lnrho) { - const Scalar lnT = LNT0 + gamma * value(ss) / cp_sound + - (gamma - Scalar(1.)) * (value(lnrho) - LNRHO0); + const Scalar lnT = lnT0 + gamma * value(ss) / cp_sound + + (gamma - Scalar(1.)) * (value(lnrho) - lnrho0); return lnT; } @@ -213,13 +216,13 @@ heat_transfer(in Vector uu, in Scalar lnrho, in Scalar tt) Vector simple_vortex_forcing(Vector a, Vector b, Scalar magnitude) { - return magnitude * cross(normalized(b - a), (Vector){0, 0, 1}); // Vortex + 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 + return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow } @@ -245,7 +248,7 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto Vector force = (Vector){ ff_re.x*real_comp_phase - ff_im.x*imag_comp_phase, ff_re.y*real_comp_phase - ff_im.y*imag_comp_phase, - ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase}; + ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase}; return force; } @@ -264,20 +267,20 @@ forcing(int3 globalVertexIdx, Scalar dt) //Placeholders until determined properly Scalar magnitude = DCONST_REAL(AC_forcing_magnitude); - Scalar phase = DCONST_REAL(AC_forcing_phase); + Scalar phase = DCONST_REAL(AC_forcing_phase); Vector k_force = (Vector){ DCONST_REAL(AC_k_forcex), DCONST_REAL(AC_k_forcey), DCONST_REAL(AC_k_forcez)}; Vector ff_re = (Vector){DCONST_REAL(AC_ff_hel_rex), DCONST_REAL(AC_ff_hel_rey), DCONST_REAL(AC_ff_hel_rez)}; Vector ff_im = (Vector){DCONST_REAL(AC_ff_hel_imx), DCONST_REAL(AC_ff_hel_imy), DCONST_REAL(AC_ff_hel_imz)}; - //Determine that forcing funtion type at this point. + //Determine that forcing funtion type at this point. //Vector force = simple_vortex_forcing(a, xx, magnitude); //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) * dt - const Scalar NN = cs*sqrt(DCONST_REAL(AC_kaver)*cs); - //MV: Like in the Pencil Code. I don't understandf the logic here. + const Scalar NN = cs*sqrt(DCONST_REAL(AC_kaver)*cs); + //MV: Like in the Pencil Code. I don't understandf the logic here. force.x = sqrt(dt)*NN*force.x; force.y = sqrt(dt)*NN*force.y; force.z = sqrt(dt)*NN*force.z; diff --git a/src/core/kernels/kernels.cuh b/src/core/kernels/kernels.cuh index faf6917..a8d4d50 100644 --- a/src/core/kernels/kernels.cuh +++ b/src/core/kernels/kernels.cuh @@ -543,11 +543,6 @@ normalized(const AcReal3& vec) return inv_len * vec; } -// Note: LNT0 and LNRHO0 must be set very carefully: if the magnitude is different that other values -// in the mesh, then we will inherently lose precision -#define LNT0 (AcReal(0.0)) -#define LNRHO0 (AcReal(0.0)) - #define H_CONST (AcReal(0.0)) #define C_CONST (AcReal(0.0)) diff --git a/src/standalone/model/model_rk3.cc b/src/standalone/model/model_rk3.cc index 227a915..71b3f2d 100644 --- a/src/standalone/model/model_rk3.cc +++ b/src/standalone/model/model_rk3.cc @@ -553,11 +553,6 @@ normalized(const ModelVector& vec) return inv_len * vec; } -// Note: LNT0 and LNRHO0 must be set very carefully: if the magnitude is different that other values -// in the mesh, then we will inherently lose precision -#define LNT0 (ModelScalar(0.0)) -#define LNRHO0 (ModelScalar(0.0)) - #define H_CONST (ModelScalar(0.0)) #define C_CONST (ModelScalar(0.0)) @@ -571,9 +566,10 @@ momentum(const ModelVectorData& uu, const ModelScalarData& lnrho { #if LENTROPY const ModelMatrix S = stress_tensor(uu); - const ModelScalar cs2 = get(AC_cs2_sound) * expl(get(AC_gamma) * value(ss) / get(AC_cp_sound) + - (get(AC_gamma) - 1) * (value(lnrho) - LNRHO0)); - const ModelVector j = (ModelScalar(1.) / get(AC_mu0)) * + const ModelScalar cs2 = get(AC_cs2_sound) * + expl(get(AC_gamma) * value(ss) / get(AC_cp_sound) + + (get(AC_gamma) - 1) * (value(lnrho) - get(AC_lnrho0))); + const ModelVector j = (ModelScalar(1.) / get(AC_mu0)) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density const ModelVector B = curl(aa); const ModelScalar inv_rho = ModelScalar(1.) / expl(value(lnrho)); @@ -593,9 +589,8 @@ momentum(const ModelVectorData& uu, const ModelScalarData& lnrho const ModelMatrix S = stress_tensor(uu); //#if LENTROPY - //const ModelScalar lnrho0 = 1; // TODO correct lnrho0 const ModelScalar cs02 = get(AC_cs2_sound); // TODO better naming - const ModelScalar cs2 = cs02;// * expl(get(AC_gamma) * value(ss) / get(AC_cp_sound) + (get(AC_gamma)-ModelScalar(1.l)) * (value(lnrho) - lnrho0)); + const ModelScalar cs2 = cs02;// * expl(get(AC_gamma) * value(ss) / get(AC_cp_sound) + (get(AC_gamma)-ModelScalar(1.l)) * (value(lnrho) - get(AC_lnrho0))); mom = -mul(gradients(uu), value(uu)) - cs2 * ((ModelScalar(1.) / get(AC_cp_sound)) * gradient(ss) + gradient(lnrho)) + @@ -642,8 +637,8 @@ induction(const ModelVectorData& uu, const ModelVectorData& aa) static inline ModelScalar lnT(const ModelScalarData& ss, const ModelScalarData& lnrho) { - const ModelScalar lnT = LNT0 + get(AC_gamma) * value(ss) / get(AC_cp_sound) + - (get(AC_gamma) - ModelScalar(1.)) * (value(lnrho) - LNRHO0); + const ModelScalar lnT = get(AC_lnT0) + get(AC_gamma) * value(ss) / get(AC_cp_sound) + + (get(AC_gamma) - ModelScalar(1.)) * (value(lnrho) - get(AC_lnrho0)); return lnT; } From e8a5579b5018b15c6b38394e90e1b9388eebb451 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Wed, 3 Jul 2019 17:25:26 +0300 Subject: [PATCH 32/32] Made the gcc error flags more lenient temporarily since there are so many float-double-float conversion errors in host_forcing.cc --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 533f38a..37225d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -129,7 +129,7 @@ set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG}\ -O0 -g") if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") - set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror -Wdouble-promotion -Wfloat-conversion") + set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror")# -Wdouble-promotion -Wfloat-conversion") elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") set (CXX_FLAGS_WARNING "-Wall -Wextra -Werror-all -Wsign-conversion") else()