diff --git a/acc/mhd_solver/stencil_definition.h b/acc/mhd_solver/stencil_definition.h index 694c043..83d240b 100644 --- a/acc/mhd_solver/stencil_definition.h +++ b/acc/mhd_solver/stencil_definition.h @@ -9,8 +9,8 @@ #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 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; diff --git a/acc/mhd_solver/stencil_process.ac b/acc/mhd_solver/stencil_process.ac index 2c9dc44..f2b96cd 100644 --- a/acc/mhd_solver/stencil_process.ac +++ b/acc/mhd_solver/stencil_process.ac @@ -1,7 +1,7 @@ #include -#include "stencil_definition.h" #include "stencil_assembly.h" +#include "stencil_definition.h" #if LUPWD Device Scalar @@ -22,42 +22,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 - 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)); +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; @@ -66,39 +68,41 @@ 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, +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; @@ -106,50 +110,50 @@ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){ // 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, +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 continuity(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt) { @@ -159,16 +163,15 @@ 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 + @@ -188,16 +191,17 @@ 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 @@ -215,11 +219,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}; @@ -240,15 +244,16 @@ 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}; @@ -339,22 +344,26 @@ heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt) #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}; } } @@ -398,41 +407,45 @@ 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); - //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 @@ -492,10 +505,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 }