update mom_flux_has_p to address with 2D spherical (#2958)

This is to address issue #2927 . I updated mom_flux_has_p to work with theta direction in 2d spherical. Also update places where it is used. I think it should just work in riemann_solver and HLL_solver, so nothing is changed there. 

For Castro_ctu_hydro.cpp and Castro_mol_hydro.cpp I will need to include ptheta just like pradial, so I'll do that in a separate pr to address that separately. I'll also address trans.cpp separately when I do the transverse flux pr

Source/driver/Castro_util.H

@@ -71,9 +71,15 @@ bool mom_flux_has_p (const int mom_dir, const int flux_dir, const int coord)
 
         if (coord != CoordSys::cartesian) {
             return false;
-        } else {
-            return true;
         }
+        return true;
+
+    } else if (mom_dir == 1 && flux_dir == 1) {
+
+        if (coord == CoordSys::SPHERICAL) {
+            return false;
+        }
+        return true;
 
     } else {
         return (mom_dir == flux_dir);

Source/hydro/Castro_ctu.cpp

@@ -42,7 +42,7 @@ Castro::consup_hydro(const Box& bx,
 #endif
         ) * volinv;
 
-   // Add the p div(u) source term to (rho e).
+    // Add the p div(u) source term to (rho e).
     if (n == UEINT) {
 
       Real pdu = (qx(i+1,j,k,GDPRES) + qx(i,j,k,GDPRES)) *
@@ -65,13 +65,19 @@ Castro::consup_hydro(const Box& bx,
 
       U_new(i,j,k,n) = U_new(i,j,k,n) - dt * pdu;
 
-    } else if (n == UMX) {
-      // Add gradp term to momentum equation -- only for axisymmetric
-      // coords (and only for the radial flux).
+    } else if (n == UMX && !mom_flux_has_p(0, 0, geomdata.Coord())) {
+        // Add gradp term to momentum equation -- only for axisymmetric
+        // coords (and only for the radial flux).
 
-      if (!mom_flux_has_p(0, 0, geomdata.Coord())) {
-        U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize()[0];
-      }
+        U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize(0);
+
+#if AMREX_SPACEDIM >= 2
+    } else if (n == UMY && !mom_flux_has_p(1, 1, geomdata.Coord())) {
+        // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry
+
+        Real r = geomdata.ProbLo(0) + (static_cast<Real>(i) + 0.5_rt) * geomdata.CellSize(0);
+        U_new(i,j,k,UMY) += - dt * (qy(i,j+1,k,GDPRES) - qy(i,j,k,GDPRES)) / (r * geomdata.CellSize(1));
+#endif
     }
   });
 }

Source/hydro/Castro_hydro.H

@@ -387,6 +387,7 @@
 /// @param area1   area of y faces
 /// @param area2   area of z faces
 /// @param q0      Godunuv state on x faces
+/// @param q1      Godunuv state on y faces
 /// @param vol     cell volume
 ///
     void mol_consup(const amrex::Box& bx,
@@ -412,6 +413,7 @@
 #endif
 #if AMREX_SPACEDIM <= 2
                     amrex::Array4<amrex::Real const> const& q0,
+                    amrex::Array4<amrex::Real const> const& q1,
 #endif
                     amrex::Array4<amrex::Real const> const& vol);
 

Source/hydro/Castro_mol.cpp

@@ -226,6 +226,7 @@ Castro::mol_consup(const Box& bx,  // NOLINT(readability-convert-member-function
 #endif
 #if AMREX_SPACEDIM <= 2
                    Array4<Real const> const& q0,
+                   Array4<Real const> const& q1,
 #endif
                    Array4<Real const> const& vol) {
 
@@ -238,6 +239,7 @@ Castro::mol_consup(const Box& bx,  // NOLINT(readability-convert-member-function
 #if AMREX_SPACEDIM <= 2
   const auto dx = geom.CellSizeArray();
   auto coord = geom.Coord();
+  auto prob_lo = geom.ProbLoArray();
 #endif
 
   amrex::ParallelFor(bx, NUM_STATE,
@@ -258,12 +260,18 @@ Castro::mol_consup(const Box& bx,  // NOLINT(readability-convert-member-function
 #endif
 
 #if AMREX_SPACEDIM <= 2
-    if (n == UMX && do_hydro == 1) {
-        // Add gradp term to momentum equation -- only for axisymmetric
-        // coords (and only for the radial flux).
+    if (do_hydro == 1) {
+        if (n == UMX && !mom_flux_has_p(0, 0, coord)) {
+            // Add gradp term to radial momentum equation -- only for axisymmetric
+            // coords.
 
-        if (!mom_flux_has_p(0, 0, coord)) {
             update(i,j,k,UMX) -= (q0(i+1,j,k,GDPRES) - q0(i,j,k,GDPRES)) / dx[0];
+
+        } else if (n == UMY && !mom_flux_has_p(1, 1, coord)) {
+            // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry
+
+            Real r = prob_lo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0];
+            update(i,j,k,UMY) -= (q1(i,j+1,k,GDPRES) - q1(i,j,k,GDPRES)) / (r * dx[1]);
         }
     }
 #endif

Source/hydro/Castro_mol_hydro.cpp

@@ -620,6 +620,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update)
 #endif
 #if AMREX_SPACEDIM <= 2
                    qe[0].array(),
+                   qe[1].array(),
 #endif
                    volume.array(mfi));