Upgrade ExoCubed to modern C++ (#1)

Merge in csh's origin Cubed Sphere code

CODEOWNERS

@@ -17,3 +17,4 @@ src/exo3 @cshsgy
 src/harp @happysky19
 src/opacity @happysky19
 src/nbody @astrophy6
+patches/??.cs_*.patch @cshsgy

cmake/athena.cmake

@@ -8,7 +8,14 @@ set(patch_command
     ${CMAKE_CURRENT_SOURCE_DIR}/patches/23.meshblock.patch
     ${CMAKE_CURRENT_SOURCE_DIR}/patches/24.time_integrator.patch
     ${CMAKE_CURRENT_SOURCE_DIR}/patches/25.constant_acc.patch
-    ${CMAKE_CURRENT_SOURCE_DIR}/patches/26.scalars_flux.patch)
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/26.scalars_flux.patch
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/25.cs_meshblock.patch
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/26.cs_adiabatic_hydro.patch
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/27.cs_coordinates.patch
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/28.cs_bvals_cc.patch
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/29.cs_bvals_var.patch
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/30.cs_bvals_base.patch
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/31.cs_calculate_fluxes.patch)
 
 FetchContent_Declare(
   athenapp

cmake/parameters.cmake

@@ -55,6 +55,12 @@ else()
   set(COORDINATE_SYSTEM "affine_coordinate")
 endif()
 
+if(NOT CUBED_SPHERE OR NOT DEFINED CUBED_SPHERE)
+  set(CUBED_SPHERE_OPTION "NOT_CUBED_SPHERE")
+else()
+  set(CUBED_SPHERE_OPTION "CUBED_SPHERE")
+endif()
+
 if(NOT NETCDF OR NOT DEFINED NETCDF)
   set(NETCDF_OPTION "NO_NETCDFOUTPUT")
 else()

configure.hpp.in

@@ -19,6 +19,9 @@ enum {
 // Affine option (AFFINE or NOT_AFFINE)
 #define @AFFINE_OPTION@
 
+// Cubed Sphere option (CUBED_SPHERE or NOT_CUBED_SPHERE)
+#define @CUBED_SPHERE_OPTION@
+
 // NetCDF output (NETCDFOUTPUT or NO_NETCDFOUTPUT)
 #define @NETCDF_OPTION@
 

patches/25.cs_meshblock.patch

@@ -0,0 +1,26 @@
+diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp
+index 29a06d19..c465d93c 100644
+--- a/src/mesh/meshblock.cpp
++++ b/src/mesh/meshblock.cpp
+@@ -45,6 +45,10 @@
+ #include "mesh_refinement.hpp"
+ #include "meshblock_tree.hpp"
+
++// exo3 injection
++#include <exo3/affine_coordinate.hpp>
++#include <exo3/gnomonic_equiangle.hpp>
++
+ //----------------------------------------------------------------------------------------
+ //! MeshBlock constructor: constructs coordinate, boundary condition, hydro, field
+ //!                        and mesh refinement objects.
+@@ -124,6 +128,10 @@ MeshBlock::MeshBlock(int igid, int ilid, LogicalLocation iloc, RegionSize input_
+     pcoord = new KerrSchild(this, pin, false);
+   } else if (std::strcmp(COORDINATE_SYSTEM, "gr_user") == 0) {
+     pcoord = new GRUser(this, pin, false);
++  } else if (std::strcmp(COORDINATE_SYSTEM, "gnomonic_equiangle") == 0) {
++    pcoord = new GnomonicEquiangle(this, pin, false);
++  } else if (std::strcmp(COORDINATE_SYSTEM, "affine_coordinate") == 0) {
++    pcoord = new AffineCoordinate(this, pin, false);
+   }
+
+

patches/26.cs_adiabatic_hydro.patch

@@ -0,0 +1,53 @@
+diff --git a/src/eos/adiabatic_hydro.cpp b/src/eos/adiabatic_hydro.cpp
+index e61ff205..49aba323 100644
+--- a/src/eos/adiabatic_hydro.cpp
++++ b/src/eos/adiabatic_hydro.cpp
+@@ -19,6 +19,10 @@
+ #include "../mesh/mesh.hpp"
+ #include "../parameter_input.hpp"
+ #include "eos.hpp"
++
++#include <configure.hpp>
++#include <impl.hpp>
++#include <exo3/cubed_sphere.hpp>
+
+ // EquationOfState constructor
+
+@@ -67,7 +68,16 @@ void EquationOfState::ConservedToPrimitive(
+         w_vy = u_m2*di;
+         w_vz = u_m3*di;
+
+-        Real e_k = 0.5*di*(SQR(u_m1) + SQR(u_m2) + SQR(u_m3));
++        Real e_k;
++#ifdef CUBED_SPHERE
++        Real U, V;
++        // Convert m2 m3 to lat-lon
++        pmy_block_->pimpl->pexo3->GetUV(&U, &V, u_m2, u_m3, k, j, i);
++        e_k = 0.5*di*(SQR(u_m1) + SQR(U) + SQR(V));
++#else
++        e_k = 0.5*di*(SQR(u_m1) + SQR(u_m2) + SQR(u_m3));
++#endif
++
+         w_p = gm1*(u_e - e_k);
+
+         // apply pressure floor, correct total energy
+@@ -107,11 +117,18 @@ void EquationOfState::PrimitiveToConserved(
+         const Real& w_vz = prim(IVZ,k,j,i);
+         const Real& w_p  = prim(IPR,k,j,i);
+
++#ifdef CUBED_SPHERE
++        Real U, V;
++        // Convert vy vz to lat-lon
++        pmy_block_->pimpl->pexo3->GetUV(&U, &V, pco, w_vy, w_vz, k, j, i);
++        u_e = w_p*igm1 + 0.5*w_d*(SQR(w_vx) + SQR(U) + SQR(V));
++#else
++        u_e = w_p*igm1 + 0.5*w_d*(SQR(w_vx) + SQR(w_vy) + SQR(w_vz));
++#endif
+         u_d = w_d;
+         u_m1 = w_vx*w_d;
+         u_m2 = w_vy*w_d;
+         u_m3 = w_vz*w_d;
+-        u_e = w_p*igm1 + 0.5*w_d*(SQR(w_vx) + SQR(w_vy) + SQR(w_vz));
+       }
+     }
+   }

patches/27.cs_coordinates.patch

@@ -0,0 +1,53 @@
+diff --git a/src/coordinates/coordinates.cpp b/src/coordinates/coordinates.cpp
+index c11dd078..d594699e 100644
+--- a/src/coordinates/coordinates.cpp
++++ b/src/coordinates/coordinates.cpp
+@@ -20,6 +20,10 @@
+ #include "../nr_radiation/radiation.hpp"
+ #include "../parameter_input.hpp"
+ #include "coordinates.hpp"
++
++#include <configure.hpp>
++#include <impl.hpp>
++#include <exo3/cubed_sphere.hpp>
+
+ //----------------------------------------------------------------------------------------
+ //! Coordinates constructor: sets coordinates and coordinate spacing of cell FACES
+@@ -207,10 +209,18 @@ Coordinates::Coordinates(MeshBlock *pmb, ParameterInput *pin, bool flag) :
+           noffset = static_cast<std::int64_t>((j-jl)*2 + lx2*block_size.nx2);
+         }
+         Real rx = ComputeMeshGeneratorX(noffset, nrootmesh, true);
++#ifdef CUBED_SPHERE
++        x2f(j) = pmy_block->pimpl->pexo3->GenerateMeshX2(rx);
++#else
+         x2f(j) = pm->MeshGenerator_[X2DIR](rx, mesh_size);
++#endif
+       }
++#ifdef CUBED_SPHERE
++      // Do nothing
++#else
+       x2f(jl) = block_size.x2min;
+       x2f(ju+1) = block_size.x2max;
++#endif
+
+       for (int j=jl-ng; j<=ju+ng; ++j) {
+         dx2f(j) = dx;
+@@ -288,10 +298,18 @@ Coordinates::Coordinates(MeshBlock *pmb, ParameterInput *pin, bool flag) :
+           noffset = static_cast<std::int64_t>((k-kl)*2 + lx3*block_size.nx3);
+         }
+         Real rx = ComputeMeshGeneratorX(noffset, nrootmesh, true);
++#ifdef CUBED_SPHERE
++        x3f(k) = pmy_block->pimpl->pexo3->GenerateMeshX3(rx);
++#else
+         x3f(k) = pm->MeshGenerator_[X3DIR](rx, mesh_size);
++#endif
+       }
++#ifdef CUBED_SPHERE
++      // Do nothing
++#else
+       x3f(kl) = block_size.x3min;
+       x3f(ku+1) = block_size.x3max;
++#endif
+
+       for (int k=kl-ng; k<=ku+ng; ++k) {
+         dx3f(k) = dx;

patches/28.cs_bvals_cc.patch

@@ -0,0 +1,27 @@
+diff --git a/src/bvals/cc/bvals_cc.cpp b/src/bvals/cc/bvals_cc.cpp
+index eb6f6039..d8b328b0 100644
+--- a/src/bvals/cc/bvals_cc.cpp
++++ b/src/bvals/cc/bvals_cc.cpp
+@@ -30,6 +30,9 @@
+ #include "../../utils/buffer_utils.hpp"
+ #include "../bvals.hpp"
+ #include "bvals_cc.hpp"
++
++#include <configure.hpp>
++#include <exo3/cubed_sphere_utility.hpp>
+
+ // MPI header
+ #ifdef MPI_PARALLEL
+@@ -297,7 +298,12 @@ int CellCenteredBoundaryVariable::LoadBoundaryBufferSameLevel(Real *buf,
+   ek = (nb.ni.ox3 < 0) ? (pmb->ks + NGHOST - 1) : pmb->ke;
+   int p = 0;
+   AthenaArray<Real> &var = *var_cc;
++#ifdef CUBED_SPHERE
++// nl_, nu_, after buf var
++  CubedSphereUtility::PackData(var, buf, nl_, nu_, si, ei, sj, ej, sk, ek, p, nb.ni.ox1, nb.ni.ox2, nb.ni.ox3, pmb->loc);
++#else
+   BufferUtility::PackData(var, buf, nl_, nu_, si, ei, sj, ej, sk, ek, p);
++#endif
+   return p;
+ }
+

patches/29.cs_bvals_var.patch

@@ -0,0 +1,81 @@
+diff --git a/src/bvals/bvals_var.cpp b/src/bvals/bvals_var.cpp
+index 6d2460f1..5dac314b 100644
+--- a/src/bvals/bvals_var.cpp
++++ b/src/bvals/bvals_var.cpp
+@@ -22,6 +22,10 @@
+ #include "../mesh/mesh.hpp"
+ #include "bvals_interfaces.hpp"
+
++#include <configure.hpp>
++#include <impl.hpp>
++#include <exo3/cubed_sphere.hpp>
++
+ // MPI header
+ #ifdef MPI_PARALLEL
+ #include <mpi.h>
+@@ -211,6 +213,65 @@ void BoundaryVariable::SetCompletedFlagSameProcess(NeighborBlock& nb) {
+
+ void BoundaryVariable::SendBoundaryBuffers() {
+   MeshBlock *pmb = pmy_block_;
++#ifdef CUBED_SPHERE
++  // For cubed sphere, synchronize the in-panel boundary buffers first
++  // In cubed sphere option, MPI is automatically enabled and each process only allow one block
++  auto& loc = pmb->loc;
++  for (int n=0; n<pbval_->nneighbor; n++) {
++    NeighborBlock& nb = pbval_->neighbor[n];
++    int ox2 = nb.ni.ox2;
++    int ox3 = nb.ni.ox3;
++    int tox2, tox3;
++    pmb->pimpl->pexo3->TransformOX(&ox2, &ox3, &tox2, &tox3);
++    LogicalLocation tloc;
++    tloc.lx1 = loc.lx1;
++    tloc.lx2 = loc.lx2 + ox2;
++    tloc.lx3 = loc.lx3 + ox3;
++    int blockID_tg = CubedSphere::FindBlockID(tloc);
++    int blockID = CubedSphere::FindBlockID(loc);
++    if (blockID != blockID_tg) // Not in the same panel
++      continue;
++    if (bd_var_.sflag[nb.bufid] == BoundaryStatus::completed) continue;
++    int ssize;
++    ssize = LoadBoundaryBufferSameLevel(bd_var_.send[nb.bufid], nb); // Cubed sphere only has same-level neighbors
++#ifdef MPI_PARALLEL
++    MPI_Start(&(bd_var_.req_send[nb.bufid]));
++#endif
++    bd_var_.sflag[nb.bufid] = BoundaryStatus::completed;
++  }
++
++  // Receive the in-panel boundary buffers
++  // Note that the BoudaryStatus flag is set to "waiting" before this is called
++  for (int n=0; n<pbval_->nneighbor; n++) {
++    NeighborBlock& nb = pbval_->neighbor[n];
++    int ox2 = nb.ni.ox2;
++    int ox3 = nb.ni.ox3;
++    int tox2, tox3;
++    pmb->pimpl->pexo3->TransformOX(&ox2, &ox3, &tox2, &tox3, loc);
++    LogicalLocation tloc;
++    tloc.lx1 = loc.lx1;
++    tloc.lx2 = loc.lx2 + ox2;
++    tloc.lx3 = loc.lx3 + ox3;
++    int blockID_tg = CubedSphere::FindBlockID(tloc);
++    int blockID = CubedSphere::FindBlockID(loc);
++    if (blockID != blockID_tg) // Not in the same panel
++      continue;
++    if (bd_var_.flag[nb.bufid] == BoundaryStatus::arrived) continue;
++    while (bd_var_.flag[nb.bufid] == BoundaryStatus::waiting) {
++        int test;
++        // probe MPI communications.  This is a bit of black magic that seems to promote
++        // communications to top of stack and gets them to complete more quickly
++#ifdef MPI_PARALLEL
++        MPI_Iprobe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &test, MPI_STATUS_IGNORE);
++        MPI_Test(&(bd_var_.req_recv[nb.bufid]), &test, MPI_STATUS_IGNORE);
++        if (!static_cast<bool>(test)) {
++          continue;
++        }
++#endif
++        bd_var_.flag[nb.bufid] = BoundaryStatus::arrived;
++    }
++  }
++#endif
+   int mylevel = pmb->loc.level;
+   for (int n=0; n<pbval_->nneighbor; n++) {
+     NeighborBlock& nb = pbval_->neighbor[n];