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

Autoformatting

Browse Source
This commit is contained in:
jpekkila
2019-10-07 17:17:58 +03:00
parent c8e0586b60
commit 16c8b1e748
2 changed files with 133 additions and 119 deletions
Download Patch File Download Diff File Expand all files Collapse all files

104
acc/mhd_solver/stencil_process.ac
View File

@@ -1,7 +1,7 @@
#include <stdderiv.h> #include <stdderiv.h>
#include "stencil_definition.h"
#include "stencil_assembly.h" #include "stencil_assembly.h"
#include "stencil_definition.h"
#if LUPWD #if LUPWD
Device Scalar Device Scalar
@@ -22,7 +22,8 @@ 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;
@@ -30,33 +31,34 @@ sink_gravity(int3 globalVertexIdx){
(globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy, (globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy,
(globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz}; (globalVertexIdx.z - DCONST(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,
AC_sink_pos_z};
const Scalar distance = length(grid_pos - sink_pos); const Scalar distance = length(grid_pos - sink_pos);
const Scalar soft = AC_soft; const Scalar soft = AC_soft;
//MV: The commit 083ff59 had AC_G_const defined wrong here in DSL making it exxessively strong. // MV: The commit 083ff59 had AC_G_const defined wrong here in DSL making it exxessively
//MV: Scalar gravity_magnitude = ... below is correct! // 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 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) /
(AC_G_const * AC_dsx * AC_dsx);
// TODO: AC_dsx will cancel out, deal with it later for optimization. // TODO: AC_dsx will cancel out, deal with it later for optimization.
Scalar accretion_rho = TJ_rho; Scalar accretion_rho = TJ_rho;
@@ -66,13 +68,12 @@ 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 - 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.y - DCONST(AC_ny_min)) * AC_dsy,
(globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz}; (globalVertexIdx.z - DCONST(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,
AC_sink_pos_z};
const Scalar profile_range = AC_accretion_range; 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;
@@ -85,7 +86,8 @@ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
// 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 = Scalar(0.5) * (Scalar(1.0) - cos(Scalar(2.0) * M_PI * window_ratio));
} else { }
else {
weight = Scalar(0.0); weight = Scalar(0.0);
} }
@@ -94,11 +96,13 @@ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
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 { }
else {
rate = Scalar(0.0); rate = Scalar(0.0);
} }
accretion_density = weight * rate; accretion_density = weight * rate;
} else { }
else {
accretion_density = Scalar(0.0); accretion_density = Scalar(0.0);
} }
return accretion_density; return accretion_density;
@@ -106,13 +110,12 @@ sink_accretion(int3 globalVertexIdx, in ScalarField lnrho, Scalar dt){
// 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 - 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.y - DCONST(AC_ny_min)) * AC_dsy,
(globalVertexIdx.z - DCONST(AC_nz_min)) * AC_dsz}; (globalVertexIdx.z - DCONST(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,
AC_sink_pos_z};
const Scalar profile_range = AC_accretion_range; 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;
@@ -128,28 +131,29 @@ sink_accretion_velocity(int3 globalVertexIdx, in VectorField uu, Scalar dt) {
// 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 = Scalar(0.5) * (Scalar(1.0) - cos(Scalar(2.0) * M_PI * window_ratio));
} else { }
else {
weight = Scalar(0.0); weight = Scalar(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)) > Scalar(0.0)) {
rate = (Scalar(1.0) / dt) * value(uu); rate = (Scalar(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)
{ {
@@ -164,11 +168,10 @@ continuity(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar
- 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 +
@@ -192,7 +195,8 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, in Scala
// 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 - //(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) sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014)
; ;
#else #else
@@ -243,7 +247,8 @@ momentum(int3 globalVertexIdx, in VectorField uu, in ScalarField lnrho, Scalar d
#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 - //(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) sink_accretion_velocity(globalVertexIdx, uu, dt) // As in Lee et al.(2014)
; ;
#else #else
@@ -339,21 +344,25 @@ heat_transfer(in VectorField uu, in ScalarField lnrho, in ScalarField tt)
#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 { }
else {
return (Vector){0, 0, 0}; 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 { }
else {
return (Vector){0, 0, 0}; return (Vector){0, 0, 0};
} }
} }
@@ -400,12 +409,12 @@ 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 = Scalar(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 - DCONST(AC_nx_min)) * AC_dsx, Vector xx = (Vector){(globalVertexIdx.x - DCONST(AC_nx_min)) * AC_dsx,
(globalVertexIdx.y - DCONST(AC_ny_min)) * AC_dsy, (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 cs2 = AC_cs2_sound;
const Scalar cs = sqrt(cs2); const Scalar cs = sqrt(cs2);
@@ -416,7 +425,6 @@ forcing(int3 globalVertexIdx, Scalar dt)
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. // Determine that forcing funtion type at this point.
// Vector force = simple_vortex_forcing(a, xx, magnitude); // Vector force = simple_vortex_forcing(a, xx, magnitude);
// Vector force = simple_outward_flow_forcing(a, xx, magnitude); // Vector force = simple_outward_flow_forcing(a, xx, magnitude);
@@ -429,9 +437,14 @@ forcing(int3 globalVertexIdx, Scalar dt)
force.y = sqrt(dt) * NN * force.y; force.y = sqrt(dt) * NN * force.y;
force.z = sqrt(dt) * NN * force.z; force.z = sqrt(dt) * NN * force.z;
if (is_valid(force)) { return force; } if (is_valid(force)) {
else { return (Vector){0, 0, 0}; } return force;
} else { }
else {
return (Vector){0, 0, 0};
}
}
else {
return (Vector){0, 0, 0}; return (Vector){0, 0, 0};
} }
} }
@@ -492,7 +505,8 @@ 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!