add P_theta flux register and scale pressure using area and volume

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
 

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<std::unique_ptr<amrex::MultiFab> > rad_fluxes;
 #endif

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

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) {

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) {

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<Real const>{};
@@ -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<Real> 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<Real> ptheta_fab = ptheta.array();
+                Array4<Real> 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

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<Real> 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<Real> pradial_fab = pradial.array();
             Array4<Real> 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<Real> ptheta_fab = ptheta.array();
+            Array4<Real> 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<Real> pradial_fab = pradial.array();
             Array4<Real> 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<Real> ptheta_fab = ptheta.array();
+            Array4<Real> 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