From b9a0114ab019c89e6be94983227fa137d14f8ac6 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 26 Sep 2024 14:20:16 -0400 Subject: [PATCH] add geometric source terms to the PPM tracing (#2473) 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. --- .../flame_wave/ci-benchmarks/grid_diag.out | 4 +- .../flame_wave/ci-benchmarks/species_diag.out | 4 +- .../ci-benchmarks/wdmerger_collision_2D.out | 46 +++++------ Source/hydro/reconstruction.H | 78 +++++++++++++++++-- Source/hydro/trace_ppm.cpp | 78 +++++++++---------- 5 files changed, 137 insertions(+), 73 deletions(-) diff --git a/Exec/science/flame_wave/ci-benchmarks/grid_diag.out b/Exec/science/flame_wave/ci-benchmarks/grid_diag.out index 973d3f2aba..3bfc78ca85 100644 --- a/Exec/science/flame_wave/ci-benchmarks/grid_diag.out +++ b/Exec/science/flame_wave/ci-benchmarks/grid_diag.out @@ -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 diff --git a/Exec/science/flame_wave/ci-benchmarks/species_diag.out b/Exec/science/flame_wave/ci-benchmarks/species_diag.out index 6baaafc600..f3c2dcfe7a 100644 --- a/Exec/science/flame_wave/ci-benchmarks/species_diag.out +++ b/Exec/science/flame_wave/ci-benchmarks/species_diag.out @@ -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 diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index 88fbaf4f88..50abd61af7 100644 --- a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out +++ b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out @@ -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 diff --git a/Source/hydro/reconstruction.H b/Source/hydro/reconstruction.H index 6ee405a39f..9a41710663 100644 --- a/Source/hydro/reconstruction.H +++ b/Source/hydro/reconstruction.H @@ -1,6 +1,8 @@ #ifndef CASTRO_RECONSTRUCTION_H #define CASTRO_RECONSTRUCTION_H +#include + namespace reconstruction { enum slope_indices { im2 = 0, @@ -14,9 +16,9 @@ namespace reconstruction { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -load_stencil(Array4 const& q_arr, const int idir, +load_stencil(amrex::Array4 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; @@ -47,9 +49,11 @@ load_stencil(Array4 const& q_arr, const int idir, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -load_passive_stencil(Array4 const& U_arr, Array4 const& rho_inv_arr, const int idir, +load_passive_stencil(amrex::Array4 const& U_arr, + amrex::Array4 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; @@ -78,5 +82,69 @@ load_passive_stencil(Array4 const& U_arr, Array4 const& } -#endif +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void +add_geometric_rho_source(amrex::Array4 const& q_arr, + amrex::Array4 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 const& q_arr, + amrex::Array4 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 const& q_arr, + amrex::Array4 const& qaux_arr, + amrex::Array4 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 diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 2595b9fe65..0795500825 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -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]; @@ -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++) { @@ -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]; @@ -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); } @@ -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); } @@ -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); } @@ -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 - }); } - -