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

Forcing function issue not yet fully resolved.

Browse Source 
Now brain hurs. No more today. Break needed.
This commit is contained in:
Miikka Vaisala
2019-06-28 19:23:18 +08:00
parent 94a25383a9
commit f04ef8e64c
2 changed files with 39 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
acc/mhd_solver/stencil_process.sps
View File

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

29
src/standalone/simulation.cc
View File

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