diff --git a/acc/mhd_solver/stencil_assembly.sas b/acc/mhd_solver/stencil_assembly.sas index e30d500..85886e6 100644 --- a/acc/mhd_solver/stencil_assembly.sas +++ b/acc/mhd_solver/stencil_assembly.sas @@ -9,9 +9,7 @@ value(in ScalarField vertex) Preprocessed Vector gradient(in ScalarField vertex) { - return (Vector){derx(vertexIdx, vertex), - dery(vertexIdx, vertex), - derz(vertexIdx, vertex)}; + return (Vector){derx(vertexIdx, vertex), dery(vertexIdx, vertex), derz(vertexIdx, vertex)}; } #if LUPWD @@ -21,14 +19,14 @@ der6x_upwd(in ScalarField vertex) { Scalar inv_ds = AC_inv_dsx; - return (Scalar){ Scalar(1.0/60.0)*inv_ds* ( - - Scalar(20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] - + Scalar(15.0)*(vertex[vertexIdx.x+1, vertexIdx.y, vertexIdx.z] - + vertex[vertexIdx.x-1, vertexIdx.y, vertexIdx.z]) - - Scalar( 6.0)*(vertex[vertexIdx.x+2, vertexIdx.y, vertexIdx.z] - + vertex[vertexIdx.x-2, vertexIdx.y, vertexIdx.z]) - + vertex[vertexIdx.x+3, vertexIdx.y, vertexIdx.z] - + vertex[vertexIdx.x-3, vertexIdx.y, vertexIdx.z])}; + return (Scalar){Scalar(1.0 / 60.0) * inv_ds * + (-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] + + Scalar(15.0) * (vertex[vertexIdx.x + 1, vertexIdx.y, vertexIdx.z] + + vertex[vertexIdx.x - 1, vertexIdx.y, vertexIdx.z]) - + Scalar(6.0) * (vertex[vertexIdx.x + 2, vertexIdx.y, vertexIdx.z] + + vertex[vertexIdx.x - 2, vertexIdx.y, vertexIdx.z]) + + vertex[vertexIdx.x + 3, vertexIdx.y, vertexIdx.z] + + vertex[vertexIdx.x - 3, vertexIdx.y, vertexIdx.z])}; } Preprocessed Scalar @@ -36,14 +34,14 @@ der6y_upwd(in ScalarField vertex) { Scalar inv_ds = AC_inv_dsy; - return (Scalar){ Scalar(1.0/60.0)*inv_ds* ( - -Scalar( 20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] - +Scalar( 15.0)*(vertex[vertexIdx.x, vertexIdx.y+1, vertexIdx.z] - + vertex[vertexIdx.x, vertexIdx.y-1, vertexIdx.z]) - -Scalar( 6.0)*(vertex[vertexIdx.x, vertexIdx.y+2, vertexIdx.z] - + vertex[vertexIdx.x, vertexIdx.y-2, vertexIdx.z]) - + vertex[vertexIdx.x, vertexIdx.y+3, vertexIdx.z] - + vertex[vertexIdx.x, vertexIdx.y-3, vertexIdx.z])}; + return (Scalar){Scalar(1.0 / 60.0) * inv_ds * + (-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] + + Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y + 1, vertexIdx.z] + + vertex[vertexIdx.x, vertexIdx.y - 1, vertexIdx.z]) - + Scalar(6.0) * (vertex[vertexIdx.x, vertexIdx.y + 2, vertexIdx.z] + + vertex[vertexIdx.x, vertexIdx.y - 2, vertexIdx.z]) + + vertex[vertexIdx.x, vertexIdx.y + 3, vertexIdx.z] + + vertex[vertexIdx.x, vertexIdx.y - 3, vertexIdx.z])}; } Preprocessed Scalar @@ -51,14 +49,14 @@ der6z_upwd(in ScalarField vertex) { Scalar inv_ds = AC_inv_dsz; - return (Scalar){ Scalar(1.0/60.0)*inv_ds* ( - -Scalar( 20.0)* vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] - +Scalar( 15.0)*(vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+1] - + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-1]) - -Scalar( 6.0)*(vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+2] - + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-2]) - + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z+3] - + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z-3])}; + return (Scalar){Scalar(1.0 / 60.0) * inv_ds * + (-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] + + Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 1] + + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 1]) - + Scalar(6.0) * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 2] + + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 2]) + + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 3] + + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 3])}; } #endif @@ -68,9 +66,10 @@ hessian(in ScalarField vertex) { Matrix hessian; - hessian.row[0] = (Vector){derxx(vertexIdx, vertex), derxy(vertexIdx, vertex), derxz(vertexIdx, vertex)}; - hessian.row[1] = (Vector){hessian.row[0].y, deryy(vertexIdx, vertex), deryz(vertexIdx, vertex)}; - hessian.row[2] = (Vector){hessian.row[0].z, hessian.row[1].z, derzz(vertexIdx, vertex)}; + hessian.row[0] = (Vector){derxx(vertexIdx, vertex), derxy(vertexIdx, vertex), + derxz(vertexIdx, vertex)}; + hessian.row[1] = (Vector){hessian.row[0].y, deryy(vertexIdx, vertex), deryz(vertexIdx, vertex)}; + hessian.row[2] = (Vector){hessian.row[0].z, hessian.row[1].z, derzz(vertexIdx, vertex)}; return hessian; } diff --git a/acc/mhd_solver/stencil_process.sps b/acc/mhd_solver/stencil_process.sps index 420193e..e605b9f 100644 --- a/acc/mhd_solver/stencil_process.sps +++ b/acc/mhd_solver/stencil_process.sps @@ -1,6 +1,5 @@ #include "stencil_definition.sdh" - Vector value(in VectorField uu) { @@ -14,7 +13,7 @@ upwd_der6(in VectorField uu, in ScalarField lnrho) Scalar uux = fabs(value(uu).x); Scalar uuy = fabs(value(uu).y); Scalar uuz = fabs(value(uu).z); - return (Scalar){uux*der6x_upwd(lnrho) + uuy*der6y_upwd(lnrho) + uuz*der6z_upwd(lnrho)}; + return (Scalar){uux * der6x_upwd(lnrho) + uuy * der6y_upwd(lnrho) + uuz * der6z_upwd(lnrho)}; } #endif @@ -25,10 +24,11 @@ gradients(in VectorField uu) } Scalar -continuity(in VectorField uu, in ScalarField lnrho) { +continuity(in VectorField uu, in ScalarField lnrho) +{ return -dot(value(uu), gradient(lnrho)) #if LUPWD - //This is a corrective hyperdiffusion term for upwinding. + // This is a corrective hyperdiffusion term for upwinding. + upwd_der6(uu, lnrho) #endif - divergence(uu); @@ -36,133 +36,136 @@ continuity(in VectorField uu, in ScalarField lnrho) { #if LENTROPY Vector -momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa) { - const Matrix S = stress_tensor(uu); - const Scalar cs2 = AC_cs2_sound * exp(AC_gamma * value(ss) / AC_cp_sound + (AC_gamma - 1) * (value(lnrho) - AC_lnrho0)); - const Vector j = (Scalar(1.) / AC_mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density +momentum(in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa) +{ + const Matrix S = stress_tensor(uu); + const Scalar cs2 = AC_cs2_sound * exp(AC_gamma * value(ss) / AC_cp_sound + + (AC_gamma - 1) * (value(lnrho) - AC_lnrho0)); + const Vector j = (Scalar(1.) / AC_mu0) * + (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density const Vector B = curl(aa); - //TODO: DOES INTHERMAL VERSTION INCLUDE THE MAGNETIC FIELD? + // TODO: DOES INTHERMAL VERSTION INCLUDE THE MAGNETIC FIELD? const Scalar inv_rho = Scalar(1.) / exp(value(lnrho)); // Regex replace CPU constants with get\(AC_([a-zA-Z_0-9]*)\) // \1 - const Vector mom = - mul(gradients(uu), value(uu)) - - cs2 * ((Scalar(1.) / AC_cp_sound) * gradient(ss) + gradient(lnrho)) - + inv_rho * cross(j, B) - + AC_nu_visc * ( - laplace_vec(uu) - + Scalar(1. / 3.) * gradient_of_divergence(uu) - + Scalar(2.) * mul(S, gradient(lnrho)) - ) - + AC_zeta * gradient_of_divergence(uu); + const Vector mom = -mul(gradients(uu), value(uu)) - + cs2 * ((Scalar(1.) / AC_cp_sound) * gradient(ss) + gradient(lnrho)) + + inv_rho * cross(j, B) + + AC_nu_visc * + (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) + + Scalar(2.) * mul(S, gradient(lnrho))) + + AC_zeta * gradient_of_divergence(uu); return mom; } #elif LTEMPERATURE Vector -momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt) { - Vector mom; +momentum(in VectorField uu, in ScalarField lnrho, in ScalarField tt) +{ + Vector mom; - const Matrix S = stress_tensor(uu); + const Matrix S = stress_tensor(uu); - const Vector pressure_term = (AC_cp_sound - AC_cv_sound) * (gradient(tt) + value(tt) * gradient(lnrho)); + const Vector pressure_term = (AC_cp_sound - AC_cv_sound) * + (gradient(tt) + value(tt) * gradient(lnrho)); - mom = -mul(gradients(uu), value(uu)) - - pressure_term + - AC_nu_visc * - (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) + - Scalar(2.) * mul(S, gradient(lnrho))) + AC_zeta * gradient_of_divergence(uu); + mom = -mul(gradients(uu), value(uu)) - pressure_term + + AC_nu_visc * (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) + + Scalar(2.) * mul(S, gradient(lnrho))) + + AC_zeta * gradient_of_divergence(uu); - #if LGRAVITY - mom = mom - (Vector){0, 0, -10.0}; - #endif +#if LGRAVITY + mom = mom - (Vector){0, 0, -10.0}; +#endif - return mom; + return mom; } #else Vector -momentum(in VectorField uu, in ScalarField lnrho) { - Vector mom; +momentum(in VectorField uu, in ScalarField lnrho) +{ + Vector mom; - const Matrix S = stress_tensor(uu); + const Matrix S = stress_tensor(uu); // Isothermal: we have constant speed of sound - mom = -mul(gradients(uu), value(uu)) - - AC_cs2_sound * gradient(lnrho) + - AC_nu_visc * - (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) + - Scalar(2.) * mul(S, gradient(lnrho))) + AC_zeta * gradient_of_divergence(uu); + mom = -mul(gradients(uu), value(uu)) - AC_cs2_sound * gradient(lnrho) + + AC_nu_visc * (laplace_vec(uu) + Scalar(1. / 3.) * gradient_of_divergence(uu) + + Scalar(2.) * mul(S, gradient(lnrho))) + + AC_zeta * gradient_of_divergence(uu); - #if LGRAVITY - mom = mom - (Vector){0, 0, -10.0}; - #endif +#if LGRAVITY + mom = mom - (Vector){0, 0, -10.0}; +#endif - return mom; + return mom; } #endif - Vector -induction(in VectorField uu, in VectorField aa) { - // Note: We do (-nabla^2 A + nabla(nabla dot A)) instead of (nabla x (nabla - // x A)) in order to avoid taking the first derivative twice (did the math, - // yes this actually works. See pg.28 in arXiv:astro-ph/0109497) - // u cross B - AC_eta * AC_mu0 * (AC_mu0^-1 * [- laplace A + grad div A ]) - const Vector B = curl(aa); - const Vector grad_div = gradient_of_divergence(aa); - const Vector lap = laplace_vec(aa); +induction(in VectorField uu, in VectorField aa) +{ + // Note: We do (-nabla^2 A + nabla(nabla dot A)) instead of (nabla x (nabla + // x A)) in order to avoid taking the first derivative twice (did the math, + // yes this actually works. See pg.28 in arXiv:astro-ph/0109497) + // u cross B - AC_eta * AC_mu0 * (AC_mu0^-1 * [- laplace A + grad div A ]) + const Vector B = curl(aa); + const Vector grad_div = gradient_of_divergence(aa); + const Vector lap = laplace_vec(aa); - // Note, AC_mu0 is cancelled out - const Vector ind = cross(value(uu), B) - AC_eta * (grad_div - lap); + // Note, AC_mu0 is cancelled out + const Vector ind = cross(value(uu), B) - AC_eta * (grad_div - lap); - return ind; + return ind; } - #if LENTROPY Scalar -lnT( in ScalarField ss, in ScalarField lnrho) { - const Scalar lnT = AC_lnT0 + AC_gamma * value(ss) / AC_cp_sound + - (AC_gamma - Scalar(1.)) * (value(lnrho) - AC_lnrho0); - return lnT; +lnT(in ScalarField ss, in ScalarField lnrho) +{ + const Scalar lnT = AC_lnT0 + AC_gamma * value(ss) / AC_cp_sound + + (AC_gamma - Scalar(1.)) * (value(lnrho) - AC_lnrho0); + return lnT; } // Nabla dot (K nabla T) / (rho T) Scalar -heat_conduction( in ScalarField ss, in ScalarField lnrho) { - const Scalar inv_AC_cp_sound = AcReal(1.) / AC_cp_sound; +heat_conduction(in ScalarField ss, in ScalarField lnrho) +{ + const Scalar inv_AC_cp_sound = AcReal(1.) / AC_cp_sound; - const Vector grad_ln_chi = - gradient(lnrho); + const Vector grad_ln_chi = -gradient(lnrho); - const Scalar first_term = AC_gamma * inv_AC_cp_sound * laplace(ss) + - (AC_gamma - AcReal(1.)) * laplace(lnrho); - const Vector second_term = AC_gamma * inv_AC_cp_sound * gradient(ss) + - (AC_gamma - AcReal(1.)) * gradient(lnrho); - const Vector third_term = AC_gamma * (inv_AC_cp_sound * gradient(ss) + - gradient(lnrho)) + grad_ln_chi; + const Scalar first_term = AC_gamma * inv_AC_cp_sound * laplace(ss) + + (AC_gamma - AcReal(1.)) * laplace(lnrho); + const Vector second_term = AC_gamma * inv_AC_cp_sound * gradient(ss) + + (AC_gamma - AcReal(1.)) * gradient(lnrho); + const Vector third_term = AC_gamma * (inv_AC_cp_sound * gradient(ss) + gradient(lnrho)) + + grad_ln_chi; - const Scalar chi = AC_THERMAL_CONDUCTIVITY / (exp(value(lnrho)) * AC_cp_sound); - return AC_cp_sound * chi * (first_term + dot(second_term, third_term)); + const Scalar chi = AC_THERMAL_CONDUCTIVITY / (exp(value(lnrho)) * AC_cp_sound); + return AC_cp_sound * chi * (first_term + dot(second_term, third_term)); } Scalar -heating(const int i, const int j, const int k) { - return 1; +heating(const int i, const int j, const int k) +{ + return 1; } Scalar -entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa) { - const Matrix S = stress_tensor(uu); +entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa) +{ + const Matrix S = stress_tensor(uu); const Scalar inv_pT = Scalar(1.) / (exp(value(lnrho)) * exp(lnT(ss, lnrho))); - const Vector j = (Scalar(1.) / AC_mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density - const Scalar RHS = H_CONST - C_CONST - + AC_eta * (AC_mu0) * dot(j, j) - + Scalar(2.) * exp(value(lnrho)) * AC_nu_visc * contract(S) - + AC_zeta * exp(value(lnrho)) * divergence(uu) * divergence(uu); + const Vector j = (Scalar(1.) / AC_mu0) * + (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density + const Scalar RHS = H_CONST - C_CONST + AC_eta * (AC_mu0)*dot(j, j) + + Scalar(2.) * exp(value(lnrho)) * AC_nu_visc * contract(S) + + AC_zeta * exp(value(lnrho)) * divergence(uu) * divergence(uu); - return - dot(value(uu), gradient(ss)) - + inv_pT * RHS - + heat_conduction(ss, lnrho); + return -dot(value(uu), gradient(ss)) + inv_pT * RHS + heat_conduction(ss, lnrho); } #endif @@ -170,14 +173,15 @@ entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorFie Scalar heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt) { - const Matrix S = stress_tensor(uu); - const Scalar heat_diffusivity_k = 0.0008; //8e-4; - return -dot(value(uu), gradient(tt)) + heat_diffusivity_k * laplace(tt) + heat_diffusivity_k * dot(gradient(lnrho), gradient(tt)) + AC_nu_visc * contract(S) * (Scalar(1.) / AC_cv_sound) - (AC_gamma - 1) * value(tt) * divergence(uu); + const Matrix S = stress_tensor(uu); + const Scalar heat_diffusivity_k = 0.0008; // 8e-4; + return -dot(value(uu), gradient(tt)) + heat_diffusivity_k * laplace(tt) + + heat_diffusivity_k * dot(gradient(lnrho), gradient(tt)) + + AC_nu_visc * contract(S) * (Scalar(1.) / AC_cv_sound) - + (AC_gamma - 1) * value(tt) * divergence(uu); } #endif - - #if LFORCING Vector simple_vortex_forcing(Vector a, Vector b, Scalar magnitude) @@ -191,8 +195,8 @@ simple_outward_flow_forcing(Vector a, Vector b, Scalar 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 Vector helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi) { @@ -206,24 +210,23 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto // MV: Good idea. No an immediate priority. // Fun related article: // https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/ - xx.x = xx.x*(2.0*M_PI/(AC_dsx*globalGridN.x)); - xx.y = xx.y*(2.0*M_PI/(AC_dsy*globalGridN.y)); - xx.z = xx.z*(2.0*M_PI/(AC_dsz*globalGridN.z)); + xx.x = xx.x * (2.0 * M_PI / (AC_dsx * globalGridN.x)); + xx.y = xx.y * (2.0 * M_PI / (AC_dsy * globalGridN.y)); + xx.z = xx.z * (2.0 * M_PI / (AC_dsz * globalGridN.z)); - Scalar cos_phi = cos(phi); - Scalar sin_phi = sin(phi); - Scalar cos_k_dot_x = cos(dot(k_force, xx)); - Scalar sin_k_dot_x = sin(dot(k_force, xx)); + Scalar cos_phi = cos(phi); + Scalar sin_phi = sin(phi); + Scalar cos_k_dot_x = cos(dot(k_force, xx)); + Scalar sin_k_dot_x = sin(dot(k_force, xx)); // Phase affect only the x-component - //Scalar real_comp = cos_k_dot_x; - //Scalar imag_comp = sin_k_dot_x; - Scalar real_comp_phase = cos_k_dot_x*cos_phi - sin_k_dot_x*sin_phi; - Scalar imag_comp_phase = cos_k_dot_x*sin_phi + sin_k_dot_x*cos_phi; + // Scalar real_comp = cos_k_dot_x; + // Scalar imag_comp = sin_k_dot_x; + Scalar real_comp_phase = cos_k_dot_x * cos_phi - sin_k_dot_x * sin_phi; + Scalar imag_comp_phase = cos_k_dot_x * sin_phi + sin_k_dot_x * cos_phi; - - Vector force = (Vector){ 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}; + Vector force = (Vector){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; } @@ -231,37 +234,39 @@ helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vecto Vector forcing(int3 globalVertexIdx, Scalar dt) { - Vector a = Scalar(.5) * (Vector){globalGridN.x * AC_dsx, - globalGridN.y * AC_dsy, + Vector a = Scalar(.5) * (Vector){globalGridN.x * AC_dsx, globalGridN.y * AC_dsy, globalGridN.z * AC_dsz}; // source (origin) - Vector xx = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx, - (globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy, - (globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz}; // sink (current index) + Vector xx = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx, + (globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy, + (globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz}; // sink (current index) const Scalar cs2 = AC_cs2_sound; - const Scalar cs = sqrt(cs2); + const Scalar cs = sqrt(cs2); - //Placeholders until determined properly + // Placeholders until determined properly Scalar magnitude = AC_forcing_magnitude; Scalar phase = AC_forcing_phase; - Vector k_force = (Vector){ AC_k_forcex, AC_k_forcey, AC_k_forcez}; + Vector k_force = (Vector){AC_k_forcex, AC_k_forcey, AC_k_forcez}; Vector ff_re = (Vector){AC_ff_hel_rex, AC_ff_hel_rey, AC_ff_hel_rez}; Vector ff_im = (Vector){AC_ff_hel_imx, AC_ff_hel_imy, AC_ff_hel_imz}; + // Determine that forcing funtion type at this point. + // Vector force = simple_vortex_forcing(a, xx, magnitude); + // Vector force = simple_outward_flow_forcing(a, xx, magnitude); + Vector force = helical_forcing(magnitude, k_force, xx, ff_re, ff_im, phase); - //Determine that forcing funtion type at this point. - //Vector force = simple_vortex_forcing(a, xx, magnitude); - //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) * dt + const Scalar NN = cs * sqrt(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; - //Scaling N = magnitude*cs*sqrt(k*cs/dt) * dt - const Scalar NN = cs*sqrt(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 (Vector){0, 0, 0}; } + if (is_valid(force)) { + return force; + } + else { + return (Vector){0, 0, 0}; + } } #endif // LFORCING @@ -271,12 +276,11 @@ in ScalarField lnrho(VTXBUF_LNRHO); out ScalarField out_lnrho(VTXBUF_LNRHO); in VectorField uu(VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); -out VectorField out_uu(VTXBUF_UUX,VTXBUF_UUY,VTXBUF_UUZ); - +out VectorField out_uu(VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ); #if LMAGNETIC -in VectorField aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ); -out VectorField out_aa(VTXBUF_AX,VTXBUF_AY,VTXBUF_AZ); +in VectorField aa(VTXBUF_AX, VTXBUF_AY, VTXBUF_AZ); +out VectorField out_aa(VTXBUF_AX, VTXBUF_AY, VTXBUF_AZ); #endif #if LENTROPY @@ -290,26 +294,27 @@ out ScalarField out_tt(VTXBUF_TEMPERATURE); #endif Kernel void -solve(Scalar dt) { +solve(Scalar dt) +{ out_lnrho = rk3(out_lnrho, lnrho, continuity(uu, lnrho), dt); - #if LMAGNETIC +#if LMAGNETIC out_aa = rk3(out_aa, aa, induction(uu, aa), dt); - #endif +#endif - #if LENTROPY - out_uu = rk3(out_uu, uu, momentum(uu, lnrho, ss, aa), dt); - out_ss = rk3(out_ss, ss, entropy(ss, uu, lnrho, aa), dt); - #elif LTEMPERATURE - out_uu = rk3(out_uu, uu, momentum(uu, lnrho, tt), dt); - out_tt = rk3(out_tt, tt, heat_transfer(uu, lnrho, tt), dt); - #else - out_uu = rk3(out_uu, uu, momentum(uu, lnrho), dt); - #endif +#if LENTROPY + out_uu = rk3(out_uu, uu, momentum(uu, lnrho, ss, aa), dt); + out_ss = rk3(out_ss, ss, entropy(ss, uu, lnrho, aa), dt); +#elif LTEMPERATURE + out_uu = rk3(out_uu, uu, momentum(uu, lnrho, tt), dt); + out_tt = rk3(out_tt, tt, heat_transfer(uu, lnrho, tt), dt); +#else + out_uu = rk3(out_uu, uu, momentum(uu, lnrho), dt); +#endif - #if LFORCING +#if LFORCING if (step_number == 2) { out_uu = out_uu + forcing(globalVertexIdx, dt); } - #endif +#endif }