#include #define LDENSITY (1) #define LHYDRO (1) #define LMAGNETIC (1) #define LENTROPY (1) #define LTEMPERATURE (0) #define LFORCING (1) #define LUPWD (1) #define LSINK (0) #define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // TODO: make an actual config parameter #define H_CONST (0) // TODO: make an actual config parameter #define C_CONST (0) // TODO: make an actual config parameter // Int params uniform int AC_max_steps; uniform int AC_save_steps; uniform int AC_bin_steps; uniform int AC_bc_type; uniform int AC_start_step; // Real params uniform Scalar AC_dt; uniform Scalar AC_max_time; // Spacing uniform Scalar AC_dsmin; // physical grid uniform Scalar AC_xlen; uniform Scalar AC_ylen; uniform Scalar AC_zlen; uniform Scalar AC_xorig; uniform Scalar AC_yorig; uniform Scalar AC_zorig; // Physical units uniform Scalar AC_unit_density; uniform Scalar AC_unit_velocity; uniform Scalar AC_unit_length; // properties of gravitating star uniform Scalar AC_star_pos_x; uniform Scalar AC_star_pos_y; uniform Scalar AC_star_pos_z; uniform Scalar AC_M_star; // properties of sink particle uniform Scalar AC_sink_pos_x; uniform Scalar AC_sink_pos_y; uniform Scalar AC_sink_pos_z; uniform Scalar AC_M_sink; uniform Scalar AC_M_sink_init; uniform Scalar AC_M_sink_Msun; uniform Scalar AC_soft; uniform Scalar AC_accretion_range; uniform Scalar AC_switch_accretion; // Run params uniform Scalar AC_cdt; uniform Scalar AC_cdtv; uniform Scalar AC_cdts; uniform Scalar AC_nu_visc; uniform Scalar AC_cs_sound; uniform Scalar AC_eta; uniform Scalar AC_mu0; uniform Scalar AC_cp_sound; uniform Scalar AC_gamma; uniform Scalar AC_cv_sound; uniform Scalar AC_lnT0; uniform Scalar AC_lnrho0; uniform Scalar AC_zeta; uniform Scalar AC_trans; // Other uniform Scalar AC_bin_save_t; // Initial condition params uniform Scalar AC_ampl_lnrho; uniform Scalar AC_ampl_uu; uniform Scalar AC_angl_uu; uniform Scalar AC_lnrho_edge; uniform Scalar AC_lnrho_out; // Forcing parameters. User configured. uniform Scalar AC_forcing_magnitude; uniform Scalar AC_relhel; uniform Scalar AC_kmin; uniform Scalar AC_kmax; // Forcing parameters. Set by the generator. uniform Scalar AC_forcing_phase; uniform Scalar AC_k_forcex; uniform Scalar AC_k_forcey; uniform Scalar AC_k_forcez; uniform Scalar AC_kaver; uniform Scalar AC_ff_hel_rex; uniform Scalar AC_ff_hel_rey; uniform Scalar AC_ff_hel_rez; uniform Scalar AC_ff_hel_imx; uniform Scalar AC_ff_hel_imy; uniform Scalar AC_ff_hel_imz; // Additional helper params // (deduced from other params do not set these directly!) uniform Scalar AC_G_const; uniform Scalar AC_GM_star; uniform Scalar AC_unit_mass; uniform Scalar AC_sq2GM_star; uniform Scalar AC_cs2_sound; /* * ============================================================================= * User-defined vertex buffers * ============================================================================= */ #if LENTROPY uniform ScalarField VTXBUF_LNRHO; uniform ScalarField VTXBUF_UUX; uniform ScalarField VTXBUF_UUY; uniform ScalarField VTXBUF_UUZ; uniform ScalarField VTXBUF_AX; uniform ScalarField VTXBUF_AY; uniform ScalarField VTXBUF_AZ; uniform ScalarField VTXBUF_ENTROPY; #elif LMAGNETIC uniform ScalarField VTXBUF_LNRHO; uniform ScalarField VTXBUF_UUX; uniform ScalarField VTXBUF_UUY; uniform ScalarField VTXBUF_UUZ; uniform ScalarField VTXBUF_AX; uniform ScalarField VTXBUF_AY; uniform ScalarField VTXBUF_AZ; #elif LHYDRO uniform ScalarField VTXBUF_LNRHO; uniform ScalarField VTXBUF_UUX; uniform ScalarField VTXBUF_UUY; uniform ScalarField VTXBUF_UUZ; #else uniform ScalarField VTXBUF_LNRHO; #endif #if LSINK uniform ScalarField VTXBUF_ACCRETION; #endif #if LUPWD Preprocessed Scalar 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])}; } Preprocessed Scalar 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])}; } Preprocessed Scalar 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])}; } #endif #if LUPWD Device Scalar 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)}; } #endif Device Matrix gradients(in VectorField uu) { return (Matrix){gradient(uu.x), gradient(uu.y), gradient(uu.z)}; } #if LSINK Device Vector sink_gravity(int3 globalVertexIdx){ int accretion_switch = int(AC_switch_accretion); if (accretion_switch == 1){ Vector force_gravity; const Vector grid_pos = (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}; const Scalar sink_mass = AC_M_sink; const Vector sink_pos = (Vector){AC_sink_pos_x, AC_sink_pos_y, AC_sink_pos_z}; const Scalar distance = length(grid_pos - sink_pos); const Scalar soft = AC_soft; //MV: The commit 083ff59 had AC_G_const defined wrong here in DSL making it exxessively strong. //MV: Scalar gravity_magnitude = ... below is correct! const Scalar gravity_magnitude = (AC_G_const * sink_mass) / pow(((distance * distance) + soft*soft), 1.5); const Vector direction = (Vector){(sink_pos.x - grid_pos.x) / distance, (sink_pos.y - grid_pos.y) / distance, (sink_pos.z - grid_pos.z) / distance}; force_gravity = gravity_magnitude * direction; return force_gravity; } else { return (Vector){0.0, 0.0, 0.0}; } } #endif #if LSINK // Give Truelove density Device Scalar truelove_density(in ScalarField lnrho){ const Scalar rho = exp(value(lnrho)); const Scalar Jeans_length_squared = (M_PI * AC_cs2_sound) / (AC_G_const * rho); const Scalar TJ_rho = ((M_PI) * ((AC_dsx * AC_dsx) / Jeans_length_squared) * AC_cs2_sound) / (AC_G_const * AC_dsx * AC_dsx); //TODO: AC_dsx will cancel out, deal with it later for optimization. Scalar accretion_rho = TJ_rho; return accretion_rho; } // This controls accretion of density/mass to the sink particle. Device Scalar sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){ const Vector grid_pos = (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}; const Vector sink_pos = (Vector){AC_sink_pos_x, AC_sink_pos_y, AC_sink_pos_z}; const Scalar profile_range = AC_accretion_range; const Scalar accretion_distance = length(grid_pos - sink_pos); int accretion_switch = AC_switch_accretion; Scalar accretion_density; Scalar weight; if (accretion_switch == 1){ if ((accretion_distance) <= profile_range){ //weight = Scalar(1.0); //Hann window function Scalar window_ratio = accretion_distance/profile_range; weight = Scalar(0.5)*(Scalar(1.0) - cos(Scalar(2.0)*M_PI*window_ratio)); } else { weight = Scalar(0.0); } //Truelove criterion is used as a kind of arbitrary density floor. const Scalar lnrho_min = log(truelove_density(lnrho)); Scalar rate; if (value(lnrho) > lnrho_min) { rate = (exp(value(lnrho)) - exp(lnrho_min)) / dt; } else { rate = Scalar(0.0); } accretion_density = weight * rate ; } else { accretion_density = Scalar(0.0); } return accretion_density; } // This controls accretion of velocity to the sink particle. Device Vector sink_accretion_velocity(int3 globalVertexIdx, in VectorField uu, Scalar dt) { const Vector grid_pos = (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}; const Vector sink_pos = (Vector){AC_sink_pos_x, AC_sink_pos_y, AC_sink_pos_z}; const Scalar profile_range = AC_accretion_range; const Scalar accretion_distance = length(grid_pos - sink_pos); int accretion_switch = AC_switch_accretion; Vector accretion_velocity; if (accretion_switch == 1){ Scalar weight; // Step function weighting // Arch of a cosine function? // Cubic spline x^3 - x in range [-0.5 , 0.5] if ((accretion_distance) <= profile_range){ //weight = Scalar(1.0); //Hann window function Scalar window_ratio = accretion_distance/profile_range; weight = Scalar(0.5)*(Scalar(1.0) - cos(Scalar(2.0)*M_PI*window_ratio)); } else { weight = Scalar(0.0); } Vector rate; // MV: Could we use divergence here ephasize velocitie which are compressive and // MV: not absorbins stuff that would not be accreted anyway? if (length(value(uu)) > Scalar(0.0)) { rate = (Scalar(1.0)/dt) * value(uu); } else { rate = (Vector){0.0, 0.0, 0.0}; } accretion_velocity = weight * rate ; } else { accretion_velocity = (Vector){0.0, 0.0, 0.0}; } return accretion_velocity; } #endif Device Scalar continuity(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt) { return -dot(value(uu), gradient(lnrho)) #if LUPWD // This is a corrective hyperdiffusion term for upwinding. + upwd_der6(uu, lnrho) #endif #if LSINK - sink_accretion(globalVertexIdx, lnrho, dt) / exp(value(lnrho)) #endif - divergence(uu); } #if LENTROPY Device Vector momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in ScalarField ss, in VectorField aa, Scalar dt) { 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.0) / AC_mu0) * (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density const Vector B = curl(aa); // TODO: DOES INTHERMAL VERSTION INCLUDE THE MAGNETIC FIELD? const Scalar inv_rho = Scalar(1.0) / 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.0) / AC_cp_sound) * gradient(ss) + gradient(lnrho)) + inv_rho * cross(j, B) + AC_nu_visc * (laplace_vec(uu) + Scalar(1.0 / 3.0) * gradient_of_divergence(uu) + Scalar(2.0) * mul(S, gradient(lnrho))) + AC_zeta * gradient_of_divergence(uu) #if LSINK //Gravity term + sink_gravity(globalVertexIdx) //Corresponding loss of momentum - //(Scalar(1.0) / Scalar( (AC_dsx*AC_dsy*AC_dsz) * exp(value(lnrho)))) * // Correction factor by unit mass sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014) ; #else ; #endif return mom; } #elif LTEMPERATURE Device Vector momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in ScalarField tt) { Vector mom; const Matrix S = stress_tensor(uu); 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.0 / 3.0) * gradient_of_divergence(uu) + Scalar(2.0) * mul(S, gradient(lnrho))) + AC_zeta * gradient_of_divergence(uu) #if LSINK + sink_gravity(globalVertexIdx); #else ; #endif #if LGRAVITY mom = mom - (Vector){0, 0, -10.0}; #endif return mom; } #else Device Vector momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt) { Vector mom; 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.0 / 3.0) * gradient_of_divergence(uu) + Scalar(2.0) * mul(S, gradient(lnrho))) + AC_zeta * gradient_of_divergence(uu) #if LSINK + sink_gravity(globalVertexIdx) //Corresponding loss of momentum - //(Scalar(1.0) / Scalar( (AC_dsx*AC_dsy*AC_dsz) * exp(value(lnrho)))) * // Correction factor by unit mass sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014) ; #else ; #endif #if LGRAVITY mom = mom - (Vector){0, 0, -10.0}; #endif return mom; } #endif Device 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); // Note, AC_mu0 is cancelled out const Vector ind = cross(value(uu), B) - AC_eta * (grad_div - lap); return ind; } #if LENTROPY Device Scalar lnT(in ScalarField ss, in ScalarField lnrho) { return AC_lnT0 + AC_gamma * value(ss) / AC_cp_sound + (AC_gamma - Scalar(1.0)) * (value(lnrho) - AC_lnrho0); } // Nabla dot (K nabla T) / (rho T) Device Scalar heat_conduction(in ScalarField ss, in ScalarField lnrho) { const Scalar inv_AC_cp_sound = AcReal(1.0) / AC_cp_sound; const Vector grad_ln_chi = -gradient(lnrho); const Scalar first_term = AC_gamma * inv_AC_cp_sound * laplace(ss) + (AC_gamma - AcReal(1.0)) * laplace(lnrho); const Vector second_term = AC_gamma * inv_AC_cp_sound * gradient(ss) + (AC_gamma - AcReal(1.0)) * 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)); } Device Scalar heating(const int i, const int j, const int k) { return 1; } Device Scalar 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.0) / (exp(value(lnrho)) * exp(lnT(ss, lnrho))); const Vector j = (Scalar(1.0) / 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.0) * 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); } #endif #if LTEMPERATURE Device 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.0) / AC_cv_sound) - (AC_gamma - 1) * value(tt) * divergence(uu); } #endif #if LFORCING Device Vector simple_vortex_forcing(Vector a, Vector b, Scalar magnitude){ int accretion_switch = AC_switch_accretion; if (accretion_switch == 0){ return magnitude * cross(normalized(b - a), (Vector){ 0, 0, 1}); // Vortex } else { return (Vector){0,0,0}; } } Device Vector simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude){ int accretion_switch = AC_switch_accretion; if (accretion_switch == 0){ return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow } else { return (Vector){0,0,0}; } } // The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable // helicity Device Vector helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi) { // JP: This looks wrong: // 1) Should it be AC_dsx * AC_nx instead of AC_dsx * AC_ny? // 2) Should you also use globalGrid.n instead of the local n? // MV: You are rigth. Made a quickfix. I did not see the error because multigpu is split // in z direction not y direction. // 3) Also final point: can we do this with vectors/quaternions instead? // Tringonometric functions are much more expensive and inaccurate/ // 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)); 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; 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; } Device Vector forcing(int3 globalVertexIdx, Scalar dt) { int accretion_switch = AC_switch_accretion; if (accretion_switch == 0){ Vector a = Scalar(0.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) const Scalar cs2 = AC_cs2_sound; const Scalar cs = sqrt(cs2); //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 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); //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}; } } else { return (Vector){0,0,0}; } } #endif // LFORCING // Declare input and output arrays using locations specified in the // array enum in astaroth.h 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); #if LMAGNETIC in VectorField aa(VTXBUF_AX, VTXBUF_AY, VTXBUF_AZ); out VectorField out_aa(VTXBUF_AX, VTXBUF_AY, VTXBUF_AZ); #endif #if LENTROPY in ScalarField ss(VTXBUF_ENTROPY); out ScalarField out_ss(VTXBUF_ENTROPY); #endif #if LTEMPERATURE in ScalarField tt(VTXBUF_TEMPERATURE); out ScalarField out_tt(VTXBUF_TEMPERATURE); #endif #if LSINK in ScalarField accretion(VTXBUF_ACCRETION); out ScalarField out_accretion(VTXBUF_ACCRETION); #endif Kernel void solve() { Scalar dt = AC_dt; out_lnrho = rk3(out_lnrho, lnrho, continuity(globalVertexIdx, uu, lnrho, dt), dt); #if LMAGNETIC out_aa = rk3(out_aa, aa, induction(uu, aa), dt); #endif #if LENTROPY out_uu = rk3(out_uu, uu, momentum(globalVertexIdx, uu, lnrho, ss, aa, dt), dt); out_ss = rk3(out_ss, ss, entropy(ss, uu, lnrho, aa), dt); #elif LTEMPERATURE out_uu = rk3(out_uu, uu, momentum(globalVertexIdx, uu, lnrho, tt, dt), dt); out_tt = rk3(out_tt, tt, heat_transfer(uu, lnrho, tt), dt); #else out_uu = rk3(out_uu, uu, momentum(globalVertexIdx, uu, lnrho, dt), dt); #endif #if LFORCING if (step_number == 2) { out_uu = out_uu + forcing(globalVertexIdx, dt); } #endif #if LSINK out_accretion = rk3(out_accretion, accretion, sink_accretion(globalVertexIdx, lnrho, dt), dt);// unit now is rho! if (step_number == 2) { out_accretion = out_accretion * AC_dsx * AC_dsy * AC_dsz;// unit is now mass! } #endif }