Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add geometric source terms to the PPM tracing #2473

Merged
merged 18 commits into from
Sep 26, 2024
Merged

add geometric source terms to the PPM tracing #2473

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
4eabeab
add geometric source terms to the PPM tracing
zingale May 30, 2023
2a44c10
finish implementation
zingale May 30, 2023
f0647a6
fix 3d compilation
zingale May 30, 2023
8516265
Merge branch 'development' into ppm_source_tracing
zingale Jun 7, 2023
d41376e
Merge branch 'development' into ppm_source_tracing
zingale Jun 17, 2023
dd08912
Merge branch 'development' into ppm_source_tracing
zingale Jun 17, 2023
12d552d
Merge branch 'development' into ppm_source_tracing
zingale Jun 18, 2023
7b0343c
Merge branch 'development' into ppm_source_tracing
zingale Jul 10, 2023
01f092c
Merge branch 'development' into ppm_source_tracing
zingale Nov 8, 2023
4c0b75d
Merge branch 'development' into ppm_source_tracing
zingale Jul 26, 2024
02e3a9e
Merge branch 'development' into ppm_source_tracing
zingale Jul 27, 2024
8828c4b
fix pressure source
zingale Jul 27, 2024
ec64d4f
Merge branch 'development' into ppm_source_tracing
zingale Aug 21, 2024
2c57bc8
remove unused var
zingale Aug 21, 2024
83e735b
Merge branch 'ppm_source_tracing' of github.com:zingale/Castro into p…
zingale Aug 21, 2024
8e914c3
update CI
zingale Aug 21, 2024
3a4662b
Merge branch 'development' into ppm_source_tracing
zingale Sep 16, 2024
267d57c
Merge branch 'development' into ppm_source_tracing
zhichen3 Sep 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions Exec/science/flame_wave/ci-benchmarks/grid_diag.out
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are we not multiplying the source by dt/2 ?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is just adding the geometric sources to the other hydrodynamic sources. Then we reconstruct and integrate under the parabola, and then the dt2 is applied when we add it to the jumps that are accumulated to the interfaces

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
Open in desktop
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

});
}