cwpearson/astaroth
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Tryin to prepare autotest for forcing.

Browse Source
This commit is contained in:
Miikka Vaisala
2019-07-02 18:24:41 +08:00
parent 03689709df
commit 4766441ffb
4 changed files with 81 additions and 174 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
acc/mhd_solver/stencil_process.sps
View File

@@ -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)

2
src/standalone/autotest.cc
View File

@@ -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);

101
src/standalone/model/model_rk3.cc
View File

@@ -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
}

150
src/standalone/simulation.cc
View File

@@ -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 <string.h>
@@ -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