add P_theta flux register and scale pressure using area and volume

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
 

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

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.

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

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

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

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<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) {
@@ -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<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
+
         }
 
 
@@ -708,6 +732,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;
 
@@ -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<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(0), 1, i, j, k, n,
+            {
+                P_theta_fab(i,j,k,0) += scale * ptheta_fab(i,j,k,0);
+            });
+
+          }
+#endif
         }
 
       } // MFIter loop