Skip to content

Commit

Permalink
add geometric source terms to the PPM tracing (#2473)
Browse files Browse the repository at this point in the history 
Add the geometric source terms to the primitive variable sources so they can be
included in the tracing.

We can't simply add them in src_to_prim, since the sources depend on the
direction of the reconstruction. We also can't just add them in the trace_ppm
to srcQ immediately at the start and then do the reconstruction since we need
the geometric source terms in the ghost cells for the reconstruction. Therefore
we add these geometric sources to the 5-point stencil that will be traced just
before the tracing. This requires some redundant calculation, but it is not too
bad.
  • Loading branch information
@zingale
zingale authored Sep 26, 2024
1 parent fb34ee1 commit b9a0114
Show file tree
Hide file tree
Showing 5 changed files with 137 additions and 73 deletions.
4 changes: 2 additions & 2 deletions Exec/science/flame_wave/ci-benchmarks/grid_diag.out
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# COLUMN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
# TIMESTEP TIME MASS XMOM YMOM ZMOM ANG. MOM. X ANG. MOM. Y ANG. MOM. Z KIN. ENERGY INT. ENERGY GAS ENERGY GRAV. ENERGY TOTAL ENERGY CENTER OF MASS X-LOC CENTER OF MASS Y-LOC CENTER OF MASS Z-LOC CENTER OF MASS X-VEL CENTER OF MASS Y-VEL CENTER OF MASS Z-VEL MAXIMUM TEMPERATURE MAXIMUM DENSITY MAXIMUM T_S / T_E
0 0.0000000000000000e+00 1.4700685690736437e+20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.9163021935541997e+37 2.9163021935541997e+37 0.0000000000000000e+00 2.9163021935541997e+37 2.7306641645125415e+04 6.8087525965685256e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.3928678114307070e+09 3.0715027784567930e+07 0.0000000000000000e+00
1 5.9806047036494849e-08 1.4700700110495903e+20 3.3564221294760178e+23 7.2876252590129013e+23 1.2612378576902532e+20 3.1318392026636616e+23 -2.7141686393328650e+24 1.9066097490832445e+28 7.8605333893026802e+29 2.9163081863790971e+37 2.9163082649844082e+37 0.0000000000000000e+00 2.9163082649844082e+37 2.7306641717925890e+04 6.8087482283157362e+02 0.0000000000000000e+00 2.2831716205676648e+03 4.9573321027137572e+03 8.5794407627549896e-01 1.3929561430929687e+09 3.0715326643214259e+07 3.3873763068343888e-04
2 1.2559269877663918e-07 1.4700717750115456e+20 7.0484895627215180e+23 1.5304335700608418e+24 5.5619062084401319e+20 1.3810901029080806e+24 -1.1969287510263666e+25 4.0039912973182075e+28 3.0853296206485770e+30 2.9163145904134766e+37 2.9163148989464290e+37 0.0000000000000000e+00 2.9163148989464290e+37 2.7306641954244220e+04 6.8087461718998884e+02 0.0000000000000000e+00 4.7946567525018709e+03 1.0410604407725754e+04 3.7834249340624542e+00 1.3933944931271181e+09 3.0715562157661017e+07 3.3857586853705914e-04
1 5.9806047036494849e-08 1.4700700110495903e+20 3.3564221217849212e+23 7.2876252590226133e+23 1.2612378556916341e+20 3.1318391966955516e+23 -2.7141686347841742e+24 1.9066097493127676e+28 7.8605324800324878e+29 2.9163081863790706e+37 2.9163082649843738e+37 0.0000000000000000e+00 2.9163082649843738e+37 2.7306641717925955e+04 6.8087482283157317e+02 0.0000000000000000e+00 2.2831716153358752e+03 4.9573321027203638e+03 8.5794407491595881e-01 1.3929561431310942e+09 3.0715326643212341e+07 3.3873763041849706e-04
2 1.2559269877663918e-07 1.4700717750115462e+20 7.0484895445039834e+23 1.5304335700823335e+24 5.5619061962888426e+20 1.3810900989697918e+24 -1.1969287484132617e+25 4.0039912979259806e+28 3.0853275723280981e+30 2.9163145904135894e+37 2.9163148989463411e+37 0.0000000000000000e+00 2.9163148989463411e+37 2.7306641954244522e+04 6.8087461718998338e+02 0.0000000000000000e+00 4.7946567401095926e+03 1.0410604407871944e+04 3.7834249257966728e+00 1.3933944939675827e+09 3.0715562157563787e+07 3.3857586854903492e-04
4 changes: 2 additions & 2 deletions Exec/science/flame_wave/ci-benchmarks/species_diag.out
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# COLUMN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
# TIMESTEP TIME Mass He4 Mass C12 Mass O16 Mass Ne20 Mass Mg24 Mass Si28 Mass S32 Mass Ar36 Mass Ca40 Mass Ti44 Mass Cr48 Mass Fe52 Mass Ni56
0 0.0000000000000000e+00 2.8142378112400895e-15 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.1117997526547418e-14
1 5.9806047036494849e-08 2.8142323856670296e-15 5.4329394513687282e-21 7.4194476167127106e-24 7.3934515289195991e-24 7.3935541950036108e-24 7.3933935990475538e-24 7.3932321468746301e-24 7.3932312019180565e-24 7.3932308328287436e-24 7.3932307869393029e-24 7.3932307850741563e-24 7.3932307849933102e-24 7.1118070045957091e-14
2 1.2559269877663918e-07 2.8142264165380905e-15 1.1401993461614438e-20 7.4929577845086585e-24 7.3944696393124482e-24 7.3939305113509452e-24 7.3935838899709814e-24 7.3932425308155010e-24 7.3932405326302037e-24 7.3932397568503172e-24 7.3932396603620439e-24 7.3932396564405536e-24 7.3932396562705448e-24 7.1118158758588142e-14
1 5.9806047036494849e-08 2.8142323856670446e-15 5.4329394398418374e-21 7.4194443455948811e-24 7.3934514829832856e-24 7.3935541949306949e-24 7.3933935989353088e-24 7.3932321468740865e-24 7.3932312019180345e-24 7.3932308328287333e-24 7.3932307869392985e-24 7.3932307850741519e-24 7.3932307849933044e-24 7.1118070045957103e-14
2 1.2559269877663918e-07 2.8142264165381050e-15 1.1401993450709496e-20 7.4929528148974785e-24 7.3944694290388374e-24 7.3939305042586412e-24 7.3935838894001731e-24 7.3932425308122493e-24 7.3932405326301655e-24 7.3932397568503098e-24 7.3932396603620498e-24 7.3932396564405566e-24 7.3932396562705492e-24 7.1118158758588142e-14
46 changes: 23 additions & 23 deletions Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out
View file
Original file line number Diff line number Diff line change
@@ -1,29 +1,29 @@
plotfile = plt00086
time = 1.25
variables minimum value maximum value
density 8.693611703e-05 19565534.299
xmom -5.4964100651e+14 1.3559128302e+14
ymom -2.5530096328e+15 2.5530122744e+15
density 8.6936468335e-05 19568007.309
xmom -5.4959005475e+14 1.3559838981e+14
ymom -2.5540206053e+15 2.5540214496e+15
zmom 0 0
rho_E 7.4982062146e+11 5.0669247219e+24
rho_e 7.1077581849e+11 5.0640768326e+24
Temp 242288.68588 1409652233.6
rho_He4 8.693611703e-17 3.599903302
rho_C12 3.4774446812e-05 7825956.6934
rho_O16 5.2161670217e-05 11739149.75
rho_Ne20 8.693611703e-17 181951.0571
rho_Mg24 8.693611703e-17 1192.7969729
rho_Si28 8.693611703e-17 6.6913702949
rho_S32 8.693611703e-17 0.00019493291655
rho_Ar36 8.693611703e-17 1.9565534609e-05
rho_Ca40 8.693611703e-17 1.9565534331e-05
rho_Ti44 8.693611703e-17 1.9565534308e-05
rho_Cr48 8.693611703e-17 1.9565534308e-05
rho_Fe52 8.693611703e-17 1.9565534308e-05
rho_Ni56 8.693611703e-17 1.9565534308e-05
phiGrav -5.8708033902e+17 -2.3375498341e+16
grav_x -684991644 -51428.243166
grav_y -739606241.84 739606820.44
rho_E 7.4987846833e+11 5.0676543192e+24
rho_e 7.1083104678e+11 5.0648015049e+24
Temp 242292.44287 1409581841.1
rho_He4 8.6936468335e-17 3.5973411848
rho_C12 3.4774587333e-05 7826946.398
rho_O16 5.2161881e-05 11740633.889
rho_Ne20 8.6936468335e-17 181819.49413
rho_Mg24 8.6936468335e-17 1190.7712386
rho_Si28 8.6936468335e-17 6.6817951193
rho_S32 8.6936468335e-17 0.00019451843176
rho_Ar36 8.6936468335e-17 1.9568007618e-05
rho_Ca40 8.6936468335e-17 1.9568007341e-05
rho_Ti44 8.6936468335e-17 1.9568007318e-05
rho_Cr48 8.6936468335e-17 1.9568007318e-05
rho_Fe52 8.6936468335e-17 1.9568007318e-05
rho_Ni56 8.6936468335e-17 1.9568007318e-05
phiGrav -5.8709462562e+17 -2.3375498549e+16
grav_x -685026429.13 -51428.265677
grav_y -739654246.49 739654206.24
grav_z 0 0
rho_enuc -2.7633982574e+12 7.6429034885e+23
rho_enuc -4.7815621457e+12 7.6360058391e+23

78 changes: 73 additions & 5 deletions Source/hydro/reconstruction.H
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
#ifndef CASTRO_RECONSTRUCTION_H
#define CASTRO_RECONSTRUCTION_H

#include <state_indices.H>

namespace reconstruction {
enum slope_indices {
im2 = 0,
Expand All @@ -14,9 +16,9 @@ namespace reconstruction {

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void
load_stencil(Array4<Real const> const& q_arr, const int idir,
load_stencil(amrex::Array4<amrex::Real const> const& q_arr, const int idir,
const int i, const int j, const int k, const int ncomp,
Real* s) {
amrex::Real* s) {

using namespace reconstruction;

Expand Down Expand Up @@ -47,9 +49,11 @@ load_stencil(Array4<Real const> const& q_arr, const int idir,

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void
load_passive_stencil(Array4<Real const> const& U_arr, Array4<Real const> const& rho_inv_arr, const int idir,
load_passive_stencil(amrex::Array4<amrex::Real const> const& U_arr,
amrex::Array4<amrex::Real const> const& rho_inv_arr,
const int idir,
const int i, const int j, const int k, const int ncomp,
Real* s) {
amrex::Real* s) {

using namespace reconstruction;

Expand Down Expand Up @@ -78,5 +82,69 @@ load_passive_stencil(Array4<Real const> const& U_arr, Array4<Real const> const&

}

#endif
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void
add_geometric_rho_source(amrex::Array4<amrex::Real const> const& q_arr,
amrex::Array4<amrex::Real const> const& dloga,
const int i, const int j, const int k,
amrex::Real* s) {

using namespace reconstruction;

// this takes the form: -alpha rho u / r
// where alpha = 1 for cylindrical and 2 for spherical

// note: this is assumed to be working only in the x-direction

s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QRHO) * q_arr(i-2,j,k,QU);
s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QRHO) * q_arr(i-1,j,k,QU);
s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QRHO) * q_arr(i,j,k,QU);
s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QRHO) * q_arr(i+1,j,k,QU);
s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QRHO) * q_arr(i+2,j,k,QU);
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void
add_geometric_rhoe_source(amrex::Array4<amrex::Real const> const& q_arr,
amrex::Array4<amrex::Real const> const& dloga,
const int i, const int j, const int k,
amrex::Real* s) {

using namespace reconstruction;

// this takes the form: -alpha (rho e + p) u / r
// where alpha = 1 for cylindrical and 2 for spherical

// note: this is assumed to be working only in the x-direction

s[im2] += -dloga(i-2,j,k) * (q_arr(i-2,j,k,QREINT) + q_arr(i-2,j,k,QPRES)) * q_arr(i-2,j,k,QU);
s[im1] += -dloga(i-1,j,k) * (q_arr(i-1,j,k,QREINT) + q_arr(i-1,j,k,QPRES)) * q_arr(i-1,j,k,QU);
s[i0] += -dloga(i,j,k) * (q_arr(i,j,k,QREINT) + q_arr(i,j,k,QPRES)) * q_arr(i,j,k,QU);
s[ip1] += -dloga(i+1,j,k) * (q_arr(i+1,j,k,QREINT) + q_arr(i+1,j,k,QPRES)) * q_arr(i+1,j,k,QU);
s[ip2] += -dloga(i+2,j,k) * (q_arr(i+2,j,k,QREINT) + q_arr(i+2,j,k,QPRES)) * q_arr(i+2,j,k,QU);
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void
add_geometric_p_source(amrex::Array4<amrex::Real const> const& q_arr,
amrex::Array4<amrex::Real const> const& qaux_arr,
amrex::Array4<amrex::Real const> const& dloga,
const int i, const int j, const int k,
amrex::Real* s) {

using namespace reconstruction;

// this takes the form: -alpha Gamma1 p u / r
// where alpha = 1 for cylindrical and 2 for spherical

// note: this is assumed to be working only in the x-direction

s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QPRES) * qaux_arr(i-2,j,k,QGAMC) * q_arr(i-2,j,k,QU);
s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QPRES) * qaux_arr(i-1,j,k,QGAMC) * q_arr(i-1,j,k,QU);
s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QPRES) * qaux_arr(i,j,k,QGAMC) * q_arr(i,j,k,QU);
s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QPRES) * qaux_arr(i+1,j,k,QGAMC) * q_arr(i+1,j,k,QU);
s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QPRES) * qaux_arr(i+2,j,k,QGAMC) * q_arr(i+2,j,k,QU);
}


#endif
78 changes: 37 additions & 41 deletions Source/hydro/trace_ppm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@ Castro::trace_ppm(const Box& bx,


const auto dx = geom.CellSizeArray();
const int coord = geom.Coord();

Real hdt = 0.5_rt * dt;
Real dtdx = dt / dx[idir];
Expand Down Expand Up @@ -82,6 +83,12 @@ Castro::trace_ppm(const Box& bx,
for (int n = 0; n < NQSRC; n++) {
do_source_trace[n] = 0;

// geometric source terms in r-direction need tracing
if (coord > 0 && idir == 0 && (n == QRHO || n == QPRES || n == QREINT)) {
do_source_trace[n] = 1;
continue;
}

for (int k = lo[2]-2*dg2; k <= hi[2]+2*dg2; k++) {
for (int j = lo[1]-2*dg1; j <= hi[1]+2*dg1; j++) {
for (int i = lo[0]-2; i <= hi[0]+2; i++) {
Expand Down Expand Up @@ -148,16 +155,9 @@ Castro::trace_ppm(const Box& bx,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{


Real cc = qaux_arr(i,j,k,QC);

#if AMREX_SPACEDIM < 3
Real csq = cc*cc;
#endif

Real un = q_arr(i,j,k,QUN);


// do the parabolic reconstruction and compute the
// integrals under the characteristic waves
Real s[nslp];
Expand Down Expand Up @@ -279,11 +279,20 @@ Castro::trace_ppm(const Box& bx,
#ifndef AMREX_USE_GPU
do_trace = do_source_trace[QRHO];
#else
do_trace = check_trace_source(srcQ, idir, i, j, k, QRHO);
if (idir == 0 && coord > 0) {
do_trace = 1;
} else {
do_trace = check_trace_source(srcQ, idir, i, j, k, QRHO);
}
#endif

if (do_trace) {
load_stencil(srcQ, idir, i, j, k, QRHO, s);
#if AMREX_SPACEDIM <= 2
if (idir == 0 && coord > 0) {
add_geometric_rho_source(q_arr, dloga, i, j, k, s);
}
#endif
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rho, Im_src_rho);
}
Expand Down Expand Up @@ -316,11 +325,20 @@ Castro::trace_ppm(const Box& bx,
#ifndef AMREX_USE_GPU
do_trace = do_source_trace[QPRES];
#else
do_trace = check_trace_source(srcQ, idir, i, j, k, QPRES);
if (idir == 0 && coord > 0) {
do_trace = 1;
} else {
do_trace = check_trace_source(srcQ, idir, i, j, k, QPRES);
}
#endif

if (do_trace) {
load_stencil(srcQ, idir, i, j, k, QPRES, s);
#if AMREX_SPACEDIM <= 2
if (idir == 0 && coord > 0) {
add_geometric_p_source(q_arr, qaux_arr, dloga, i, j, k, s);
}
#endif
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_p, Im_src_p);
}
Expand All @@ -333,11 +351,20 @@ Castro::trace_ppm(const Box& bx,
#ifndef AMREX_USE_GPU
do_trace = do_source_trace[QREINT];
#else
do_trace = check_trace_source(srcQ, idir, i, j, k, QREINT);
if (idir == 0 && coord > 0) {
do_trace = 1;
} else {
do_trace = check_trace_source(srcQ, idir, i, j, k, QREINT);
}
#endif

if (do_trace) {
load_stencil(srcQ, idir, i, j, k, QREINT, s);
#if AMREX_SPACEDIM <= 2
if (idir == 0 && coord > 0) {
add_geometric_rhoe_source(q_arr, dloga, i, j, k, s);
}
#endif
ppm_reconstruct(s, flat, sm, sp);
ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rhoe, Im_src_rhoe);
}
Expand Down Expand Up @@ -604,36 +631,5 @@ Castro::trace_ppm(const Box& bx,

}

// geometry source terms
#if (AMREX_SPACEDIM < 3)
// these only apply for x states (idir = 0)
if (idir == 0 && dloga(i,j,k) != 0.0_rt) {
Real rho = q_arr(i,j,k,QRHO);

Real courn = dt/dx[0]*(cc+std::abs(un));
Real eta = (1.0_rt - courn)/(cc*dt*std::abs(dloga(i,j,k)));
Real dlogatmp = amrex::min(eta, 1.0_rt)*dloga(i,j,k);
Real sourcr = -0.5_rt*dt*rho*dlogatmp*un;
Real sourcp = sourcr*csq;
Real source = sourcp*((q_arr(i,j,k,QPRES) + q_arr(i,j,k,QREINT))/rho)/csq;

if (i <= vhi[0]) {
qm(i+1,j,k,QRHO) = qm(i+1,j,k,QRHO) + sourcr;
qm(i+1,j,k,QRHO) = amrex::max(qm(i+1,j,k,QRHO), lsmall_dens);
qm(i+1,j,k,QPRES) = qm(i+1,j,k,QPRES) + sourcp;
qm(i+1,j,k,QREINT) = qm(i+1,j,k,QREINT) + source;
}

if (i >= vlo[0]) {
qp(i,j,k,QRHO) = qp(i,j,k,QRHO) + sourcr;
qp(i,j,k,QRHO) = amrex::max(qp(i,j,k,QRHO), lsmall_dens);
qp(i,j,k,QPRES) = qp(i,j,k,QPRES) + sourcp;
qp(i,j,k,QREINT) = qp(i,j,k,QREINT) + source;
}
}
#endif

});
}


0 comments on commit b9a0114

Please sign in to comment.