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

Corrected a bug in the timestep and some scaling problems.

Browse Source 
Now I can reach a saturated stated in forcing without crashing the code.
This commit is contained in:
Miikka Vaisala
2019-07-01 14:19:56 +08:00
parent 9f0be0d9ff
commit 0600790f41
3 changed files with 10 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
acc/mhd_solver/stencil_process.sps
View File

@@ -236,15 +236,15 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto
Scalar cos_k_dox_x = cos(dot(k_force, xx)); Scalar cos_k_dox_x = cos(dot(k_force, xx));
Scalar sin_k_dox_x = sin(dot(k_force, xx)); Scalar sin_k_dox_x = sin(dot(k_force, xx));
// Phase affect only the x-component // Phase affect only the x-component
Scalar real_comp = cos_k_dox_x; //Scalar real_comp = cos_k_dox_x;
Scalar imag_comp = sin_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 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; Scalar imag_comp_phase = cos_k_dox_x*sin_phi + sin_k_dox_x*cos_phi;
Vector force = (Vector){ ff_re.x*real_comp_phase - ff_im.x*imag_comp_phase, 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.y*real_comp_phase - ff_im.y*imag_comp_phase,
ff_re.z*real_comp - ff_im.z*imag_comp}; ff_re.z*real_comp_phase - ff_im.z*imag_comp_phase};
return force; return force;
} }
@@ -278,8 +278,8 @@ forcing(int3 globalVertexIdx, Scalar dt)
const Scalar NN = cs*sqrt(DCONST_REAL(AC_kaver)*cs); const Scalar NN = cs*sqrt(DCONST_REAL(AC_kaver)*cs);
//MV: Like in the Pencil Code. I don't understandf the logic here. //MV: Like in the Pencil Code. I don't understandf the logic here.
force.x = sqrt(dt)*NN*force.x; force.x = sqrt(dt)*NN*force.x;
force.y = force.y; force.y = sqrt(dt)*NN*force.y;
force.z = force.z; force.z = sqrt(dt)*NN*force.z;
if (is_valid(force)) { return force; } if (is_valid(force)) { return force; }
else { return (Vector){0, 0, 0}; } else { return (Vector){0, 0, 0}; }

3
src/standalone/model/host_timestep.cc
View File

@@ -50,7 +50,8 @@ host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info)
// New, closer to the actual Courant timestep // New, closer to the actual Courant timestep
// See Pencil Code user manual p. 38 (timestep section) // See Pencil Code user manual p. 38 (timestep section)
const long double uu_dt = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + 0.0l)); const long double uu_dt = cdt * dsmin / (fabsl(umax) + sqrtl(cs2_sound + 0.0l));
const long double visc_dt = cdtv * dsmin * dsmin / max(max(nu_visc, eta), max(gamma, chi)) + 1; // TODO NOTE: comment the +1 out to get scientifically accurate results const long double visc_dt = cdtv * dsmin * dsmin / max(max(nu_visc, eta), max(gamma, chi));// + 1; // TODO NOTE: comment the +1 out to get scientifically accurate results
//MV: White the +1? It was messing up my computations!
const long double dt = min(uu_dt, visc_dt); const long double dt = min(uu_dt, visc_dt);
return AcReal(timescale) * AcReal(dt); return AcReal(timescale) * AcReal(dt);

4
src/standalone/simulation.cc
View File

@@ -386,8 +386,8 @@ run_simulation(void)
AcReal bin_crit_t = bin_save_t; AcReal bin_crit_t = bin_save_t;
//Placeholders until determined properly //Placeholders until determined properly
AcReal magnitude = 0.05; AcReal magnitude = 1e-5;
AcReal relhel = 0.5; AcReal relhel = 0.0;
AcReal kmin = 0.8; AcReal kmin = 0.8;
AcReal kmax = 1.2; AcReal kmax = 1.2;