cwpearson/astaroth
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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.

Browse Source
This commit is contained in:
jpekkila
2019-10-07 18:24:32 +03:00
parent aa6c2b23d9
commit ff12332f06
4 changed files with 184 additions and 171 deletions
Download Patch File Download Diff File Expand all files Collapse all files

331
acc/mhd_solver/stencil_kernel.ac
View File

@@ -9,9 +9,9 @@
#define LUPWD (1) #define LUPWD (1)
#define LSINK (0) #define LSINK (0)
#define AC_THERMAL_CONDUCTIVITY (AcReal(0.001)) // 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 H_CONST (0) // TODO: make an actual config parameter
#define C_CONST (0) // TODO: make an actual config parameter #define C_CONST (0) // TODO: make an actual config parameter
// Int params // Int params
uniform int AC_max_steps; uniform int AC_max_steps;
@@ -140,12 +140,12 @@ der6x_upwd(in ScalarField vertex)
{ {
Scalar inv_ds = AC_inv_dsx; Scalar inv_ds = AC_inv_dsx;
return (Scalar){Scalar(1.0 / 60.0) * inv_ds * return (Scalar){(1.0 / 60.0) * inv_ds *
(-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] + (-20.0 * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
Scalar(15.0) * (vertex[vertexIdx.x + 1, vertexIdx.y, vertexIdx.z] + 15.0 * (vertex[vertexIdx.x + 1, vertexIdx.y, vertexIdx.z] +
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] + 6.0 * (vertex[vertexIdx.x + 2, vertexIdx.y, vertexIdx.z] +
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] +
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; Scalar inv_ds = AC_inv_dsy;
return (Scalar){Scalar(1.0 / 60.0) * inv_ds * return (Scalar){(1.0 / 60.0) * inv_ds *
(-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] + (-20.0 * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y + 1, vertexIdx.z] + 15.0 * (vertex[vertexIdx.x, vertexIdx.y + 1, vertexIdx.z] +
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] + 6.0 * (vertex[vertexIdx.x, vertexIdx.y + 2, vertexIdx.z] +
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] +
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; Scalar inv_ds = AC_inv_dsz;
return (Scalar){Scalar(1.0 / 60.0) * inv_ds * return (Scalar){(1.0 / 60.0) * inv_ds *
(-Scalar(20.0) * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] + (-20.0 * vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z] +
Scalar(15.0) * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 1] + 15.0 * (vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 1] +
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] + 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 - 2]) +
vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 3] + vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z + 3] +
vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 3])}; vertex[vertexIdx.x, vertexIdx.y, vertexIdx.z - 3])};
} }
#endif #endif
#if LUPWD #if LUPWD
Device Scalar Device Scalar
upwd_der6(in VectorField uu, in ScalarField lnrho) upwd_der6(in VectorField uu, in ScalarField lnrho)
@@ -202,42 +201,44 @@ gradients(in VectorField uu)
#if LSINK #if LSINK
Device Vector Device Vector
sink_gravity(int3 globalVertexIdx){ sink_gravity(int3 globalVertexIdx)
{
int accretion_switch = int(AC_switch_accretion); int accretion_switch = int(AC_switch_accretion);
if (accretion_switch == 1){ if (accretion_switch == 1) {
Vector force_gravity; 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.y - AC_ny_min) * AC_dsy,
(globalVertexIdx.z - AC_nz_min) * AC_dsz}; (globalVertexIdx.z - AC_nz_min) * AC_dsz};
const Scalar sink_mass = AC_M_sink; const Scalar sink_mass = AC_M_sink;
const Vector sink_pos = (Vector){AC_sink_pos_x, const Vector sink_pos = (Vector){AC_sink_pos_x, AC_sink_pos_y, AC_sink_pos_z};
AC_sink_pos_y, const Scalar distance = length(grid_pos - sink_pos);
AC_sink_pos_z}; const Scalar soft = AC_soft;
const Scalar distance = length(grid_pos - sink_pos); // MV: The commit 083ff59 had AC_G_const defined wrong here in DSL making it exxessively
const Scalar soft = AC_soft; // strong. MV: Scalar gravity_magnitude = ... below is correct!
//MV: The commit 083ff59 had AC_G_const defined wrong here in DSL making it exxessively strong. const Scalar gravity_magnitude = (AC_G_const * sink_mass) /
//MV: Scalar gravity_magnitude = ... below is correct! pow(((distance * distance) + soft * soft), 1.5);
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, const Vector direction = (Vector){(sink_pos.x - grid_pos.x) / distance,
(sink_pos.y - grid_pos.y) / distance, (sink_pos.y - grid_pos.y) / distance,
(sink_pos.z - grid_pos.z) / distance}; (sink_pos.z - grid_pos.z) / distance};
force_gravity = gravity_magnitude * direction; force_gravity = gravity_magnitude * direction;
return force_gravity; return force_gravity;
} else { }
else {
return (Vector){0.0, 0.0, 0.0}; return (Vector){0.0, 0.0, 0.0};
} }
} }
#endif #endif
#if LSINK #if LSINK
// Give Truelove density // Give Truelove density
Device Scalar Device Scalar
truelove_density(in ScalarField lnrho){ truelove_density(in ScalarField lnrho)
const Scalar rho = exp(value(lnrho)); {
const Scalar rho = exp(value(lnrho));
const Scalar Jeans_length_squared = (M_PI * AC_cs2_sound) / (AC_G_const * rho); 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); const Scalar TJ_rho = ((M_PI) * ((AC_dsx * AC_dsx) / Jeans_length_squared) * AC_cs2_sound) /
//TODO: AC_dsx will cancel out, deal with it later for optimization. (AC_G_const * AC_dsx * AC_dsx);
// TODO: AC_dsx will cancel out, deal with it later for optimization.
Scalar accretion_rho = TJ_rho; Scalar accretion_rho = TJ_rho;
@@ -246,90 +247,92 @@ truelove_density(in ScalarField lnrho){
// This controls accretion of density/mass to the sink particle. // This controls accretion of density/mass to the sink particle.
Device Scalar Device Scalar
sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt)
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.y - AC_ny_min) * AC_dsy,
(globalVertexIdx.z - AC_nz_min) * AC_dsz}; (globalVertexIdx.z - AC_nz_min) * AC_dsz};
const Vector sink_pos = (Vector){AC_sink_pos_x, const Vector sink_pos = (Vector){AC_sink_pos_x, AC_sink_pos_y, AC_sink_pos_z};
AC_sink_pos_y, const Scalar profile_range = AC_accretion_range;
AC_sink_pos_z};
const Scalar profile_range = AC_accretion_range;
const Scalar accretion_distance = length(grid_pos - sink_pos); 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 accretion_density;
Scalar weight; Scalar weight;
if (accretion_switch == 1){ if (accretion_switch == 1) {
if ((accretion_distance) <= profile_range){ if ((accretion_distance) <= profile_range) {
//weight = Scalar(1.0); // weight = 1.0;
//Hann window function // Hann window function
Scalar window_ratio = accretion_distance/profile_range; Scalar window_ratio = accretion_distance / profile_range;
weight = Scalar(0.5)*(Scalar(1.0) - cos(Scalar(2.0)*M_PI*window_ratio)); weight = 0.5 * (1.0 - cos(2.0 * M_PI * window_ratio));
} else { }
weight = Scalar(0.0); 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)); const Scalar lnrho_min = log(truelove_density(lnrho));
Scalar rate; Scalar rate;
if (value(lnrho) > lnrho_min) { if (value(lnrho) > lnrho_min) {
rate = (exp(value(lnrho)) - exp(lnrho_min)) / dt; rate = (exp(value(lnrho)) - exp(lnrho_min)) / dt;
} else {
rate = Scalar(0.0);
} }
accretion_density = weight * rate ; else {
} else { rate = 0.0;
accretion_density = Scalar(0.0); }
accretion_density = weight * rate;
}
else {
accretion_density = 0.0;
} }
return accretion_density; return accretion_density;
} }
// This controls accretion of velocity to the sink particle. // This controls accretion of velocity to the sink particle.
Device Vector Device Vector
sink_accretion_velocity(int3 globalVertexIdx, in VectorField uu, Scalar dt) { sink_accretion_velocity(int3 globalVertexIdx, in VectorField uu, Scalar dt)
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.y - AC_ny_min) * AC_dsy,
(globalVertexIdx.z - AC_nz_min) * AC_dsz}; (globalVertexIdx.z - AC_nz_min) * AC_dsz};
const Vector sink_pos = (Vector){AC_sink_pos_x, const Vector sink_pos = (Vector){AC_sink_pos_x, AC_sink_pos_y, AC_sink_pos_z};
AC_sink_pos_y, const Scalar profile_range = AC_accretion_range;
AC_sink_pos_z};
const Scalar profile_range = AC_accretion_range;
const Scalar accretion_distance = length(grid_pos - sink_pos); const Scalar accretion_distance = length(grid_pos - sink_pos);
int accretion_switch = AC_switch_accretion; int accretion_switch = AC_switch_accretion;
Vector accretion_velocity; Vector accretion_velocity;
if (accretion_switch == 1){ if (accretion_switch == 1) {
Scalar weight; Scalar weight;
// Step function weighting // Step function weighting
// Arch of a cosine function? // Arch of a cosine function?
// Cubic spline x^3 - x in range [-0.5 , 0.5] // Cubic spline x^3 - x in range [-0.5 , 0.5]
if ((accretion_distance) <= profile_range){ if ((accretion_distance) <= profile_range) {
//weight = Scalar(1.0); // weight = 1.0;
//Hann window function // Hann window function
Scalar window_ratio = accretion_distance/profile_range; Scalar window_ratio = accretion_distance / profile_range;
weight = Scalar(0.5)*(Scalar(1.0) - cos(Scalar(2.0)*M_PI*window_ratio)); weight = 0.5 * (1.0 - cos(2.0 * M_PI * window_ratio));
} else { }
weight = Scalar(0.0); else {
weight = 0.0;
} }
Vector rate; Vector rate;
// MV: Could we use divergence here ephasize velocitie which are compressive and // MV: Could we use divergence here ephasize velocitie which are compressive and
// MV: not absorbins stuff that would not be accreted anyway? // MV: not absorbins stuff that would not be accreted anyway?
if (length(value(uu)) > Scalar(0.0)) { if (length(value(uu)) > 0.0) {
rate = (Scalar(1.0)/dt) * value(uu); rate = (1.0 / dt) * value(uu);
} else { }
else {
rate = (Vector){0.0, 0.0, 0.0}; rate = (Vector){0.0, 0.0, 0.0};
} }
accretion_velocity = weight * rate ; accretion_velocity = weight * rate;
} else { }
else {
accretion_velocity = (Vector){0.0, 0.0, 0.0}; accretion_velocity = (Vector){0.0, 0.0, 0.0};
} }
return accretion_velocity; return accretion_velocity;
} }
#endif #endif
Device Scalar Device Scalar
continuity(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar dt) 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) + upwd_der6(uu, lnrho)
#endif #endif
#if LSINK #if LSINK
- sink_accretion(globalVertexIdx, lnrho, dt) / exp(value(lnrho)) - sink_accretion(globalVertexIdx, lnrho, dt) / exp(value(lnrho))
#endif #endif
- divergence(uu); - divergence(uu);
} }
#if LENTROPY #if LENTROPY
Device Vector 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 Matrix S = stress_tensor(uu);
const Scalar cs2 = AC_cs2_sound * exp(AC_gamma * value(ss) / AC_cp_sound + const Scalar cs2 = AC_cs2_sound * exp(AC_gamma * value(ss) / AC_cp_sound +
(AC_gamma - 1) * (value(lnrho) - AC_lnrho0)); (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 (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
const Vector B = curl(aa); const Vector B = curl(aa);
// TODO: DOES INTHERMAL VERSTION INCLUDE THE MAGNETIC FIELD? // 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]*)\) // Regex replace CPU constants with get\(AC_([a-zA-Z_0-9]*)\)
// \1 // \1
const Vector mom = -mul(gradients(uu), value(uu)) - 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) + inv_rho * cross(j, B) +
AC_nu_visc * AC_nu_visc * (laplace_vec(uu) + (1.0 / 3.0) * gradient_of_divergence(uu) +
(laplace_vec(uu) + Scalar(1.0 / 3.0) * gradient_of_divergence(uu) + 2.0 * mul(S, gradient(lnrho))) +
Scalar(2.0) * mul(S, gradient(lnrho))) +
AC_zeta * gradient_of_divergence(uu) AC_zeta * gradient_of_divergence(uu)
#if LSINK #if LSINK
//Gravity term // Gravity term
+ sink_gravity(globalVertexIdx) + sink_gravity(globalVertexIdx)
//Corresponding loss of momentum // Corresponding loss of momentum
- //(Scalar(1.0) / Scalar( (AC_dsx*AC_dsy*AC_dsz) * exp(value(lnrho)))) * // Correction factor by unit mass - //(1.0 / ( (AC_dsx*AC_dsy*AC_dsz) * exp(value(lnrho)))) * //
sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014) // Correction factor by unit mass
; sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014)
#else ;
; #else
#endif ;
#endif
return mom; return mom;
} }
#elif LTEMPERATURE #elif LTEMPERATURE
@@ -392,14 +394,14 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in Scala
(gradient(tt) + value(tt) * gradient(lnrho)); (gradient(tt) + value(tt) * gradient(lnrho));
mom = -mul(gradients(uu), value(uu)) - pressure_term + mom = -mul(gradients(uu), value(uu)) - pressure_term +
AC_nu_visc * (laplace_vec(uu) + Scalar(1.0 / 3.0) * gradient_of_divergence(uu) + AC_nu_visc * (laplace_vec(uu) + (1.0 / 3.0) * gradient_of_divergence(uu) +
Scalar(2.0) * mul(S, gradient(lnrho))) + 2.0 * mul(S, gradient(lnrho))) +
AC_zeta * gradient_of_divergence(uu) AC_zeta * gradient_of_divergence(uu)
#if LSINK #if LSINK
+ sink_gravity(globalVertexIdx); + sink_gravity(globalVertexIdx);
#else #else
; ;
#endif #endif
#if LGRAVITY #if LGRAVITY
mom = mom - (Vector){0, 0, -10.0}; 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 // Isothermal: we have constant speed of sound
mom = -mul(gradients(uu), value(uu)) - AC_cs2_sound * gradient(lnrho) + 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) + AC_nu_visc * (laplace_vec(uu) + (1.0 / 3.0) * gradient_of_divergence(uu) +
Scalar(2.0) * mul(S, gradient(lnrho))) + 2.0 * mul(S, gradient(lnrho))) +
AC_zeta * gradient_of_divergence(uu) AC_zeta * gradient_of_divergence(uu)
#if LSINK #if LSINK
+ sink_gravity(globalVertexIdx) + sink_gravity(globalVertexIdx)
//Corresponding loss of momentum // Corresponding loss of momentum
- //(Scalar(1.0) / Scalar( (AC_dsx*AC_dsy*AC_dsz) * exp(value(lnrho)))) * // Correction factor by unit mass - //(1.0 / ( (AC_dsx*AC_dsy*AC_dsz) * exp(value(lnrho)))) * // Correction
sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014) // factor by unit mass
; sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014)
#else ;
; #else
#endif ;
#endif
#if LGRAVITY #if LGRAVITY
mom = mom - (Vector){0, 0, -10.0}; mom = mom - (Vector){0, 0, -10.0};
@@ -460,21 +463,21 @@ Device Scalar
lnT(in ScalarField ss, in ScalarField lnrho) lnT(in ScalarField ss, in ScalarField lnrho)
{ {
return AC_lnT0 + AC_gamma * value(ss) / AC_cp_sound + 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) // Nabla dot (K nabla T) / (rho T)
Device Scalar Device Scalar
heat_conduction(in ScalarField ss, in ScalarField lnrho) 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 Vector grad_ln_chi = -gradient(lnrho);
const Scalar first_term = AC_gamma * inv_AC_cp_sound * laplace(ss) + 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) + 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)) + const Vector third_term = AC_gamma * (inv_AC_cp_sound * gradient(ss) + gradient(lnrho)) +
grad_ln_chi; grad_ln_chi;
@@ -492,11 +495,11 @@ Device Scalar
entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa) entropy(in ScalarField ss, in VectorField uu, in ScalarField lnrho, in VectorField aa)
{ {
const Matrix S = stress_tensor(uu); const Matrix S = stress_tensor(uu);
const Scalar inv_pT = Scalar(1.0) / (exp(value(lnrho)) * exp(lnT(ss, lnrho))); const Scalar inv_pT = 1.0 / (exp(value(lnrho)) * exp(lnT(ss, lnrho)));
const Vector j = (Scalar(1.0) / AC_mu0) * const Vector j = (1.0 / AC_mu0) *
(gradient_of_divergence(aa) - laplace_vec(aa)); // Current density (gradient_of_divergence(aa) - laplace_vec(aa)); // Current density
const Scalar RHS = H_CONST - C_CONST + AC_eta * (AC_mu0)*dot(j, j) + 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); AC_zeta * exp(value(lnrho)) * divergence(uu) * divergence(uu);
return -dot(value(uu), gradient(ss)) + inv_pT * RHS + heat_conduction(ss, lnrho); 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; const Scalar heat_diffusivity_k = 0.0008; // 8e-4;
return -dot(value(uu), gradient(tt)) + heat_diffusivity_k * laplace(tt) + return -dot(value(uu), gradient(tt)) + heat_diffusivity_k * laplace(tt) +
heat_diffusivity_k * dot(gradient(lnrho), gradient(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); (AC_gamma - 1) * value(tt) * divergence(uu);
} }
#endif #endif
#if LFORCING #if LFORCING
Device Vector 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; int accretion_switch = AC_switch_accretion;
if (accretion_switch == 0){ if (accretion_switch == 0) {
return magnitude * cross(normalized(b - a), (Vector){ 0, 0, 1}); // Vortex return magnitude * cross(normalized(b - a), (Vector){0, 0, 1}); // Vortex
} else { }
return (Vector){0,0,0}; else {
return (Vector){0, 0, 0};
} }
} }
Device Vector 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; 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 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) forcing(int3 globalVertexIdx, Scalar dt)
{ {
int accretion_switch = AC_switch_accretion; int accretion_switch = AC_switch_accretion;
if (accretion_switch == 0){ if (accretion_switch == 0) {
Vector a = Scalar(0.5) * (Vector){globalGridN.x * AC_dsx, Vector a = 0.5 * (Vector){globalGridN.x * AC_dsx, globalGridN.y * AC_dsy,
globalGridN.y * AC_dsy, globalGridN.z * AC_dsz}; // source (origin)
globalGridN.z * AC_dsz}; // source (origin) Vector xx = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx,
Vector xx = (Vector){(globalVertexIdx.x - AC_nx_min) * AC_dsx,
(globalVertexIdx.y - AC_ny_min) * AC_dsy, (globalVertexIdx.y - AC_ny_min) * AC_dsy,
(globalVertexIdx.z - AC_nz_min) * AC_dsz}; // sink (current index) (globalVertexIdx.z - AC_nz_min) * AC_dsz}; // sink (current index)
const Scalar cs2 = AC_cs2_sound; 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 magnitude = AC_forcing_magnitude;
Scalar phase = AC_forcing_phase; 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_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}; 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. // Scaling N = magnitude*cs*sqrt(k*cs/dt) * dt
//Vector force = simple_vortex_forcing(a, xx, magnitude); const Scalar NN = cs * sqrt(AC_kaver * cs);
//Vector force = simple_outward_flow_forcing(a, xx, magnitude); // MV: Like in the Pencil Code. I don't understandf the logic here.
Vector force = helical_forcing(magnitude, k_force, xx, ff_re,ff_im, phase); 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 if (is_valid(force)) {
const Scalar NN = cs*sqrt(AC_kaver*cs); return force;
//MV: Like in the Pencil Code. I don't understandf the logic here. }
force.x = sqrt(dt)*NN*force.x; else {
force.y = sqrt(dt)*NN*force.y; return (Vector){0, 0, 0};
force.z = sqrt(dt)*NN*force.z; }
}
if (is_valid(force)) { return force; } else {
else { return (Vector){0, 0, 0}; } return (Vector){0, 0, 0};
} else {
return (Vector){0,0,0};
} }
} }
#endif // LFORCING #endif // LFORCING
@@ -671,10 +681,11 @@ solve()
#endif #endif
#if LSINK #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) { 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 #endif
} }

3
acc/src/acc.l
View File

@@ -39,7 +39,8 @@ L [a-zA-Z_]
"return" { return RETURN; } "return" { return RETURN; }
{D}+"."{D}+ { return REAL_NUMBER; } /* Literals */ {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; } {D}+[lu]* { return NUMBER; }
{L}({L}|{D})* { return IDENTIFIER; } {L}({L}|{D})* { return IDENTIFIER; }
\"(.)*\" { return IDENTIFIER; } /* String */ \"(.)*\" { return IDENTIFIER; } /* String */

3
acc/src/acc.y
View File

@@ -14,7 +14,7 @@ int yyget_lineno();
%} %}
%token CONSTANT IN OUT UNIFORM %token CONSTANT IN OUT UNIFORM
%token IDENTIFIER NUMBER REAL_NUMBER %token IDENTIFIER NUMBER REAL_NUMBER DOUBLE_NUMBER
%token RETURN %token RETURN
%token SCALAR VECTOR MATRIX SCALARFIELD SCALARARRAY %token SCALAR VECTOR MATRIX SCALARFIELD SCALARARRAY
%token VOID INT INT3 COMPLEX %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, $$); } 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, $$); } | NUMBER { $$ = astnode_create(NODE_UNKNOWN, NULL, NULL); astnode_set_buffer(yytext, $$); }
; ;

18
acc/stdlib/stdderiv.h
View File

@@ -302,18 +302,18 @@ stress_tensor(in VectorField vec)
{ {
Matrix S; Matrix S;
S.row[0].x = Scalar(2.0 / 3.0) * gradient(vec.x).x - S.row[0].x = (2.0 / 3.0) * gradient(vec.x).x -
Scalar(1.0 / 3.0) * (gradient(vec.y).y + gradient(vec.z).z); (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].y = (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].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 - S.row[1].y = (2.0 / 3.0) * gradient(vec.y).y -
Scalar(1.0 / 3.0) * (gradient(vec.x).x + gradient(vec.z).z); (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 - S.row[2].z = (2.0 / 3.0) * gradient(vec.z).z -
Scalar(1.0 / 3.0) * (gradient(vec.x).x + gradient(vec.y).y); (1.0 / 3.0) * (gradient(vec.x).x + gradient(vec.y).y);
S.row[1].x = S.row[0].y; S.row[1].x = S.row[0].y;
S.row[2].x = S.row[0].z; S.row[2].x = S.row[0].z;