From ff12332f06823e2898b49cfdf100e3e52f32a968 Mon Sep 17 00:00:00 2001 From: jpekkila Date: Mon, 7 Oct 2019 18:24:32 +0300 Subject: [PATCH] Clarified the syntax for real number literals. 1.0 is the same precision as AcReal, 1.0f is an explicit float and 1.0d is an explicit double. --- acc/mhd_solver/stencil_kernel.ac | 331 ++++++++++++++++--------------- acc/src/acc.l | 3 +- acc/src/acc.y | 3 +- acc/stdlib/stdderiv.h | 18 +- 4 files changed, 184 insertions(+), 171 deletions(-) diff --git a/acc/mhd_solver/stencil_kernel.ac b/acc/mhd_solver/stencil_kernel.ac index 12774b3..1fe57a9 100644 --- a/acc/mhd_solver/stencil_kernel.ac +++ b/acc/mhd_solver/stencil_kernel.ac @@ -9,9 +9,9 @@ #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 +#define AC_THERMAL_CONDUCTIVITY (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; @@ -140,12 +140,12 @@ 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]) + + return (Scalar){(1.0 / 60.0) * inv_ds * + (-20.0 * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] + + 15.0 * (vertex[vertexIdx.x + 1, vertexIdx.y, vertexIdx.z] + + vertex[vertexIdx.x - 1, vertexIdx.y, vertexIdx.z]) - + 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])}; } @@ -155,12 +155,12 @@ 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]) + + return (Scalar){(1.0 / 60.0) * inv_ds * + (-20.0 * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] + + 15.0 * (vertex[vertexIdx.x, vertexIdx.y + 1, vertexIdx.z] + + vertex[vertexIdx.x, vertexIdx.y - 1, vertexIdx.z]) - + 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])}; } @@ -170,19 +170,18 @@ 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]) + + return (Scalar){(1.0 / 60.0) * inv_ds * + (-20.0 * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] + + 15.0 * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 1] + + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 1]) - + 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) @@ -202,42 +201,44 @@ gradients(in VectorField uu) #if LSINK Device Vector -sink_gravity(int3 globalVertexIdx){ +sink_gravity(int3 globalVertexIdx) +{ int accretion_switch = int(AC_switch_accretion); - if (accretion_switch == 1){ + if (accretion_switch == 1) { Vector force_gravity; - const Vector grid_pos = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx, + const Vector grid_pos = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx, (globalVertexIdx.y - AC_ny_min) * AC_dsy, (globalVertexIdx.z - 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 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; + force_gravity = gravity_magnitude * direction; return force_gravity; - } else { + } + 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)); +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. + 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; @@ -246,90 +247,92 @@ truelove_density(in ScalarField lnrho){ // 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 - AC_nx_min) * AC_dsx, +sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt) +{ + const Vector grid_pos = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx, (globalVertexIdx.y - AC_ny_min) * AC_dsy, (globalVertexIdx.z - 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 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; + 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); + if (accretion_switch == 1) { + if ((accretion_distance) <= profile_range) { + // weight = 1.0; + // Hann window function + Scalar window_ratio = accretion_distance / profile_range; + weight = 0.5 * (1.0 - cos(2.0 * M_PI * window_ratio)); + } + else { + weight = 0.0; } - //Truelove criterion is used as a kind of arbitrary density floor. + // 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); + else { + rate = 0.0; + } + accretion_density = weight * rate; + } + else { + accretion_density = 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 - AC_nx_min) * AC_dsx, +sink_accretion_velocity(int3 globalVertexIdx, in VectorField uu, Scalar dt) +{ + const Vector grid_pos = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx, (globalVertexIdx.y - AC_ny_min) * AC_dsy, (globalVertexIdx.z - 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 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; + int accretion_switch = AC_switch_accretion; Vector accretion_velocity; - if (accretion_switch == 1){ + 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); + if ((accretion_distance) <= profile_range) { + // weight = 1.0; + // Hann window function + Scalar window_ratio = accretion_distance / profile_range; + weight = 0.5 * (1.0 - cos(2.0 * M_PI * window_ratio)); + } + else { + weight = 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 { + if (length(value(uu)) > 0.0) { + rate = (1.0 / dt) * value(uu); + } + else { rate = (Vector){0.0, 0.0, 0.0}; } - accretion_velocity = weight * rate ; - } else { + 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) { @@ -339,45 +342,44 @@ continuity(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar + upwd_der6(uu, lnrho) #endif #if LSINK - - sink_accretion(globalVertexIdx, lnrho, dt) / exp(value(lnrho)) + - 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) +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) * + const Vector j = (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)); + const Scalar inv_rho = 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)) + + cs2 * ((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_nu_visc * (laplace_vec(uu) + (1.0 / 3.0) * gradient_of_divergence(uu) + + 2.0 * mul(S, gradient(lnrho))) + AC_zeta * gradient_of_divergence(uu) - #if LSINK - //Gravity term +#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 + // Corresponding loss of momentum + - //(1.0 / ( (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 @@ -392,14 +394,14 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in Scala (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_nu_visc * (laplace_vec(uu) + (1.0 / 3.0) * gradient_of_divergence(uu) + + 2.0 * mul(S, gradient(lnrho))) + AC_zeta * gradient_of_divergence(uu) - #if LSINK +#if LSINK + sink_gravity(globalVertexIdx); - #else - ; - #endif +#else + ; +#endif #if LGRAVITY mom = mom - (Vector){0, 0, -10.0}; @@ -417,18 +419,19 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar d // 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_nu_visc * (laplace_vec(uu) + (1.0 / 3.0) * gradient_of_divergence(uu) + + 2.0 * mul(S, gradient(lnrho))) + AC_zeta * gradient_of_divergence(uu) - #if LSINK +#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 + // Corresponding loss of momentum + - //(1.0 / ( (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}; @@ -460,21 +463,21 @@ 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); + (AC_gamma - 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 Scalar inv_AC_cp_sound = 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); + (AC_gamma - 1.0) * laplace(lnrho); const Vector second_term = AC_gamma * inv_AC_cp_sound * gradient(ss) + - (AC_gamma - AcReal(1.0)) * gradient(lnrho); + (AC_gamma - 1.0) * gradient(lnrho); const Vector third_term = AC_gamma * (inv_AC_cp_sound * gradient(ss) + gradient(lnrho)) + grad_ln_chi; @@ -492,11 +495,11 @@ 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) * + const Scalar inv_pT = 1.0 / (exp(value(lnrho)) * exp(lnT(ss, lnrho))); + const Vector j = (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) + + 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); @@ -511,29 +514,33 @@ heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt) 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_nu_visc * contract(S) * (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){ +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}; + 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){ +simple_outward_flow_forcing(Vector a, Vector b, Scalar magnitude) +{ int accretion_switch = AC_switch_accretion; - if (accretion_switch == 0){ + if (accretion_switch == 0) { return magnitude * (1 / length(b - a)) * normalized(b - a); // Outward flow - } else { - return (Vector){0,0,0}; + } + else { + return (Vector){0, 0, 0}; } } @@ -577,41 +584,44 @@ Device Vector forcing(int3 globalVertexIdx, Scalar dt) { int accretion_switch = AC_switch_accretion; - if (accretion_switch == 0){ + 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 - AC_nx_min) * AC_dsx, + Vector a = 0.5 * (Vector){globalGridN.x * AC_dsx, globalGridN.y * AC_dsy, + globalGridN.z * AC_dsz}; // source (origin) + Vector xx = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx, (globalVertexIdx.y - AC_ny_min) * AC_dsy, (globalVertexIdx.z - 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}; } - } else { - return (Vector){0,0,0}; + if (is_valid(force)) { + return force; + } + else { + return (Vector){0, 0, 0}; + } + } + else { + return (Vector){0, 0, 0}; } } #endif // LFORCING @@ -671,10 +681,11 @@ solve() #endif #if LSINK - out_accretion = rk3(out_accretion, accretion, sink_accretion(globalVertexIdx, lnrho, dt), dt);// unit now is rho! + 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! + out_accretion = out_accretion * AC_dsx * AC_dsy * AC_dsz; // unit is now mass! } #endif } diff --git a/acc/src/acc.l b/acc/src/acc.l index 92d51ad..c2628f8 100644 --- a/acc/src/acc.l +++ b/acc/src/acc.l @@ -39,7 +39,8 @@ L [a-zA-Z_] "return" { return RETURN; } {D}+"."{D}+ { return REAL_NUMBER; } /* Literals */ -{D}+"."{D}+[fd]+ { return NUMBER; } +{D}+"."{D}+d+ { return DOUBLE_NUMBER; } +{D}+"."{D}+f+ { return NUMBER; } {D}+[lu]* { return NUMBER; } {L}({L}|{D})* { return IDENTIFIER; } \"(.)*\" { return IDENTIFIER; } /* String */ diff --git a/acc/src/acc.y b/acc/src/acc.y index 7c41f02..7bed905 100644 --- a/acc/src/acc.y +++ b/acc/src/acc.y @@ -14,7 +14,7 @@ int yyget_lineno(); %} %token CONSTANT IN OUT UNIFORM -%token IDENTIFIER NUMBER REAL_NUMBER +%token IDENTIFIER NUMBER REAL_NUMBER DOUBLE_NUMBER %token RETURN %token SCALAR VECTOR MATRIX SCALARFIELD SCALARARRAY %token VOID INT INT3 COMPLEX @@ -222,6 +222,7 @@ identifier: IDENTIFIER ; number: REAL_NUMBER { $$ = astnode_create(NODE_REAL_NUMBER, NULL, NULL); astnode_set_buffer(yytext, $$); } + | DOUBLE_NUMBER { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); astnode_set_buffer(yytext, $$); $$->buffer[strlen($$->buffer) - 1] = '\0'; } | NUMBER { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); astnode_set_buffer(yytext, $$); } ; diff --git a/acc/stdlib/stdderiv.h b/acc/stdlib/stdderiv.h index 4dc20ea..f7e34fd 100644 --- a/acc/stdlib/stdderiv.h +++ b/acc/stdlib/stdderiv.h @@ -302,18 +302,18 @@ stress_tensor(in VectorField vec) { Matrix S; - S.row[0].x = Scalar(2.0 / 3.0) * gradient(vec.x).x - - Scalar(1.0 / 3.0) * (gradient(vec.y).y + gradient(vec.z).z); - S.row[0].y = Scalar(1.0 / 2.0) * (gradient(vec.x).y + gradient(vec.y).x); - S.row[0].z = Scalar(1.0 / 2.0) * (gradient(vec.x).z + gradient(vec.z).x); + S.row[0].x = (2.0 / 3.0) * gradient(vec.x).x - + (1.0 / 3.0) * (gradient(vec.y).y + gradient(vec.z).z); + S.row[0].y = (1.0 / 2.0) * (gradient(vec.x).y + gradient(vec.y).x); + S.row[0].z = (1.0 / 2.0) * (gradient(vec.x).z + gradient(vec.z).x); - S.row[1].y = Scalar(2.0 / 3.0) * gradient(vec.y).y - - Scalar(1.0 / 3.0) * (gradient(vec.x).x + gradient(vec.z).z); + S.row[1].y = (2.0 / 3.0) * gradient(vec.y).y - + (1.0 / 3.0) * (gradient(vec.x).x + gradient(vec.z).z); - S.row[1].z = Scalar(1.0 / 2.0) * (gradient(vec.y).z + gradient(vec.z).y); + S.row[1].z = (1.0 / 2.0) * (gradient(vec.y).z + gradient(vec.z).y); - S.row[2].z = Scalar(2.0 / 3.0) * gradient(vec.z).z - - Scalar(1.0 / 3.0) * (gradient(vec.x).x + gradient(vec.y).y); + S.row[2].z = (2.0 / 3.0) * gradient(vec.z).z - + (1.0 / 3.0) * (gradient(vec.x).x + gradient(vec.y).y); S.row[1].x = S.row[0].y; S.row[2].x = S.row[0].z;