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 - }); } - -