From 028638564f7be0694b9898f8d4088cdbf9a6f9f5 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Thu, 12 Sep 2024 17:41:36 -0400 Subject: [PATCH] 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 | 78 ++++++++++++++++++++++++++++++++++++- Src/Base/AMReX_Geometry.H | 15 ++++++- 2 files changed, 90 insertions(+), 3 deletions(-) diff --git a/Src/Base/AMReX_CoordSys.cpp b/Src/Base/AMReX_CoordSys.cpp index 757af532cc5..a958d33a5b6 100644 --- a/Src/Base/AMReX_CoordSys.cpp +++ b/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& 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& vc, vc[i] = 0.5_rt*r*r; } } + else if (c_sys == SPHERICAL) + { + if (dir == 0) + { + int len = static_cast(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(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& 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& vc, vc[i] = 0.5_rt*r*r; } } + else if (c_sys == SPHERICAL) + { + if (dir == 0) + { + int len = static_cast(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(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(vc.size()); @@ -462,6 +513,9 @@ CoordSys::Volume (const Real xlo[AMREX_SPACEDIM], #if (AMREX_SPACEDIM==2) case RZ: return static_cast(0.5*TWOPI)*(xhi[1]-xlo[1])*(xhi[0]*xhi[0]-xlo[0]*xlo[0]); + case SPHERICAL: + return static_cast(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(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(0.5*TWOPI); + default: + AMREX_ASSERT(0); + } + return 0._rt; // to silent compiler warning default: AMREX_ASSERT(0); } diff --git a/Src/Base/AMReX_Geometry.H b/Src/Base/AMReX_Geometry.H index 6231fbd1f9b..1b7bae1a7d3 100644 --- a/Src/Base/AMReX_Geometry.H +++ b/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(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(point[0]) * dx[0]; + Real r_r = geomdata.ProbLo()[0] + static_cast(point[0]+1) * dx[0]; + + Real theta_l = geomdata.ProbLo()[1] + static_cast(point[1]) * dx[1]; + Real theta_r = geomdata.ProbLo()[1] + static_cast(point[1]+1) * dx[1]; + + constexpr Real twoThirdsPi = static_cast(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