diff --git a/acc/mhd_solver/stencil_assembly.h b/acc/mhd_solver/stencil_assembly.sas similarity index 75% rename from acc/mhd_solver/stencil_assembly.h rename to acc/mhd_solver/stencil_assembly.sas index d1ffb17..85886e6 100644 --- a/acc/mhd_solver/stencil_assembly.h +++ b/acc/mhd_solver/stencil_assembly.sas @@ -1,4 +1,17 @@ -#pragma once +#include "stencil_definition.sdh" + +Preprocessed Scalar +value(in ScalarField vertex) +{ + return vertex[vertexIdx]; +} + +Preprocessed Vector +gradient(in ScalarField vertex) +{ + return (Vector){derx(vertexIdx, vertex), dery(vertexIdx, vertex), derz(vertexIdx, vertex)}; +} + #if LUPWD Preprocessed Scalar @@ -47,3 +60,16 @@ der6z_upwd(in ScalarField vertex) } #endif + +Preprocessed Matrix +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)}; + + return hessian; +} diff --git a/acc/mhd_solver/stencil_definition.h b/acc/mhd_solver/stencil_definition.sdh similarity index 94% rename from acc/mhd_solver/stencil_definition.h rename to acc/mhd_solver/stencil_definition.sdh index 83d240b..47ea8ab 100644 --- a/acc/mhd_solver/stencil_definition.h +++ b/acc/mhd_solver/stencil_definition.sdh @@ -1,4 +1,3 @@ -#pragma once #define LDENSITY (1) #define LHYDRO (1) #define LMAGNETIC (1) @@ -9,8 +8,6 @@ #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; @@ -23,6 +20,9 @@ uniform int AC_start_step; uniform Scalar AC_dt; uniform Scalar AC_max_time; // Spacing +uniform Scalar AC_dsx; +uniform Scalar AC_dsy; +uniform Scalar AC_dsz; uniform Scalar AC_dsmin; // physical grid uniform Scalar AC_xlen; @@ -96,6 +96,9 @@ uniform Scalar AC_GM_star; uniform Scalar AC_unit_mass; uniform Scalar AC_sq2GM_star; uniform Scalar AC_cs2_sound; +uniform Scalar AC_inv_dsx; +uniform Scalar AC_inv_dsy; +uniform Scalar AC_inv_dsz; /* * ============================================================================= @@ -131,3 +134,4 @@ uniform ScalarField VTXBUF_LNRHO; #if LSINK uniform ScalarField VTXBUF_ACCRETION; #endif + diff --git a/acc/mhd_solver/stencil_process.ac b/acc/mhd_solver/stencil_process.sps similarity index 69% rename from acc/mhd_solver/stencil_process.ac rename to acc/mhd_solver/stencil_process.sps index f2b96cd..417c85d 100644 --- a/acc/mhd_solver/stencil_process.ac +++ b/acc/mhd_solver/stencil_process.sps @@ -1,10 +1,13 @@ -#include +#include "stencil_definition.sdh" -#include "stencil_assembly.h" -#include "stencil_definition.h" +Vector +value(in VectorField uu) +{ + return (Vector){value(uu.x), value(uu.y), value(uu.z)}; +} #if LUPWD -Device Scalar +Scalar upwd_der6(in VectorField uu, in ScalarField lnrho) { Scalar uux = fabs(value(uu).x); @@ -14,52 +17,50 @@ upwd_der6(in VectorField uu, in ScalarField lnrho) } #endif -Device Matrix +Matrix gradients(in VectorField uu) { return (Matrix){gradient(uu.x), gradient(uu.y), gradient(uu.z)}; } #if LSINK -Device Vector -sink_gravity(int3 globalVertexIdx) -{ +Vector +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 - DCONST(AC_nx_min)) * AC_dsx, + 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 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)); +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. + 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; @@ -67,94 +68,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 - DCONST(AC_nx_min)) * AC_dsx, +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 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 { + 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. + //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 { + } else { rate = Scalar(0.0); } - accretion_density = weight * rate; - } - else { + 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, +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 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 { + 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 = (Scalar(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 + +Scalar continuity(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt) { return -dot(value(uu), gradient(lnrho)) @@ -163,15 +162,16 @@ 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) +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 + @@ -191,21 +191,20 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in Scala (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 + #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 + - //(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 +Vector momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in ScalarField tt) { Vector mom; @@ -219,11 +218,11 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in Scala 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 + #if LSINK + sink_gravity(globalVertexIdx); -#else - ; -#endif + #else + ; + #endif #if LGRAVITY mom = mom - (Vector){0, 0, -10.0}; @@ -231,7 +230,7 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in Scala return mom; } #else -Device Vector +Vector momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt) { Vector mom; @@ -244,16 +243,15 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar d 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 + #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 + - //(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}; @@ -263,7 +261,7 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar d } #endif -Device Vector +Vector induction(in VectorField uu, in VectorField aa) { // Note: We do (-nabla^2 A + nabla(nabla dot A)) instead of (nabla x (nabla @@ -281,7 +279,7 @@ induction(in VectorField uu, in VectorField aa) } #if LENTROPY -Device Scalar +Scalar lnT(in ScalarField ss, in ScalarField lnrho) { const Scalar lnT = AC_lnT0 + AC_gamma * value(ss) / AC_cp_sound + @@ -290,7 +288,7 @@ lnT(in ScalarField ss, in ScalarField lnrho) } // Nabla dot (K nabla T) / (rho T) -Device Scalar +Scalar heat_conduction(in ScalarField ss, in ScalarField lnrho) { const Scalar inv_AC_cp_sound = AcReal(1.0) / AC_cp_sound; @@ -308,13 +306,13 @@ heat_conduction(in ScalarField ss, in ScalarField lnrho) return AC_cp_sound * chi * (first_term + dot(second_term, third_term)); } -Device Scalar +Scalar heating(const int i, const int j, const int k) { return 1; } -Device Scalar +Scalar entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa) { const Matrix S = stress_tensor(uu); @@ -330,7 +328,7 @@ entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorFie #endif #if LTEMPERATURE -Device Scalar +Scalar heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt) { const Matrix S = stress_tensor(uu); @@ -343,33 +341,29 @@ heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt) #endif #if LFORCING -Device Vector -simple_vortex_forcing(Vector a, Vector b, Scalar magnitude) -{ +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}; + 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) -{ +Vector + 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}; } } // The Pencil Code forcing_hel_noshear(), manual Eq. 222, inspired forcing function with adjustable // helicity -Device Vector +Vector helical_forcing(Scalar magnitude, Vector k_force, Vector xx, Vector ff_re, Vector ff_im, Scalar phi) { // JP: This looks wrong: @@ -407,45 +401,41 @@ 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 - DCONST(AC_nx_min)) * AC_dsx, + 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) + (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); - // 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; + //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); - if (is_valid(force)) { - return force; - } - else { - return (Vector){0, 0, 0}; - } - } - else { - return (Vector){0, 0, 0}; + //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 @@ -505,11 +495,10 @@ 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 }