Add polar vortex intruder test (#77)

.github/workflows/mac.yml

@@ -153,9 +153,6 @@ jobs:
         with:
           lfs: false
 
-      - name: set up git lfs
-        run: git lfs install
-
       - name: fetch the needed lfs file
         run: git lfs pull -I examples/2020-Li-polar-vortex/swxy_output.nc
 
@@ -172,7 +169,7 @@ jobs:
 
       - name: Build
         # Build your program with the given configuration
-        run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}} -- -j4
+        run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}} -- -j3
 
       - name: run simulation
         run: cd ${{github.workspace}}/build/bin && mpiexec -n 4 ./shallow_xy.release -i shallow_water.inp
@@ -182,6 +179,44 @@ jobs:
 
       - name: compare result
         run: >
-          python3 ${{github.workspace}}/examples/2020-Li-polar-vortex/test_shallow.py
+          python3 ${{github.workspace}}/examples/2020-Li-polar-vortex/test_shallow_xy.py
           ${{github.workspace}}/build/bin/sw-test-main.nc
           ${{github.workspace}}/examples/2020-Li-polar-vortex/swxy_output.nc
+
+  polar_vortex:
+    if: github.event.pull_request.draft == false
+    runs-on: macos-13
+    steps:
+      - uses: actions/checkout@v3
+        with:
+          lfs: false
+
+      - name: fetch the needed lfs file
+        run: git lfs pull -I examples/2020-Li-polar-vortex/intruder_output.nc
+
+      - name: set up libraries
+        run: brew bundle
+
+      - name: set up python libraries
+        run: pip3 install -r requirements.txt
+
+      - 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 -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} -DCMAKE_Fortran_COMPILER=gfortran-12 -DTASK=vortex
+
+      - name: Build
+        # Build your program with the given configuration
+        run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}} -- -j3
+
+      - name: run simulation
+        run: cd ${{github.workspace}}/build/bin && mpiexec -n 3 ./polar_vortex.release -i fig_g_intruder.inp
+
+      - name: combile outputs
+        run: cd ${{github.workspace}}/build/bin && ./combine.py -o test
+
+      - name: compare result
+        run: >
+          python3 ${{github.workspace}}/examples/2020-Li-polar-vortex/test_polar_vortex.py
+          ${{github.workspace}}/build/bin/intruder-test-main.nc
+          ${{github.workspace}}/examples/2020-Li-polar-vortex/intruder_output.nc

.github/workflows/main.yml

@@ -380,6 +380,39 @@ jobs:
 
       - name: compare result
         run: >
-          python3 ${{github.workspace}}/examples/2020-Li-polar-vortex/test_shallow.py
+          python3 ${{github.workspace}}/examples/2020-Li-polar-vortex/test_shallow_xy.py
           ${{github.workspace}}/build/bin/sw-test-main.nc
           ${{github.workspace}}/examples/2020-Li-polar-vortex/swxy_output.nc
+
+  polar_vortex:
+    if: github.event.pull_request.draft == false
+    runs-on: [self-hosted, Linux, 16core]
+    needs: build-non-hydrostatic
+    steps:
+      - uses: actions/checkout@v3
+        with:
+          lfs: false
+
+      - name: fetch the needed lfs file
+        run: git lfs pull -I examples/2020-Li-polar-vortex/intruder_output.nc
+
+      - 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 -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} -DTASK=vortex
+
+      - name: Build
+        # Build your program with the given configuration
+        run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}} -- -j8
+
+      - name: run simulation
+        run: cd ${{github.workspace}}/build/bin && mpiexec -n 12 ./polar_vortex.release -i fig_g_intruder.inp
+
+      - name: combile outputs
+        run: cd ${{github.workspace}}/build/bin && ./combine.py -o test
+
+      - name: compare result
+        run: >
+          python3 ${{github.workspace}}/examples/2020-Li-polar-vortex/test_polar_vortex.py
+          ${{github.workspace}}/build/bin/intruder-test-main.nc
+          ${{github.workspace}}/examples/2020-Li-polar-vortex/intruder_output.nc

Brewfile

@@ -16,3 +16,4 @@ brew "wget"
 brew "node"
 brew "cfitsio"
 brew "gh"
+brew "boost"

README.md

@@ -24,13 +24,13 @@ brew bundle
 ### Ubuntu Linux Installation Guide
 On a Ubuntu linux system, use `apt` to install
 ```
-sudo apt install clang-format cmake nco libnetcdf-dev libpnetcdf-dev
+sudo apt install clang-format cmake nco libnetcdf-dev libpnetcdf-dev libboost-all-dev
 ```
 
 ### Redhat Linux Installation Guide
 On a Redhat linux system, use `yum` to install
 ```
-sudo yum install clang-tools-extra cmake nco netcdf-devel
+sudo yum install clang-tools-extra cmake nco netcdf-devel boost boost-devel
 ```
 
 ## Install python libraries

cmake/vortex.cmake

@@ -0,0 +1,14 @@
+# configuration for polar vortex simulation
+
+macro(SET_IF_EMPTY _variable)
+  if("${${_variable}}" STREQUAL "")
+    set(${_variable} ${ARGN})
+  endif()
+endmacro()
+
+# athena variables
+set_if_empty(NUMBER_GHOST_CELLS 3)
+
+set(EOS shallow_xy)
+set(PNETCDF ON)
+set(MPI ON)

doc/README-mods.md

@@ -0,0 +1,5 @@
+# Custom configurations
+
+- `PNETCDF (OFF)` : turn on pnetcdf output
+- `NETCDF (OFF)` : turn on netcdf output
+- `USER_MODS` : user module file

examples/2020-Li-polar-vortex/.gitattributes

@@ -1 +1,2 @@
 swxy_output.nc filter=lfs diff=lfs merge=lfs -text
+intruder_output.nc filter=lfs diff=lfs merge=lfs -text

examples/2020-Li-polar-vortex/CMakeLists.txt

@@ -5,7 +5,9 @@ string(TOUPPER ${CMAKE_BUILD_TYPE} buildu)
 
 #set(namel straka)
 
-foreach(namel shallow_xy)
+foreach(namel shallow_xy polar_vortex)
+  find_package(Boost REQUIRED)
+
   add_executable(${namel}.${buildl} ${namel}.cpp
     ${CMAKE_BINARY_DIR}/main_${TASK}.cpp)
 

examples/2020-Li-polar-vortex/fig_g_intruder.inp

@@ -4,7 +4,7 @@ reference =
 configure = --prob=polar_vortex --eos=shallow_water --nghost=3 -mpi -pnetcdf --pnetcdf_path=${HOME}/opt --include=src/tools/
 
 <job>
-problem_id = data/intruder # problem ID: basename of output filenames
+problem_id = intruder # problem ID: basename of output filenames
 
 <output0>
 file_type  = hst          # History data dump
@@ -32,13 +32,13 @@ xorder     = 5            # horizontal reconstruction order
 integrator = rk3          # integration method
 
 <mesh>
-nx1        = 256          # Number of zones in X1-direction (longitude)
+nx1        = 240          # Number of zones in X1-direction (longitude)
 x1min      = -3.5E7       # minimum value of X1
 x1max      = 3.5E7        # maximum value of X1, 120 degree
 ix1_bc     = periodic     # inner-X1 boundary flag
 ox1_bc     = periodic     # inner-X1 boundary flag
 
-nx2        = 256          # Number of zones in X2-direction (latitude)
+nx2        = 240          # Number of zones in X2-direction (latitude)
 x2min      = -3.5E7       # minimum value of X2
 x2max      = 3.5E7        # maximum value of X2, 70 degree
 ix2_bc     = periodic     # inner-X2 boundary flag
@@ -54,8 +54,8 @@ ox3_bc     = periodic     # inner-X3 boundary flag
 gamma      = 1.4
 
 <meshblock>
-nx1       = 64
-nx2       = 64
+nx1       = 60
+nx2       = 80
 nx3       = 1
 
 <problem>

examples/2020-Li-polar-vortex/polar_vortex.cpp

@@ -1,28 +1,34 @@
 // changes: no time loop, exp v2, vphi = vmax (positive #), coriolis not
 // approximated, damping model--Model #2--ingersoll's equation n=2
-//! \file vortex.cpp
+//! \file polar_vortex.cpp
 //  \brief jupiter polar vortex model
 
-// C++ headers
-#include <boost/math/special_functions/gamma.hpp>
+// C/C++
 #include <cmath>
 #include <random>
 #include <sstream>
 #include <stdexcept>
 #include <vector>
 
-// Athena++ headers
-#include "../athena.hpp"
-#include "../athena_arrays.hpp"
-#include "../coordinates/coordinates.hpp"
-#include "../eos/eos.hpp"
-#include "../field/field.hpp"
-#include "../globals.hpp"
-#include "../hydro/hydro.hpp"
-#include "../math/core.h"  // sqr
-#include "../mesh/mesh.hpp"
-#include "../parameter_input.hpp"
-#include "../utils/utils.hpp"  // ReadDataTable
+// boost
+#include <boost/math/special_functions/gamma.hpp>
+
+// Athena++
+#include <athena/athena.hpp>
+#include <athena/athena_arrays.hpp>
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/eos/eos.hpp>
+#include <athena/field/field.hpp>
+#include <athena/globals.hpp>
+#include <athena/hydro/hydro.hpp>
+#include <athena/mesh/mesh.hpp>
+#include <athena/parameter_input.hpp>
+
+// climath
+#include <climath/core.h>  // sqr
+
+// utils
+#include <utils/fileio.hpp>  // ReadDataTable
 
 // global parameters
 Real phi0, phi1, lambda, alpha, vphi, vrad, vis, polarity, skewness, omega,
@@ -55,8 +61,9 @@ Real Phib(Coordinates *pcoord, int j,
 }
 
 void PolarVortexForcing(MeshBlock *pmb, const Real time, const Real dt,
-                        const AthenaArray<Real> &w,
-                        AthenaArray<Real> const &bcc, AthenaArray<Real> &u) {
+                        AthenaArray<Real> const &w, AthenaArray<Real> const &r,
+                        AthenaArray<Real> const &bcc, AthenaArray<Real> &u,
+                        AthenaArray<Real> &s) {
   for (int j = pmb->js; j <= pmb->je; ++j)
     for (int i = pmb->is; i <= pmb->ie; ++i) {
       Real x1 = pmb->pcoord->x1v(i);
@@ -219,7 +226,7 @@ void MeshBlock::ProblemGenerator(
     ParameterInput *pin)  // sets initial conditions
 {
   AthenaArray<Real> vpos;
-  ReadDataTable(vpos, pin->GetString("problem", "vfile"));
+  ReadDataTable(&vpos, pin->GetString("problem", "vfile"));
   int vnum = vpos.GetDim2();  // input file allows positions
   // int vnum = 6;
   for (int i = 0; i < vnum; ++i) {

examples/2020-Li-polar-vortex/shallow_xy.cpp

@@ -1,4 +1,4 @@
-//! \file shallow_water.cpp
+//! \file shallow_xy.cpp
 //  \brief shallow water test model
 
 // C/C++

examples/2020-Li-polar-vortex/test_polar_vortex.py

@@ -0,0 +1,20 @@
+#! /usr/bin/env python3
+from netCDF4 import Dataset
+from numpy import sum, sqrt
+import sys
+
+case1 = sys.argv[1]
+case2 = sys.argv[2]
+
+data1 = Dataset(case1, "r")
+var1 = data1["vort"][-1, :, :, :]
+
+data2 = Dataset(case2, "r")
+var2 = data2["vort"][-1, :, :, :]
+
+diff = sqrt(sum((var2 - var1) * (var2 - var1)))
+
+if diff < 1.0e-6:
+    print("### polar vortex case passed. ###")
+else:
+    raise ValueError("ERROR: polar vortex failed. L2-norm is %.2g" % diff)

examples/2024-Li-plume/plume_pgen.cpp

@@ -67,28 +67,10 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
       }
 }
 
-void CondensateGravity(MeshBlock *pmb, Real const time, Real const dt,
-                       AthenaArray<Real> const &w, AthenaArray<Real> const &r,
-                       AthenaArray<Real> const &bcc, AthenaArray<Real> &u,
-                       AthenaArray<Real> &s) {
-  // acceleration in 1-direction
-  for (int k = pmb->ks; k <= pmb->ke; ++k)
-    for (int j = pmb->js; j <= pmb->je; ++j)
-      for (int i = pmb->is; i <= pmb->ie; ++i) {
-        Real qd = 1.;
-#pragma omp simd reduction(+ : qd)
-        for (int n = 1; n <= NVAPOR; ++n) qd += -w(n, k, j, i);
-
-        Real rho_dry = w(IDN, k, j, i) * qd;
-        Real rho_cloud = 0.;
-#pragma omp simd
-        for (int n = 0; n < NCLOUD; ++n) rho_cloud += rho_dry * r(n, k, j, i);
-
-        Real src = -dt * rho_cloud * grav;
-        u(IM1, k, j, i) += src;
-        if (NON_BAROTROPIC_EOS) u(IEN, k, j, i) += src * w(IVX, k, j, i);
-      }
-
+void SurfacePlumeSource(MeshBlock *pmb, Real const time, Real const dt,
+                        AthenaArray<Real> const &w, AthenaArray<Real> const &r,
+                        AthenaArray<Real> const &bcc, AthenaArray<Real> &u,
+                        AthenaArray<Real> &s) {
   // surface plume source
   Real x2mid =
       (pmb->pmy_mesh->mesh_size.x2min + pmb->pmy_mesh->mesh_size.x2max) / 2.;
@@ -129,7 +111,7 @@ void Mesh::InitUserMeshData(ParameterInput *pin) {
   flux_heat = pin->GetReal("problem", "flux_heat");
   flux_width = pin->GetReal("problem", "flux_width");
 
-  EnrollUserExplicitSourceFunction(CondensateGravity);
+  EnrollUserExplicitSourceFunction(SurfacePlumeSource);
 }
 
 void MeshBlock::ProblemGenerator(ParameterInput *pin) {