add geometric terms for spherical 2D support. (#4141)

## Summary
This adds the appropriate geometric terms (area and volume) for
spherical 2D geometry. We assume dx[0]=dr, and dx[1]=d$`\theta`$.

## Additional background
This intended to deal with issue #3670

Src/Base/AMReX_CoordSys.cpp

@@ -164,7 +164,7 @@ CoordSys::UpperIndex(const Real* point) const noexcept
     IntVect ix;
     for (int k = 0; k < AMREX_SPACEDIM; k++)
     {
-        ix[k] = (int) ((point[k]-offset[k])/dx[k]);
+        ix[k] = (int) ((point[k]-offset[k])/dx[k]) + 1;
     }
     return ix;
 }
@@ -330,6 +330,8 @@ CoordSys::GetEdgeVolCoord (Vector<Real>& vc,
     GetEdgeLoc(vc,region,dir);
     //
     // In R direction of RZ, vol coord = (r^2)/2
+    // In R direction of SPHERICAL, vol coord = (r^3)/3
+    // In theta direction of SPHERICAL, vol coord = -cos(theta)
     //
 #if (AMREX_SPACEDIM == 2)
     if (dir == 0 && c_sys == RZ)
@@ -342,6 +344,29 @@ CoordSys::GetEdgeVolCoord (Vector<Real>& vc,
             vc[i] = 0.5_rt*r*r;
         }
     }
+    else if (c_sys == SPHERICAL)
+    {
+        if (dir == 0)
+        {
+            int len = static_cast<int>(vc.size());
+            AMREX_PRAGMA_SIMD
+            for (int i = 0; i < len; i++)
+            {
+                Real r = vc[i];
+                vc[i] = r*r*r/3.0_rt;
+            }
+        }
+        else
+        {
+            int len = static_cast<int>(vc.size());
+            AMREX_PRAGMA_SIMD
+            for (int i = 0; i < len; i++)
+            {
+                Real theta = vc[i];
+                vc[i] = -std::cos(theta);
+            }
+        }
+    }
 #elif (AMREX_SPACEDIM == 1)
     if (c_sys == SPHERICAL)
     {
@@ -365,8 +390,11 @@ CoordSys::GetCellVolCoord (Vector<Real>& vc,
     // are identical to physical distance from axis.
     //
     GetCellLoc(vc,region,dir);
+
     //
-    // In R direction of RZ, vol coord = (r^2)/2.
+    // In R direction of RZ, vol coord = (r^2)/2
+    // In R direction of SPHERICAL, vol coord = (r^3)/3
+    // In theta direction of SPHERICAL, vol coord = -cos(theta)
     //
 #if (AMREX_SPACEDIM == 2)
     if (dir == 0 && c_sys == RZ)
@@ -379,6 +407,29 @@ CoordSys::GetCellVolCoord (Vector<Real>& vc,
             vc[i] = 0.5_rt*r*r;
         }
     }
+    else if (c_sys == SPHERICAL)
+    {
+        if (dir == 0)
+        {
+            int len = static_cast<int>(vc.size());
+            AMREX_PRAGMA_SIMD
+            for (int i = 0; i < len; i++)
+            {
+                Real r = vc[i];
+                vc[i] = r*r*r/3.0_rt;
+            }
+        }
+        else
+        {
+            int len = static_cast<int>(vc.size());
+            AMREX_PRAGMA_SIMD
+            for (int i = 0; i < len; i++)
+            {
+                Real theta = vc[i];
+                vc[i] = -std::cos(theta);
+            }
+        }
+    }
 #elif (AMREX_SPACEDIM == 1)
     if (c_sys == SPHERICAL) {
         int len = static_cast<int>(vc.size());
@@ -462,6 +513,9 @@ CoordSys::Volume (const Real xlo[AMREX_SPACEDIM],
 #if (AMREX_SPACEDIM==2)
     case RZ:
         return static_cast<Real>(0.5*TWOPI)*(xhi[1]-xlo[1])*(xhi[0]*xhi[0]-xlo[0]*xlo[0]);
+    case SPHERICAL:
+        return static_cast<Real>(TWOPI/3.)*(std::cos(xlo[1])-std::cos(xhi[1])) *
+            (xhi[0]-xlo[0])*(xhi[0]*xhi[0]+xhi[0]*xlo[0]+xlo[0]*xlo[0]);
 #endif
     default:
         AMREX_ASSERT(0);
@@ -496,6 +550,16 @@ CoordSys::AreaLo (const IntVect& point, int dir) const noexcept // NOLINT(readab
             AMREX_ASSERT(0);
         }
         return 0._rt; // to silent compiler warning
+    case SPHERICAL:
+        LoNode(point,xlo);
+        switch (dir)
+        {
+        case 0: return Real(TWOPI)*xlo[0]*xlo[0]*(std::cos(xlo[1]) - std::cos(xlo[1]+dx[1]));
+        case 1: return (xlo[0]+xlo[0]+dx[0])*dx[0]*std::sin(xlo[1])*static_cast<Real>(0.5*TWOPI);
+        default:
+            AMREX_ASSERT(0);
+        }
+        return 0._rt; // to silent compiler warning
     default:
         AMREX_ASSERT(0);
     }
@@ -540,6 +604,16 @@ CoordSys::AreaHi (const IntVect& point, int dir) const noexcept // NOLINT(readab
             AMREX_ASSERT(0);
         }
         return 0._rt; // to silent compiler warning
+    case SPHERICAL:
+        HiNode(point,xhi);
+        switch (dir)
+        {
+        case 0: return Real(TWOPI)*xhi[0]*xhi[0]*(std::cos(xhi[1]-dx[1]) - std::cos(xhi[1]));
+        case 1: return (xhi[0]+xhi[0]-dx[0])*dx[0]*std::sin(xhi[1])*static_cast<Real>(0.5*TWOPI);
+        default:
+            AMREX_ASSERT(0);
+        }
+        return 0._rt; // to silent compiler warning
     default:
         AMREX_ASSERT(0);
     }

Src/Base/AMReX_Geometry.H

@@ -264,7 +264,7 @@ public:
 
             vol = dx[0] * dx[1];
         }
-        else {
+        else if (coord == CoordSys::RZ) {
             // Cylindrical
 
             Real r_l = geomdata.ProbLo()[0] + static_cast<Real>(point[0]) * dx[0];
@@ -273,6 +273,19 @@ public:
             constexpr Real pi = Real(3.1415926535897932);
             vol = pi * (r_l + r_r) * dx[0] * dx[1];
         }
+        else {
+            // Spherical
+
+            Real r_l = geomdata.ProbLo()[0] + static_cast<Real>(point[0]) * dx[0];
+            Real r_r = geomdata.ProbLo()[0] + static_cast<Real>(point[0]+1) * dx[0];
+
+            Real theta_l = geomdata.ProbLo()[1] + static_cast<Real>(point[1]) * dx[1];
+            Real theta_r = geomdata.ProbLo()[1] + static_cast<Real>(point[1]+1) * dx[1];
+
+            constexpr Real twoThirdsPi = static_cast<Real>(2.0 * 3.1415926535897932 / 3.0);
+            vol = twoThirdsPi * (std::cos(theta_l) - std::cos(theta_r)) * dx[0] *
+                (r_r*r_r + r_r*r_l + r_l*r_l);
+        }
 
 #else