diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index 86ff0bb723..1cf5235930 100644 --- a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out +++ b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out @@ -1,30 +1,30 @@ plotfile = plt00086 time = 1.25 variables minimum value maximum value - density 8.6940338039e-05 19441641.404 - xmom -5.4953770558e+14 1.3594264808e+14 - ymom -2.4933244206e+15 2.4933251729e+15 + density 8.7357118639e-05 19430545.7 + xmom -5.4899611972e+14 1.3502124975e+14 + ymom -2.4942903676e+15 2.4942903672e+15 zmom 0 0 - rho_E 7.4973602186e+11 5.0768249302e+24 - rho_e 7.1068648972e+11 5.0744784597e+24 - Temp 242282.60874 1404450065.3 - rho_He4 8.6940338039e-17 3.3980876296 - rho_C12 3.4776135215e-05 7775851.0286 - rho_O16 5.2164202823e-05 11664450.983 - rho_Ne20 8.6940338039e-17 172486.25459 - rho_Mg24 8.6940338039e-17 1043.0854443 - rho_Si28 8.6940338039e-17 5.9870930966 - rho_S32 8.6940338039e-17 0.00016460178455 - rho_Ar36 8.6940338039e-17 1.9441643904e-05 - rho_Ca40 8.6940338039e-17 1.9441641599e-05 - rho_Ti44 8.6940338039e-17 1.9441641468e-05 - rho_Cr48 8.6940338039e-17 1.9441641465e-05 - rho_Fe52 8.6940338039e-17 1.9441641465e-05 - rho_Ni56 8.6940338039e-17 1.9441641465e-05 + rho_E 1.1708526936e+12 5.0796765656e+24 + rho_e 1.1071097588e+12 5.0773720951e+24 + Temp 265286.90326 1411724449.2 + rho_He4 8.7357118639e-17 3.721331829 + rho_C12 3.4942847455e-05 7771349.199 + rho_O16 5.2414271183e-05 11657751.712 + rho_Ne20 8.7357118639e-17 189502.87202 + rho_Mg24 8.7357118639e-17 1284.4110655 + rho_Si28 8.7357118639e-17 7.1480863104 + rho_S32 8.7357118639e-17 0.0002130320139 + rho_Ar36 8.7357118639e-17 1.9430548918e-05 + rho_Ca40 8.7357118639e-17 1.9430545919e-05 + rho_Ti44 8.7357118639e-17 1.9430545769e-05 + rho_Cr48 8.7357118639e-17 1.9430545765e-05 + rho_Fe52 8.7357118639e-17 1.9430545765e-05 + rho_Ni56 8.7357118639e-17 1.9430545765e-05 Shock 0 1 - phiGrav -5.8707431189e+17 -2.337549858e+16 - grav_x -685044085.35 -51428.268861 - grav_y -739591083.9 739591039.39 + phiGrav -5.870522222e+17 -2.3375498723e+16 + grav_x -685017818.43 -51428.283425 + grav_y -739492913.39 739492913.53 grav_z 0 0 - rho_enuc 0 7.1503919723e+23 + rho_enuc 0 8.1320260393e+23 diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H index 64c664db9a..63c8ead297 100644 --- a/Source/driver/Castro.H +++ b/Source/driver/Castro.H @@ -1282,6 +1282,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 e2b50f3a8c..ef74c14d48 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -898,6 +898,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); @@ -2619,6 +2625,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) { @@ -2649,6 +2661,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) { @@ -2884,14 +2902,22 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) #if (AMREX_SPACEDIM <= 2) if (!Geom().IsCartesian()) { + // Get pressure flux register of this level. + reg = &getLevel(lev).pres_reg; - MultiFab dr(crse_lev.grids, crse_lev.dmap, 1, 0); - dr.setVal(crse_lev.geom.CellSize(0)); + // Clear out flux at the internal borders, + // i.e. not the borders between coarse and fine boxes. reg->ClearInternalBorders(crse_lev.geom); - reg->Reflux(crse_state, dr, 1.0, 0, UMX, 1, crse_lev.geom); + // Perform the reflux + // i.e. U^{c, new} = U^{c, old} + (Σ p^f A^f dt^f - p^c A^c dt^c) / V^c + // Note that the pressure flux register holds what's inside the parenthesis + // And this is only done for U = UMX for radial pressure flux register + // which is stored in the 0-dir of pres_reg. + + reg->Reflux(crse_state, crse_lev.volume, 0, 1.0, 0, UMX, 1, crse_lev.geom); if (update_sources_after_reflux || !in_post_timestep) { @@ -2913,6 +2939,39 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) } +#if (AMREX_SPACEDIM == 2) + // Now deal with theta pressure flux register with 2d spherical geometry + + if (Geom().IsSPHERICAL()) { + + // Do reflux, but note theta pressure flux register is stored + // in the 1-dir of pres_reg. And it is only applied for U=UMY. + + 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); } diff --git a/Source/driver/Castro_advance.cpp b/Source/driver/Castro_advance.cpp index 482d8a1876..525f61fe1c 100644 --- a/Source/driver/Castro_advance.cpp +++ b/Source/driver/Castro_advance.cpp @@ -574,6 +574,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 af1c4e1671..c7a2dabd6f 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()) { + getLevel(lev).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 dddfb6a2ef..6fab90914e 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 REACTIONS auto sdc_src_arr = castro::time_integration_method == SimplifiedSpectralDeferredCorrections ? SDC_react_source.array(mfi) : Array4{}; @@ -1275,10 +1285,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. @@ -1322,7 +1347,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(1), + [=] 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 c50f1154c7..a74c84540f 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()); @@ -649,8 +652,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) { @@ -666,15 +673,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 + } @@ -703,6 +727,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; @@ -713,6 +738,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(1), 1, i, j, k, n, + { + P_theta_fab(i,j,k,0) += scale * ptheta_fab(i,j,k,0); + }); + + } +#endif } } // MFIter loop