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/driver/timestep.cpp b/Source/driver/timestep.cpp index a7b5936e4e..fd7e029a4b 100644 --- a/Source/driver/timestep.cpp +++ b/Source/driver/timestep.cpp @@ -35,6 +35,9 @@ Castro::estdt_cfl (int is_new) // Courant-condition limited timestep const auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); + const auto coord = geom.Coord(); + amrex::ignore_unused(problo, coord); const MultiFab& stateMF = is_new ? get_new_data(State_Type) : get_old_data(State_Type); @@ -82,6 +85,11 @@ Castro::estdt_cfl (int is_new) Real dt2; #if AMREX_SPACEDIM >= 2 dt2 = dx[1]/(c + std::abs(uy)); + if (coord == 2) { + // dx[1] in Spherical2D is just dtheta, need rdtheta for physical length + // so just multiply by the smallest r + dt2 *= problo[0] + 0.5_rt * dx[0]; + } #else dt2 = dt1; #endif @@ -127,6 +135,9 @@ Castro::estdt_mhd (int is_new) // MHD timestep limiter const auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); + const auto coord = geom.Coord(); + amrex::ignore_unused(problo, coord); const MultiFab& U_state = is_new ? get_new_data(State_Type) : get_old_data(State_Type); @@ -206,6 +217,9 @@ Castro::estdt_mhd (int is_new) Real dt2; #if AMREX_SPACEDIM >= 2 dt2 = dx[1]/(cy + std::abs(uy)); + if (coord == 2) { + dt2 *= problo[0] + 0.5_rt * dx[0]; + } #else dt2 = dt1; #endif @@ -238,6 +252,9 @@ Castro::estdt_temp_diffusion (int is_new) // where D = k/(rho c_v), and k is the conductivity const auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); + const auto coord = geom.Coord(); + amrex::ignore_unused(problo, coord); const MultiFab& stateMF = is_new ? get_new_data(State_Type) : get_old_data(State_Type); @@ -287,6 +304,10 @@ Castro::estdt_temp_diffusion (int is_new) Real dt2; #if AMREX_SPACEDIM >= 2 dt2 = 0.5_rt * dx[1]*dx[1] / D; + if (coord == 2) { + Real r = problo[0] + 0.5_rt * dx[0]; + dt2 *= r * r; + } #else dt2 = dt1; #endif @@ -320,6 +341,9 @@ Castro::estdt_burning (int is_new) } const auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); + const auto coord = geom.Coord(); + amrex::ignore_unused(problo, coord); MultiFab& stateMF = is_new ? get_new_data(State_Type) : get_old_data(State_Type); @@ -368,7 +392,13 @@ Castro::estdt_burning (int is_new) #if AMREX_SPACEDIM == 1 burn_state.dx = dx[0]; #else - burn_state.dx = amrex::min(AMREX_D_DECL(dx[0], dx[1], dx[2])); + Real dx1 = dx[1]; +#if AMREX_SPACEDIM >= 2 + if (coord == 2) { + dx1 *= problo[0] + 0.5_rt * dx[0]; + } +#endif + burn_state.dx = amrex::min(AMREX_D_DECL(dx[0], dx1, dx[2])); #endif burn_state.rho = S(i,j,k,URHO); @@ -464,6 +494,9 @@ Real Castro::estdt_rad (int is_new) { auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); + const auto coord = geom.Coord(); + amrex::ignore_unused(problo, coord); const MultiFab& stateMF = is_new ? get_new_data(State_Type) : get_old_data(State_Type); const MultiFab& radMF = is_new ? get_new_data(Rad_Type) : get_old_data(Rad_Type); @@ -523,6 +556,9 @@ Castro::estdt_rad (int is_new) Real dt1 = dx[0] / (c + std::abs(ux)); #if AMREX_SPACEDIM >= 2 Real dt2 = dx[1] / (c + std::abs(uy)); + if (coord == 2) { + dt2 *= problo[0] + 0.5_rt * dx[0]; + } #else Real dt2 = std::numeric_limits::max(); #endif diff --git a/Source/hydro/Castro_hydro.cpp b/Source/hydro/Castro_hydro.cpp index 0dedbc17aa..0e2250ed20 100644 --- a/Source/hydro/Castro_hydro.cpp +++ b/Source/hydro/Castro_hydro.cpp @@ -237,12 +237,20 @@ Castro::check_for_cfl_violation(const MultiFab& State, const Real dt) int cfl_violation = 0; auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); + const auto coord = geom.Coord(); + amrex::ignore_unused(problo, coord); Real dtdx = dt / dx[0]; Real dtdy = 0.0_rt; if (AMREX_SPACEDIM >= 2) { dtdy = dt / dx[1]; + if (coord == 2) { + // dx[1] in Spherical2D is just rdtheta, need rdtheta for physical length + // Just choose to divide by the smallest r + dtdy /= problo[0] + 0.5_rt * dx[0]; + } } Real dtdz = 0.0_rt; diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index 869ec477b8..e17e4d0326 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -113,7 +113,6 @@ Castro::shock(const Box& bx, #endif #endif -#if AMREX_SPACEDIM == 1 } else if (coord_type == 2) { // 1-d spherical @@ -122,6 +121,15 @@ Castro::shock(const Box& bx, Real rp = (i + 1 + 0.5_rt) * dx[0]; div_u += 0.5_rt * (rp * rp * q_arr(i+1,j,k,QU) - rm * rm * q_arr(i-1,j,k,QU)) / (rc * rc * dx[0]); +#if AMREX_SPACEDIM == 2 + + Real thetac = (j + 0.5_rt) * dx[1]; + Real thetam = (j - 1 + 0.5_rt) * dx[1]; + Real thetap = (j + 1 + 0.5_rt) * dx[1]; + + div_u += 0.5_rt * (std::sin(thetap) * q_arr(i,j+1,k,QV) - + std::sin(thetam) * q_arr(i,j-1,k,QV)) / + (rc * sin(thetac) * dx[1]); #endif #ifndef AMREX_USE_GPU @@ -134,7 +142,14 @@ Castro::shock(const Box& bx, // now compute (grad P - rho g) . dx // We subtract off the hydrostatic force, since the pressure that - // balances that is not available to make a shock. + // balances that is not available to make a shock. We compute this + // as: + // + // P'_{i+1} = P_{i+1} - [ P_i + \int_{x_i}^{x_{i+1}} rho g dx ] + // + // where the term in the [ ] is the hydrostatic pressure in i+1 + // computed by integrating from x_i to x_{i+1} + // // We'll use a centered diff for the pressure gradient. Real dP_x = 0.5_rt * (q_arr(i+1,j,k,QPRES) - q_arr(i-1,j,k,QPRES)); if (shock_detection_include_sources == 1) { @@ -145,7 +160,13 @@ Castro::shock(const Box& bx, #if AMREX_SPACEDIM >= 2 dP_y = 0.5_rt * (q_arr(i,j+1,k,QPRES) - q_arr(i,j-1,k,QPRES)); if (shock_detection_include_sources == 1) { - dP_y += -0.25_rt * dx[1] * (U_src_arr(i,j+1,k,UMY) + 2.0_rt * U_src_arr(i,j,k,UMY) + U_src_arr(i,j-1,k,UMY)); + Real dy{dx[1]}; + if (coord_type == 2) { + // dx[1] is just dtheta + Real rc = (i + 0.5_rt) * dx[0]; + dy *= rc; + } + dP_y += -0.25_rt * dy * (U_src_arr(i,j+1,k,UMY) + 2.0_rt * U_src_arr(i,j,k,UMY) + U_src_arr(i,j-1,k,UMY)); } #endif #if AMREX_SPACEDIM == 3 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 70ea810611..4387f9f7c6 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]; @@ -297,11 +297,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); } @@ -334,11 +343,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); } @@ -351,11 +369,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); } @@ -627,34 +654,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) = std::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) = std::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 - }); } diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index acbc4947df..18dd50c2e0 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -101,6 +101,8 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member #if AMREX_SPACEDIM == 2 int coord = geom.Coord(); + const auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); #endif bool reset_density = transverse_reset_density; @@ -292,9 +294,17 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member // // in cylindrical coords -- note that the p term is not // in a divergence for UMX in the x-direction, so there - // are no area factors. For this geometry, we do not + // are no area factors. + // + // Similarly for Spherical2D geometry, we have: + // d(rho u)/dt + d(rho u v)/(rdtheta) = -1/r^2 d(r^2 rho u u)/dr - dp/dr + // d(rho v)/dt + d(rho u v)/dr = -1/(r sin(theta)) d(sin(theta) rho v v)/dtheta - 1/r dp/dtheta + // + // For these non-cartesian geometries, we do not // include p in our definition of the flux in the - // x-direction, for we need to fix this now. + // x-direction for Cylindrical2D or both x- and y-direction for spherical 2D + // So we need to fix this now. + Real runewn = run - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMX) - area_t(il,jl,kl) * flux_t(il,jl,kl,UMX)) * volinv; if (idir_t == 0 && !mom_flux_has_p(0, idir_t, coord)) { @@ -302,6 +312,10 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member } Real rvnewn = rvn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMY) - area_t(il,jl,kl) * flux_t(il,jl,kl,UMY)) * volinv; + if (idir_t == 1 && !mom_flux_has_p(1, idir_t, coord)) { + Real r = problo[0] + static_cast(il + 0.5_rt) * dx[0]; + rvnewn = rvnewn - cdtdx / r * (pgp - pgm); + } Real rwnewn = rwn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMZ) - area_t(il,jl,kl) * flux_t(il,jl,kl,UMZ)) * volinv; Real renewn = ren - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UEDEN) -