From 9f0be0d9ff333556c1c9c518bc3eab132d396472 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Mon, 1 Jul 2019 11:06:42 +0800 Subject: [PATCH] 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);