diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index f476921..39527e0 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -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 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 = 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_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}; + 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; } @@ -278,8 +278,8 @@ forcing(int3 globalVertexIdx, Scalar dt) const Scalar NN = cs*sqrt(DCONST_REAL(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 = force.y; - force.z = force.z; + force.y = sqrt(dt)*NN*force.y; + force.z = sqrt(dt)*NN*force.z; if (is_valid(force)) { return force; } else { return (Vector){0, 0, 0}; } diff --git a/src/standalone/model/host_timestep.cc b/src/standalone/model/host_timestep.cc index 9c51d83..54f47cc 100644 --- a/src/standalone/model/host_timestep.cc +++ b/src/standalone/model/host_timestep.cc @@ -50,7 +50,8 @@ host_timestep(const AcReal& umax, const AcMeshInfo& mesh_info) // New, closer to the actual Courant timestep // 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 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); return AcReal(timescale) * AcReal(dt); diff --git a/src/standalone/simulation.cc b/src/standalone/simulation.cc index dc90dc0..659e382 100644 --- a/src/standalone/simulation.cc +++ b/src/standalone/simulation.cc @@ -386,8 +386,8 @@ run_simulation(void) AcReal bin_crit_t = bin_save_t; //Placeholders until determined properly - AcReal magnitude = 0.05; - AcReal relhel = 0.5; + AcReal magnitude = 1e-5; + AcReal relhel = 0.0; AcReal kmin = 0.8; AcReal kmax = 1.2;