Fixing Cubed Sphere Scalar Tracer (#164)

Mac/CI passes Linux/CI fails on bryan --------- Co-authored-by: mac/cli <[email protected]> Co-authored-by: Sihe Chen <[email protected]>
chengcli · Jul 12, 2024 · 1b9b64f · 1b9b64f
1 parent 5c04146
commit 1b9b64f
Show file tree

Hide file tree

Showing 12 changed files with 135 additions and 53 deletions.
diff --git a/cmake/athena.cmake b/cmake/athena.cmake
@@ -11,10 +11,12 @@ set(patch_command
     ${CMAKE_CURRENT_SOURCE_DIR}/patches/26.scalars_flux.patch
     ${CMAKE_CURRENT_SOURCE_DIR}/patches/26.adiabatic_hydro.patch
     ${CMAKE_CURRENT_SOURCE_DIR}/patches/27.coordinates.patch
-    ${CMAKE_CURRENT_SOURCE_DIR}/patches/28.bvals_cc.patch
     ${CMAKE_CURRENT_SOURCE_DIR}/patches/29.mesh.patch
     ${CMAKE_CURRENT_SOURCE_DIR}/patches/30.bvals_base.patch
-    ${CMAKE_CURRENT_SOURCE_DIR}/patches/31.task_list.patch)
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/31.task_list.patch
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/32.scalars.patch
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/33.bvals_cc_cpp.patch
+    ${CMAKE_CURRENT_SOURCE_DIR}/patches/34.bvals_cc_hpp.patch)
 
 FetchContent_Declare(
   athenapp

diff --git a/cmake/examples/exo2.cmake b/cmake/examples/exo2.cmake
@@ -8,6 +8,7 @@ endmacro()
 
 # athena variables
 set_if_empty(NUMBER_GHOST_CELLS 3)
+set_if_empty (NTRACER 3)
 
 # canoe configure
 set(CUBED_SPHERE ON)

diff --git a/cmake/examples/exo3.cmake b/cmake/examples/exo3.cmake
@@ -15,3 +15,5 @@ set(NETCDF ON)
 set(MPI ON)
 set(PNETCDF ON)
 set(RSOLVER hllc_transform)
+
+set (NTRACER 3)
diff --git a/cmake/examples/hjupiter.cmake b/cmake/examples/hjupiter.cmake
@@ -15,6 +15,7 @@ set(PNETCDF ON)
 set(MPI ON)
 set(DISORT ON)
 set(PYTHON_BINDINGS ON)
+set(NCHEMISTRY 1)
 set_if_empty(RSOLVER hllc_transform)
 # set(GLOG ON)
 

diff --git a/examples/2023-Chen-exo3/steady_zonal.cpp b/examples/2023-Chen-exo3/steady_zonal.cpp
@@ -4,6 +4,7 @@
 #include <athena/hydro/hydro.hpp>
 #include <athena/mesh/mesh.hpp>
 #include <athena/parameter_input.hpp>
+#include <athena/scalars/scalars.hpp>
 
 // application
 #include <application/application.hpp>
@@ -80,6 +81,11 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
         pexo3->GetLatLon(&lat, &lon, k, j, i);
         phydro->w(IDN, k, j, i) =
             g * h0 - (a * om_earth * u0 + u0 * u0 / 2) * sin(lat) * sin(lat);
+        if ((lat > PI / 2.0 * 0.9) && (lat < PI / 2.0))
+            pscalars->s(2, k, j, i) = 1.0;
+        if ((lat > PI / 2.0 * 0.5) && (lat < PI / 2.0 * 0.6))
+            pscalars->s(1, k, j, i) = 1.0;
+
         Real U = u0 * cos(lat);
         Real V = 0.0;
         Real Vy, Vz;
@@ -118,16 +124,20 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
         user_out_var(2, k, j, i) = U;
         user_out_var(3, k, j, i) = V;
         user_out_var(4, k, j, i) = delta / (C * D);
+        user_out_var(5, k, j, i) = pscalars->s(1, k, j, i);
+        user_out_var(6, k, j, i) = pscalars->s(2, k, j, i);
       }
 }
 
 void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
-  AllocateUserOutputVariables(5);
+  AllocateUserOutputVariables(7);
   SetUserOutputVariableName(0, "lat");
   SetUserOutputVariableName(1, "lon");
   SetUserOutputVariableName(2, "U");
   SetUserOutputVariableName(3, "V");
   SetUserOutputVariableName(4, "sqrtg");
+  SetUserOutputVariableName(5, "s1");
+  SetUserOutputVariableName(6, "s2");
 }
 
 void Mesh::InitUserMeshData(ParameterInput *pin) {

diff --git a/patches/24.time_integrator.patch b/patches/24.time_integrator.patch
@@ -1,5 +1,5 @@
 diff --git a/src/task_list/time_integrator.cpp b/src/task_list/time_integrator.cpp
-index 2bfc5178..e0fca29d 100644
+index 2bfc5178..173b2566 100644
 --- a/src/task_list/time_integrator.cpp
 +++ b/src/task_list/time_integrator.cpp
 @@ -36,6 +36,13 @@
@@ -16,13 +16,15 @@ index 2bfc5178..e0fca29d 100644
  //----------------------------------------------------------------------------------------
  //! TimeIntegratorTaskList constructor
 
-@@ -998,6 +1005,14 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
+@@ -998,6 +1005,16 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
        AddTask(SEND_HYD,src_aterm);
        AddTask(RECV_HYD,NONE);
        AddTask(SETB_HYD,(RECV_HYD|SRC_TERM));
 +
 +#ifdef CUBED_SPHERE
-+      AddTask(CLEAR_ALLBND2, SETB_HYD);
++      if (NSCALARS == 0) {
++        AddTask(CLEAR_ALLBND2, SETB_HYD);
++      }
 +      AddTask(RESTART_RECV, CLEAR_ALLBND2);
 +      AddTask(SEND_HYD2, CLEAR_ALLBND2);
 +      AddTask(RECV_HYD2, RESTART_RECV);
@@ -31,12 +33,25 @@ index 2bfc5178..e0fca29d 100644
      }
 
      if (SHEAR_PERIODIC) {
-@@ -1157,9 +1172,17 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
+@@ -1035,6 +1052,12 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
+         AddTask(SEND_SCLR,SRC_TERM);
+         AddTask(RECV_SCLR,NONE);
+         AddTask(SETB_SCLR,(RECV_SCLR|SRC_TERM));
++#ifdef CUBED_SPHERE
++        AddTask(CLEAR_ALLBND2, (SETB_HYD|SETB_SCLR));
++        AddTask(SEND_SCLR2, CLEAR_ALLBND2);
++        AddTask(RECV_SCLR2, RESTART_RECV);
++        AddTask(SETB_SCLR2, RECV_SCLR2);
++#endif
+       }
+       if (SHEAR_PERIODIC) {
+         AddTask(SEND_SCLRSH,SETB_SCLR);
+@@ -1157,9 +1180,17 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
            }
          } else {
            if (NSCALARS > 0) {
 +#ifdef CUBED_SPHERE
-+            AddTask(CONS2PRIM,(SETB_HYD2|SETB_SCLR));
++            AddTask(CONS2PRIM,(SETB_HYD2|SETB_SCLR2));
 +#else
              AddTask(CONS2PRIM,(SETB_HYD|SETB_SCLR));
 +#endif  // CUBED_SPHERE
@@ -49,7 +64,7 @@ index 2bfc5178..e0fca29d 100644
            }
          }
        }
-@@ -1225,11 +1248,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
+@@ -1225,11 +1256,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
    //!   SRC_TERM and SourceTerms(), USERWORK, PHY_BVAL, PROLONG, CONS2PRIM,
    //!   ... Although, AMR_FLAG = "flag blocks for AMR" should be FLAG_AMR in VERB_OBJECT
    using namespace HydroIntegratorTaskNames; // NOLINT (build/namespace)
@@ -67,7 +82,7 @@ index 2bfc5178..e0fca29d 100644
    } else if (id == CALC_HYDFLX) {
      task_list_[ntasks].TaskFunc=
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
-@@ -1275,7 +1303,7 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
+@@ -1275,7 +1311,7 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::AddSourceTerms);
      task_list_[ntasks].lb_time = true;
@@ -76,7 +91,7 @@ index 2bfc5178..e0fca29d 100644
      task_list_[ntasks].TaskFunc=
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::SendHydro);
-@@ -1285,17 +1313,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
+@@ -1285,17 +1321,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::SendField);
      task_list_[ntasks].lb_time = true;
@@ -96,7 +111,28 @@ index 2bfc5178..e0fca29d 100644
      task_list_[ntasks].TaskFunc=
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::SetBoundariesHydro);
-@@ -1688,6 +1715,19 @@ TaskStatus TimeIntegratorTaskList::ClearAllBoundary(MeshBlock *pmb, int stage) {
+@@ -1405,17 +1440,17 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
+         static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
+         (&TimeIntegratorTaskList::IntegrateScalars);
+     task_list_[ntasks].lb_time = true;
+-  } else if (id == SEND_SCLR) {
++  } else if (id == SEND_SCLR or id == SEND_SCLR2) {
+     task_list_[ntasks].TaskFunc=
+         static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
+         (&TimeIntegratorTaskList::SendScalars);
+     task_list_[ntasks].lb_time = true;
+-  } else if (id == RECV_SCLR) {
++  } else if (id == RECV_SCLR or id == RECV_SCLR2) {
+     task_list_[ntasks].TaskFunc=
+         static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
+         (&TimeIntegratorTaskList::ReceiveScalars);
+     task_list_[ntasks].lb_time = false;
+-  } else if (id == SETB_SCLR) {
++  } else if (id == SETB_SCLR or id == SETB_SCLR2) {
+     task_list_[ntasks].TaskFunc=
+         static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
+         (&TimeIntegratorTaskList::SetBoundariesScalars);
+@@ -1688,6 +1723,19 @@ TaskStatus TimeIntegratorTaskList::ClearAllBoundary(MeshBlock *pmb, int stage) {
    return TaskStatus::success;
  }
 
@@ -116,7 +152,7 @@ index 2bfc5178..e0fca29d 100644
  //----------------------------------------------------------------------------------------
  // Functions to calculates Hydro fluxes
 
-@@ -2305,6 +2345,12 @@ TaskStatus TimeIntegratorTaskList::PhysicalBoundary(MeshBlock *pmb, int stage) {
+@@ -2305,6 +2353,12 @@ TaskStatus TimeIntegratorTaskList::PhysicalBoundary(MeshBlock *pmb, int stage) {
  TaskStatus TimeIntegratorTaskList::UserWork(MeshBlock *pmb, int stage) {
    if (stage != nstages) return TaskStatus::success; // only do on last stage
 
@@ -129,7 +165,7 @@ index 2bfc5178..e0fca29d 100644
    pmb->UserWorkInLoop();
    return TaskStatus::success;
  }
-@@ -2395,6 +2441,7 @@ TaskStatus TimeIntegratorTaskList::IntegrateScalars(MeshBlock *pmb, int stage) {
+@@ -2395,6 +2449,7 @@ TaskStatus TimeIntegratorTaskList::IntegrateScalars(MeshBlock *pmb, int stage) {
        ave_wghts[2] = stage_wghts[stage-1].gamma_3;
        if (ave_wghts[0] == 0.0 && ave_wghts[1] == 1.0 && ave_wghts[2] == 0.0) {
          ps->s.SwapAthenaArray(ps->s1);

diff --git a/patches/31.task_list.patch b/patches/31.task_list.patch
@@ -1,5 +1,5 @@
 diff --git a/src/task_list/task_list.hpp b/src/task_list/task_list.hpp
-index d71ae3ea..ab14ee30 100644
+index d71ae3ea..931b17c2 100644
 --- a/src/task_list/task_list.hpp
 +++ b/src/task_list/task_list.hpp
 @@ -146,6 +146,7 @@ class TimeIntegratorTaskList : public TaskList {
@@ -10,7 +10,7 @@ index d71ae3ea..ab14ee30 100644
 
    TaskStatus CalculateHydroFlux(MeshBlock *pmb, int stage);
    TaskStatus CalculateEMF(MeshBlock *pmb, int stage);
-@@ -387,5 +388,11 @@ const TaskID RECV_RADSH(73);
+@@ -387,5 +388,15 @@ const TaskID RECV_RADSH(73);
 
  const TaskID SRCTERM_IMRAD(74);
 
@@ -19,6 +19,10 @@ index d71ae3ea..ab14ee30 100644
 +const TaskID SETB_HYD2(103);
 +const TaskID CLEAR_ALLBND2(104);
 +const TaskID RESTART_RECV(105);
++
++const TaskID SEND_SCLR2(106);
++const TaskID RECV_SCLR2(107);
++const TaskID SETB_SCLR2(108);
 +
  }  // namespace HydroIntegratorTaskNames
  #endif  // TASK_LIST_TASK_LIST_HPP_
diff --git a/patches/32.scalars.patch b/patches/32.scalars.patch
@@ -0,0 +1,12 @@
+diff --git a/src/scalars/scalars.cpp b/src/scalars/scalars.cpp
+index 8471fb20..2f3ad1ce 100644
+--- a/src/scalars/scalars.cpp
++++ b/src/scalars/scalars.cpp
+@@ -84,6 +84,7 @@ PassiveScalars::PassiveScalars(MeshBlock *pmb, ParameterInput *pin)  :
+
+   // enroll CellCenteredBoundaryVariable object
+   sbvar.bvar_index = pmb->pbval->bvars.size();
++  sbvar.TypeFlag = 1; // PassiveScalars
+   pmb->pbval->bvars.push_back(&sbvar);
+   pmb->pbval->bvars_main_int.push_back(&sbvar);
+   if (STS_ENABLED) {
diff --git a/patches/28.bvals_cc.patch → patches/33.bvals_cc_cpp.patch b/patches/28.bvals_cc.patch → patches/33.bvals_cc_cpp.patch
@@ -1,24 +1,24 @@
 diff --git a/src/bvals/cc/bvals_cc.cpp b/src/bvals/cc/bvals_cc.cpp
-index eb6f6039..d8b328b0 100644
+index 223962de..74b3ab79 100644
 --- a/src/bvals/cc/bvals_cc.cpp
 +++ b/src/bvals/cc/bvals_cc.cpp
-@@ -30,6 +30,9 @@
- #include "../../utils/buffer_utils.hpp"
+@@ -31,6 +31,9 @@
  #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,
+ #include <mpi.h>
+@@ -297,7 +300,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);
++  CubedSphereUtility::PackData(var, buf, nl_, nu_, si, ei, sj, ej, sk, ek, p, nb.ni.ox1, nb.ni.ox2, nb.ni.ox3, pmb->loc, TypeFlag);
 +#else
    BufferUtility::PackData(var, buf, nl_, nu_, si, ei, sj, ej, sk, ek, p);
 +#endif

diff --git a/patches/34.bvals_cc_hpp.patch b/patches/34.bvals_cc_hpp.patch
@@ -0,0 +1,14 @@
+diff --git a/src/bvals/cc/bvals_cc.hpp b/src/bvals/cc/bvals_cc.hpp
+index d167d729..ce2f8dcc 100644
+--- a/src/bvals/cc/bvals_cc.hpp
++++ b/src/bvals/cc/bvals_cc.hpp
+@@ -39,6 +39,9 @@ class CellCenteredBoundaryVariable : public BoundaryVariable {
+                                AthenaArray<Real> *var_flux, bool fflux, int flag);
+   ~CellCenteredBoundaryVariable();
+
++  // Type information
++  int TypeFlag = 2; // Hydro by default
++
+   //! \note
+   //! may want to rebind var_cc to u,u1,u2,w,w1, etc. registers for time integrator logic.
+   //! Also, derived class HydroBoundaryVariable needs to keep switching var and coarse_var