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

Changed #if 0 to #if LFORCING instead to get the code to compile if forcing is used (even though autotesting does not support it yet). Also more autoformatting. Maybe I should disable it or then everyone should start using it to avoid cluttering commits with these superficial changes

Browse Source
This commit is contained in:
jpekkila
2019-07-03 17:49:34 +03:00
parent 609cfaea14
commit b4eea4b6b6
1 changed files with 48 additions and 42 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -687,7 +687,7 @@ is_valid(const ModelVector& a)
return is_valid(a.x) && is_valid(a.y) && is_valid(a.z);
}
#if 0
#if LFORCING
// FORCING NOT SUPPORTED FOR AUTOTEST
static inline ModelVector
@@ -702,10 +702,11 @@ 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
// 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)
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))));
@@ -720,7 +721,6 @@ helical_forcing(ModelScalar magnitude, ModelVector k_force, ModelVector xx, Mode
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};
@@ -731,10 +731,13 @@ helical_forcing(ModelScalar magnitude, ModelVector k_force, ModelVector xx, Mode
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),
*/
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);
@@ -747,7 +750,6 @@ forcing(int3 globalVertexIdx, ModelScalar dt)
ModelVector ff_re = (ModelVector){get(AC_ff_hel_rex), get(AC_ff_hel_rey), get(AC_ff_hel_rez)};
ModelVector ff_im = (ModelVector){get(AC_ff_hel_imx), get(AC_ff_hel_imy), get(AC_ff_hel_imz)};
// 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);
@@ -760,8 +762,12 @@ forcing(int3 globalVertexIdx, ModelScalar dt)
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}; }
if (is_valid(force)) {
return force;
}
else {
return (ModelVector){0, 0, 0};
}
}
#endif