include pressure flux during transverse correction term

Source/hydro/trans.cpp

@@ -292,16 +292,28 @@ Castro::actual_trans_single(const Box& bx,  // NOLINT(readability-convert-member
         //
         // in cylindrical coords -- note that the p term is not
         // in a divergence for UMX in the x-direction, so there
-        // are no area factors.  For this geometry, we do not
+        // are no area factors.
+        //
+        // Similarly for Spherical2D geometry, we have:
+        // d(rho u)/dt + d(rho u v)/(rdtheta) = -1/r^2 d(r^2 rho u u)/dr - dp/dr
+        // d(rho v)/dt + d(rho u v)/dr = -1/(r sin(theta)) d(sin(theta) rho v v)/dtheta - 1/r dp/dtheta
+        //
+        // For these non-cartesian geometries, we do not
         // include p in our definition of the flux in the
-        // x-direction, for we need to fix this now.
+        // x-direction for Cylindrical2D or both x- and y-direction for spherical 2D
+        // So we need to fix this now.
+
         Real runewn = run - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMX) -
                                    area_t(il,jl,kl) * flux_t(il,jl,kl,UMX)) * volinv;
         if (idir_t == 0 && !mom_flux_has_p(0, idir_t, coord)) {
             runewn = runewn - cdtdx * (pgp - pgm);
         }
         Real rvnewn = rvn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMY) -
                                    area_t(il,jl,kl) * flux_t(il,jl,kl,UMY)) * volinv;
+        if (idir_t == 1 && !mom_flux_has_p(1, idir_t, coord)) {
+            Real r = geom.ProbLo(0) + static_cast<Real>(il + 0.5_rt) * geom.CellSize(0);
+            rvnewn = rvnewn - cdtdy / r * (pgp - pgm);
+        }
         Real rwnewn = rwn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMZ) -
                                    area_t(il,jl,kl) * flux_t(il,jl,kl,UMZ)) * volinv;
         Real renewn = ren - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UEDEN) -