From 4766441ffbbaa1aeb5d7cfb1cc5e7b1187c0f4c2 Mon Sep 17 00:00:00 2001 From: Miikka Vaisala Date: Tue, 2 Jul 2019 18:24:41 +0800 Subject: [PATCH] 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