diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index 54c364e26e..d951913be1 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.6940338039e-05 19441641.375 - xmom -5.4953770416e+14 1.3594264808e+14 - ymom -2.4933243003e+15 2.4933250525e+15 + density 8.6944011764e-05 19444371.155 + xmom -5.4917303392e+14 1.3509526313e+14 + ymom -2.4930150346e+15 2.4930157613e+15 zmom 0 0 - rho_E 7.4973602186e+11 5.0768248379e+24 - rho_e 7.1068648973e+11 5.0744783673e+24 - Temp 242282.60874 1404450633.1 - rho_He4 8.6940338039e-17 3.398107373 - rho_C12 3.4776135215e-05 7775850.9371 - rho_O16 5.2164202823e-05 11664450.012 - rho_Ne20 8.6940338039e-17 172485.537 - rho_Mg24 8.6940338039e-17 1043.054267 - rho_Si28 8.6940338039e-17 5.9869391361 - rho_S32 8.6940338039e-17 0.00016459247232 - rho_Ar36 8.6940338039e-17 1.9441643669e-05 - rho_Ca40 8.6940338039e-17 1.9441641397e-05 - rho_Ti44 8.6940338039e-17 1.9441641384e-05 - rho_Cr48 8.6940338039e-17 1.9441641384e-05 - rho_Fe52 8.6940338039e-17 1.9441641384e-05 - rho_Ni56 8.6940338039e-17 1.9441641384e-05 - phiGrav -5.870743119e+17 -2.337549858e+16 - grav_x -685044085.4 -51428.268861 - grav_y -739591083.78 739591039.26 + rho_E 7.4445974314e+11 5.0751963974e+24 + rho_e 7.0573031626e+11 5.0728530525e+24 + Temp 241944.6157 1402008261.1 + rho_He4 8.6944011764e-17 3.3044554839 + rho_C12 3.4777604705e-05 7776963.7255 + rho_O16 5.2166407058e-05 11666101.692 + rho_Ne20 8.6944011764e-17 167487.84482 + rho_Mg24 8.6944011764e-17 977.92440606 + rho_Si28 8.6944011764e-17 5.6766266345 + rho_S32 8.6944011764e-17 0.00015247530548 + rho_Ar36 8.6944011764e-17 1.9444373276e-05 + rho_Ca40 8.6944011764e-17 1.9444371178e-05 + rho_Ti44 8.6944011764e-17 1.9444371165e-05 + rho_Cr48 8.6944011764e-17 1.9444371164e-05 + rho_Fe52 8.6944011764e-17 1.9444371164e-05 + rho_Ni56 8.6944011764e-17 1.9444371164e-05 + phiGrav -5.8708220182e+17 -2.3375498709e+16 + grav_x -685053820.85 -51428.28215 + grav_y -739593455 739593428.31 grav_z 0 0 - rho_enuc -7.5385126121e+12 7.1503781318e+23 + rho_enuc -9.0785617027e+12 6.8835762986e+23 diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H index 0a77218ee6..741687b690 100644 --- a/Source/driver/Castro.H +++ b/Source/driver/Castro.H @@ -1301,6 +1301,9 @@ protected: #if (AMREX_SPACEDIM <= 2) amrex::MultiFab P_radial; #endif +#if (AMREX_SPACEDIM == 2) + amrex::MultiFab P_theta; +#endif #ifdef RADIATION amrex::Vector > rad_fluxes; #endif diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 0133c68bfe..064914bb98 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -895,6 +895,12 @@ Castro::initMFs() } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + P_theta.define(getEdgeBoxArray(1), dmap, 1, 0); + } +#endif + #ifdef RADIATION if (Radiation::rad_hydro_combined) { rad_fluxes.resize(AMREX_SPACEDIM); @@ -2599,6 +2605,12 @@ Castro::FluxRegCrseInit() { } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + fine_level.pres_reg.CrseInit(P_theta, 1, 0, 0, 1, pres_crse_scale); + } +#endif + #ifdef RADIATION if (Radiation::rad_hydro_combined) { for (int i = 0; i < AMREX_SPACEDIM; ++i) { @@ -2629,6 +2641,12 @@ Castro::FluxRegFineAdd() { } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + getLevel(level).pres_reg.FineAdd(P_theta, 1, 0, 0, 1, pres_fine_scale); + } +#endif + #ifdef RADIATION if (Radiation::rad_hydro_combined) { for (int i = 0; i < AMREX_SPACEDIM; ++i) { @@ -2866,12 +2884,9 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) reg = &getLevel(lev).pres_reg; - MultiFab dr(crse_lev.grids, crse_lev.dmap, 1, 0); - dr.setVal(crse_lev.geom.CellSize(0)); - reg->ClearInternalBorders(crse_lev.geom); - reg->Reflux(crse_state, dr, 1.0, 0, UMX, 1, crse_lev.geom); + reg->Reflux(crse_state, crse_lev.volume, 0, 1.0, 0, UMX, 1, crse_lev.geom); if (update_sources_after_reflux || !in_post_timestep) { @@ -2893,11 +2908,40 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) } +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + + reg->Reflux(crse_state, crse_lev.volume, 1, 1.0, 0, UMY, 1, crse_lev.geom); + + if (update_sources_after_reflux || !in_post_timestep) { + + MultiFab tmp_fluxes(crse_lev.P_theta.boxArray(), + crse_lev.P_theta.DistributionMap(), + crse_lev.P_theta.nComp(), crse_lev.P_theta.nGrow()); + + tmp_fluxes.setVal(0.0); + + for (OrientationIter fi; fi.isValid(); ++fi) + { + const FabSet& fs = (*reg)[fi()]; + if (fi().coordDir() == 1) { + fs.copyTo(tmp_fluxes, 0, 0, 0, tmp_fluxes.nComp()); + } + } + + MultiFab::Add(crse_lev.P_theta, tmp_fluxes, 0, 0, crse_lev.P_theta.nComp(), 0); + + } + + } +#endif + reg->setVal(0.0); } #endif + #ifdef RADIATION // This follows the same logic as the pure hydro fluxes; see above for details. diff --git a/Source/driver/Castro_advance.cpp b/Source/driver/Castro_advance.cpp index 6f9f0b9d07..b3eb3cdf59 100644 --- a/Source/driver/Castro_advance.cpp +++ b/Source/driver/Castro_advance.cpp @@ -551,6 +551,12 @@ Castro::initialize_advance(Real time, Real dt, int amr_iteration) } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + P_theta.setVal(0.0); + } +#endif + #ifdef RADIATION if (Radiation::rad_hydro_combined) { for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { diff --git a/Source/driver/Castro_advance_ctu.cpp b/Source/driver/Castro_advance_ctu.cpp index af387a92d2..d2fe4fa1da 100644 --- a/Source/driver/Castro_advance_ctu.cpp +++ b/Source/driver/Castro_advance_ctu.cpp @@ -204,6 +204,12 @@ Castro::retry_advance_ctu(Real dt, const advance_status& status) } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + P_theta.setVal(0.0); + } +#endif + #ifdef RADIATION if (Radiation::rad_hydro_combined) { for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 9791280851..c9f5af5bcc 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -179,6 +179,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co #if AMREX_SPACEDIM <= 2 FArrayBox pradial(The_Async_Arena()); #endif +#if AMREX_SPACEDIM == 2 + FArrayBox ptheta(The_Async_Arena()); +#endif #if AMREX_SPACEDIM == 3 FArrayBox qmyx(The_Async_Arena()), qpyx(The_Async_Arena()); FArrayBox qmzx(The_Async_Arena()), qpzx(The_Async_Arena()); @@ -456,6 +459,13 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co fab_size += pradial.nBytes(); #endif +#if AMREX_SPACEDIM == 2 + if (Geom().IsSPHERICAL()) { + ptheta.resize(ybx, 1); + } + fab_size += ptheta.nBytes(); +#endif + #ifdef SIMPLIFIED_SDC #ifdef REACTIONS Array4 const sdc_src_arr = SDC_react_source.array(mfi); @@ -1287,10 +1297,25 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - pradial_fab(i,j,k) = qex_arr(i,j,k,GDPRES) * dt; + pradial_fab(i,j,k) = area_arr(i,j,k) * qex_arr(i,j,k,GDPRES) * dt; }); } #endif + +#if AMREX_SPACEDIM == 2 + // get the scaled pressure in the theta direction + + if (idir == 1 && !mom_flux_has_p(1, 1, coord)) { + Array4 ptheta_fab = ptheta.array(); + + amrex::ParallelFor(nbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + ptheta_fab(i,j,k) = area_arr(i,j,k) * qey_arr(i,j,k,GDPRES) * dt; + }); + } +#endif + // Store the fluxes from this advance. For simplified SDC integration we // only need to do this on the last iteration. @@ -1334,7 +1359,19 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co P_radial_fab(i,j,k,0) += pradial_fab(i,j,k,0); }); } +#endif + +#if AMREX_SPACEDIM == 2 + if (idir == 1 && !mom_flux_has_p(1, 1, coord)) { + Array4 ptheta_fab = ptheta.array(); + Array4 P_theta_fab = P_theta.array(mfi); + amrex::ParallelFor(mfi.nodaltilebox(0), + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + P_theta_fab(i,j,k,0) += ptheta_fab(i,j,k,0); + }); + } #endif } // add_fluxes diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 8b00cd958b..0772b71b79 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -84,6 +84,9 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) } #if AMREX_SPACEDIM <= 2 FArrayBox pradial(The_Async_Arena()); +#endif +#if AMREX_SPACEDIM == 2 + FArrayBox ptheta(The_Async_Arena()); #endif FArrayBox avis(The_Async_Arena()); @@ -654,8 +657,12 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) if (!Geom().IsCartesian()) { pradial.resize(xbx, 1); } +#endif - Array4 pradial_fab = pradial.array(); +#if AMREX_SPACEDIM == 2 + if (Geom().IsSPHERICAL()) { + ptheta.resize(ybx, 1); + } #endif for (int idir = 0; idir < AMREX_SPACEDIM; ++idir) { @@ -671,15 +678,32 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) // get the scaled radial pressure -- we need to treat this specially if (idir == 0 && !mom_flux_has_p(0, 0, coord)) { + Array4 pradial_fab = pradial.array(); Array4 const qex_arr = qe[idir].array(); amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - pradial_fab(i,j,k) = qex_arr(i,j,k,GDPRES) * dt; + pradial_fab(i,j,k) = area_arr(i,j,k) * qex_arr(i,j,k,GDPRES) * dt; }); + } #endif + +#if AMREX_SPACEDIM == 2 + // get the scaled pressure in the theta direction + + if (idir == 1 && !mom_flux_has_p(1, 1, coord)) { + Array4 ptheta_fab = ptheta.array(); + Array4 const qey_arr = qe[idir].array(); + + amrex::ParallelFor(nbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + ptheta_fab(i,j,k) = area_arr(i,j,k) * qey_arr(i,j,k,GDPRES) * dt; + }); } +#endif + } @@ -708,6 +732,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) #if AMREX_SPACEDIM <= 2 if (!Geom().IsCartesian()) { + Array4 pradial_fab = pradial.array(); Array4 P_radial_fab = P_radial.array(mfi); const Real scale = stage_weight; @@ -718,6 +743,21 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) } #endif + +#if AMREX_SPACEDIM == 2 + if (Geom().IsSPHERICAL()) { + + Array4 ptheta_fab = ptheta.array(); + Array4 P_theta_fab = P_theta.array(mfi); + const Real scale = stage_weight; + + AMREX_HOST_DEVICE_FOR_4D(mfi.nodaltilebox(0), 1, i, j, k, n, + { + P_theta_fab(i,j,k,0) += scale * ptheta_fab(i,j,k,0); + }); + + } +#endif } } // MFIter loop