Add basic diagnostics (#126)

- add diagnostics that does not depend on communication

.github/workflows/main.yml

@@ -34,6 +34,54 @@ jobs:
       - name: Run cpplint
         run: ./tools/check_style.sh
 
+  build-cubed-sphere-2d:
+    if: github.event.pull_request.draft == false
+    runs-on: ubuntu-latest
+    needs: style-checker
+
+    steps:
+      - name: install dependencies
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y lcov libnetcdf-dev libpnetcdf-dev libeigen3-dev libmpich-dev
+
+      - uses: actions/checkout@v3
+        with:
+          lfs: false
+
+      - name: create build directory
+        # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
+        # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
+        run: cmake -B ${{github.workspace}}/build -DTASK=exo2 -DCMAKE_CXX_FLAGS=${{ env.BUILD_FLAG }} -DCMAKE_C_FLAGS=${{ env.BUILD_FLAG }}
+
+      - name: Build
+        # Build your program with the given configuration
+        run: cmake --build ${{github.workspace}}/build --config Debug -- -j2
+
+  build-cubed-sphere-3d:
+    if: github.event.pull_request.draft == false
+    runs-on: ubuntu-latest
+    needs: style-checker
+
+    steps:
+      - name: install dependencies
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y lcov libnetcdf-dev libpnetcdf-dev libeigen3-dev libmpich-dev
+
+      - uses: actions/checkout@v3
+        with:
+          lfs: false
+
+      - name: create build directory
+        # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make.
+        # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type
+        run: cmake -B ${{github.workspace}}/build -DTASK=exo3 -DCMAKE_CXX_FLAGS=${{ env.BUILD_FLAG }} -DCMAKE_C_FLAGS=${{ env.BUILD_FLAG }}
+
+      - name: Build
+        # Build your program with the given configuration
+        run: cmake --build ${{github.workspace}}/build --config Debug -- -j2
+
   build-non-hydrostatic:
     # The CMake configure and build commands are platform agnostic and should work equally well on Windows or Mac.
     # You can convert this to a matrix build if you need cross-platform coverage.

examples/2019-Li-snap/straka.inp

@@ -21,6 +21,11 @@ file_type  = netcdf
 variable   = uov
 dt         = 300.
 
+<output4>
+file_type  = netcdf
+variable   = diag
+dt         = 300.
+
 <time>
 cfl_number = 0.9          # The Courant, Friedrichs, & Lewy (CFL) Number
 nlim       = -1           # cycle limit
@@ -64,3 +69,5 @@ zr    = 2.E3
 p0    = 1.E5
 Ts    = 300.
 K     = 75.
+
+diagnostics = div,curl

patches/24.time_integrator.patch

@@ -1,19 +1,22 @@
 diff --git a/src/task_list/time_integrator.cpp b/src/task_list/time_integrator.cpp
-index 097afa53..225ab61f 100644
+index 2bfc5178..e0fca29d 100644
 --- a/src/task_list/time_integrator.cpp
 +++ b/src/task_list/time_integrator.cpp
-@@ -36,6 +36,10 @@
+@@ -36,6 +36,13 @@
  #include "../scalars/scalars.hpp"
  #include "task_list.hpp"
 
 +// canoe injection
 +#include <configure.hpp>
 +#include <impl.hpp>
++
++// diagnostics
++#include <diagnostics/diagnostics.hpp>
 +
  //----------------------------------------------------------------------------------------
  //! TimeIntegratorTaskList constructor
 
-@@ -998,6 +1002,14 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
+@@ -998,6 +1005,14 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
        AddTask(SEND_HYD,src_aterm);
        AddTask(RECV_HYD,NONE);
        AddTask(SETB_HYD,(RECV_HYD|SRC_TERM));
@@ -28,7 +31,7 @@ index 097afa53..225ab61f 100644
      }
 
      if (SHEAR_PERIODIC) {
-@@ -1156,9 +1168,17 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
+@@ -1157,9 +1172,17 @@ TimeIntegratorTaskList::TimeIntegratorTaskList(ParameterInput *pin, Mesh *pm) {
            }
          } else {
            if (NSCALARS > 0) {
@@ -46,7 +49,7 @@ index 097afa53..225ab61f 100644
            }
          }
        }
-@@ -1224,11 +1244,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
+@@ -1225,11 +1248,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)
@@ -64,7 +67,7 @@ index 097afa53..225ab61f 100644
    } else if (id == CALC_HYDFLX) {
      task_list_[ntasks].TaskFunc=
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
-@@ -1274,7 +1299,7 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
+@@ -1275,7 +1303,7 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::AddSourceTerms);
      task_list_[ntasks].lb_time = true;
@@ -73,7 +76,7 @@ index 097afa53..225ab61f 100644
      task_list_[ntasks].TaskFunc=
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::SendHydro);
-@@ -1284,17 +1309,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
+@@ -1285,17 +1313,16 @@ void TimeIntegratorTaskList::AddTask(const TaskID& id, const TaskID& dep) {
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::SendField);
      task_list_[ntasks].lb_time = true;
@@ -93,7 +96,7 @@ index 097afa53..225ab61f 100644
      task_list_[ntasks].TaskFunc=
          static_cast<TaskStatus (TaskList::*)(MeshBlock*,int)>
          (&TimeIntegratorTaskList::SetBoundariesHydro);
-@@ -1683,6 +1707,19 @@ TaskStatus TimeIntegratorTaskList::ClearAllBoundary(MeshBlock *pmb, int stage) {
+@@ -1688,6 +1715,19 @@ TaskStatus TimeIntegratorTaskList::ClearAllBoundary(MeshBlock *pmb, int stage) {
    return TaskStatus::success;
  }
 
@@ -113,7 +116,20 @@ index 097afa53..225ab61f 100644
  //----------------------------------------------------------------------------------------
  // Functions to calculates Hydro fluxes
 
-@@ -2388,6 +2425,7 @@ TaskStatus TimeIntegratorTaskList::IntegrateScalars(MeshBlock *pmb, int stage) {
+@@ -2305,6 +2345,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
+
++  auto all_diags = pmb->pimpl->all_diags;
++
++  for (auto& diag : all_diags) {
++    diag->Progress(pmb->phydro->w);
++  }
++
+   pmb->UserWorkInLoop();
+   return TaskStatus::success;
+ }
+@@ -2395,6 +2441,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);

src/CMakeLists.txt

@@ -53,6 +53,7 @@ add_subdirectory(exchanger)
 add_subdirectory(nbody)
 add_subdirectory(special)
 add_subdirectory(single_column)
+add_subdirectory(diagnostics)
 
 set(CANOE_LIBRARY_${buildu}
   "canoe_${buildl}"
@@ -74,8 +75,8 @@ set(CANOE_LIBRARY_${buildu}
   "exchanger_${buildl}"
   "nbody_${buildl}"
   "scm_${buildl}"
+  "diagnostics_${buildl}"
   "yaml-cpp"
-  #${CDISORT_LIBRARY_${buildu}}
   ${CPPDISORT_LIBRARY_${buildu}}
   ${PYTHON_LIBRARY_RELEASE}
   ${NETCDF_LIBRARIES}

src/astro/celestrial_body.cpp

@@ -8,18 +8,19 @@
 
 // athena
 #include <athena/athena.hpp>
+#include <athena/mesh/mesh.hpp>
 #include <athena/parameter_input.hpp>
 
 // climath
 #include <climath/core.h>
 #include <climath/interpolation.h>
 
+// harp
+#include <harp/spectral_grid.hpp>
+
 // utils
 #include <utils/fileio.hpp>
 
-// canoe
-#include <common.hpp>
-
 // astro
 #include "celestrial_body.hpp"
 
@@ -161,7 +162,8 @@ Real CelestrialBody::ParentDistanceInAu(Real time) const {
                                             equinox - M_PI / 2.)));
 }
 
-CelestrialBodyPtr PlanetFactory::Create(ParameterInput *pin) {
+CelestrialBodyPtr PlanetFactory::CreateFrom(MeshBlock *pmb,
+                                            ParameterInput *pin) {
   if (pin->DoesParameterExist("astronomy", "planet")) {
     std::string name = pin->GetString("astronomy", "planet");
     return std::make_shared<CelestrialBody>(pin, name);

src/astro/celestrial_body.hpp

@@ -14,6 +14,7 @@
 class ParameterInput;
 struct float_triplet;
 struct Direction;
+class Meshblock;
 
 class CelestrialBody : public NamedGroup {
  public:
@@ -67,7 +68,7 @@ Real GetGravity(char const *name, Real pclat);
 
 class PlanetFactory {
  public:
-  static CelestrialBodyPtr Create(ParameterInput *pin);
+  static CelestrialBodyPtr CreateFrom(MeshBlock *pmb, ParameterInput *pin);
 };
 
 #endif  // SRC_ASTRO_CELESTRIAL_BODY_HPP_

src/diagnostics/CMakeLists.txt

@@ -0,0 +1,19 @@
+set(namel diagnostics)
+string(TOUPPER ${namel} nameu)
+
+file(GLOB src_files
+    *.cpp
+    )
+
+string(TOLOWER ${CMAKE_BUILD_TYPE} buildl)
+string(TOUPPER ${CMAKE_BUILD_TYPE} buildu)
+
+add_library(${namel}_${buildl}
+    OBJECT
+    ${src_files}
+    )
+
+set_target_properties(${namel}_${buildl}
+    PROPERTIES
+    COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}}
+    )

src/snap/diagnostics/buoyancy.cpp → src/diagnostics/buoyancy.cpp

@@ -1,13 +1,14 @@
-#include "../coordinates/coordinates.hpp"
-#include "../hydro/hydro.hpp"
-#include "../hydro/srcterms/hydro_srcterms.hpp"
-#include "../reconstruct/interpolation.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/hydro/hydro.hpp>
+#include <athena/hydro/srcterms/hydro_srcterms.hpp>
+#include <athena/reconstruct/interpolation.hpp>
+
+// diagnostics
 #include "diagnostics.hpp"
 
 Buoyancy::Buoyancy(MeshBlock* pmb) : Diagnostics(pmb, "b") {
   type = "SCALARS";
-  long_name = "buoyancy";
-  units = "m/s^2";
   data.NewAthenaArray(ncells3_, ncells2_, ncells1_);
   pf_.NewAthenaArray(ncells3_, ncells2_, ncells1_ + 1);
   grav_ = pmb->phydro->hsrc.GetG1();

src/snap/diagnostics/curl.cpp → src/diagnostics/curl.cpp

@@ -1,10 +1,11 @@
-#include "../coordinates/coordinates.hpp"
-#include "../reconstruct/interpolation.hpp"
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/reconstruct/interpolation.hpp>
+
+// canoe
 #include "diagnostics.hpp"
 
 Curl::Curl(MeshBlock *pmb) : Diagnostics(pmb, "curl") {
-  long_name = "curl";
-  units = "1/s";
   if (pmb->block_size.nx3 > 1) {
     type = "VECTORS";
     data.NewAthenaArray(3, ncells3_, ncells2_, ncells1_);

src/diagnostics/diagnostics.cpp

@@ -0,0 +1,46 @@
+// C/C++ headers
+#include <cstring>
+#include <sstream>
+#include <stdexcept>
+
+// application
+#include <application/application.hpp>
+
+// Athena++ headers
+#include <athena/mesh/mesh.hpp>
+
+// diagnostics
+#include "diagnostics.hpp"
+
+Diagnostics::Diagnostics(MeshBlock *pmb, std::string name)
+    : NamedGroup(name), pmy_block_(pmb), ncycle_(0) {
+  Application::Logger app("main");
+  app->Log("Initialize Diagnostics");
+
+  ncells1_ = pmb->block_size.nx1 + 2 * (NGHOST);
+  ncells2_ = 1;
+  ncells3_ = 1;
+  if (pmb->pmy_mesh->f2) ncells2_ = pmb->block_size.nx2 + 2 * (NGHOST);
+  if (pmb->pmy_mesh->f3) ncells3_ = pmb->block_size.nx3 + 2 * (NGHOST);
+
+  x1edge_.NewAthenaArray(ncells1_ + 1);
+  x1edge_p1_.NewAthenaArray(ncells1_);
+  x2edge_.NewAthenaArray(ncells1_ + 1);
+  x2edge_p1_.NewAthenaArray(ncells1_);
+  x3edge_.NewAthenaArray(ncells1_ + 1);
+  x3edge_p1_.NewAthenaArray(ncells1_);
+
+  x1area_.NewAthenaArray(ncells1_ + 1);
+  x2area_.NewAthenaArray(ncells1_);
+  x2area_p1_.NewAthenaArray(ncells1_);
+  x3area_.NewAthenaArray(ncells1_);
+  x3area_p1_.NewAthenaArray(ncells1_);
+
+  vol_.NewAthenaArray(ncells1_);
+  total_vol_.NewAthenaArray(ncells1_);
+}
+
+Diagnostics::~Diagnostics() {
+  Application::Logger app("main");
+  app->Log("Destroy Diagnostics");
+}