From 3fe7b62d3edfe75b4fa126db3c5512dffd408a84 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Mon, 7 Oct 2019 17:37:09 +0300 Subject: [PATCH] Removed the old accrevision directory --- acc/accrevision/stencil_kernel.ac | 680 ------------------------------ 1 file changed, 680 deletions(-) delete mode 100644 acc/accrevision/stencil_kernel.ac diff --git a/acc/accrevision/stencil_kernel.ac b/acc/accrevision/stencil_kernel.ac deleted file mode 100644 index abe1053..0000000 --- a/acc/accrevision/stencil_kernel.ac +++ /dev/null @@ -1,680 +0,0 @@ -#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 -}