Fix exo3 in new canoe

CMakeLists.txt

@@ -10,6 +10,25 @@ project(canoe
 # search for <Package>_ROOT
 cmake_policy(SET CMP0074 NEW)
 
+find_program(Python3_EXECUTABLE NAMES python)
+find_package(Python3 REQUIRED COMPONENTS Interpreter Development.Module)
+execute_process(
+  COMMAND ${Python3_EXECUTABLE} -c "import torch; print(int(torch._C._GLIBCXX_USE_CXX11_ABI))"
+  RESULT_VARIABLE _ABI_RESULT
+  OUTPUT_VARIABLE _ABI_OUTPUT
+  OUTPUT_STRIP_TRAILING_WHITESPACE
+)
+
+if(NOT _ABI_RESULT EQUAL 0)
+  message(FATAL_ERROR "Failed to detect Torch GLIBCXX ABI")
+endif()
+
+# Set the _GLIBCXX_USE_CXX11_ABI flag
+add_compile_definitions(_GLIBCXX_USE_CXX11_ABI=${_ABI_OUTPUT})
+
+# Set MacOS deployment target
+set(CMAKE_OSX_DEPLOYMENT_TARGET "15.0")
+
 string(ASCII 27 Esc)
 set(ColorReset "${Esc}[m")
 set(Green      "${Esc}[32m")

cmake/compilers.cmake

@@ -16,12 +16,6 @@ if(CMAKE_CXX_COMPILER_ID MATCHES "GNU")
   set(CMAKE_CXX_FLAGS_DEBUG "-g3")
   set(CMAKE_C_FLAGS_DEBUG "-g3")
 
-  set(CMAKE_CXX_FLAGS_RELEASE
-      "${CMAKE_CXX_FLAGS_RELEASE} -D_GLIBCXX_USE_CXX11_ABI=0")
-  set(CMAKE_CXX_FLAGS_DEBUG
-      "${CMAKE_CXX_FLAGS_DEBUG} -D_GLIBCXX_USE_CXX11_ABI=0")
-
-  # set(CMAKE_Fortran_FLAGS_RELEASE "-O3" )
   set(KNOWN_COMPILER TRUE)
 endif()
 
@@ -32,7 +26,6 @@ if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
   set(CMAKE_CXX_FLAGS_DEBUG "-g3")
   set(CMAKE_C_FLAGS_DEBUG "-g3")
 
-  # set(CMAKE_Fortran_FLAGS_RELEASE "-O3" )
   set(KNOWN_COMPILER TRUE)
 endif()
 

cmake/examples/exo3.cmake

@@ -16,4 +16,4 @@ set(MPI ON)
 set(PNETCDF ON)
 set(RSOLVER hllc_transform)
 
-set (NTRACER 3)
+# set (NTRACER 3)

examples/2019-Li-snap/CMakeLists.txt

@@ -16,5 +16,5 @@ endif()
 # 4. Copy input files to run directory
 file(GLOB inputs *.inp *.yaml)
 foreach(input ${inputs})
-  execute_process(COMMAND ln -sf ${input} ${CMAKE_BINARY_DIR}/bin/${inp})
+  execute_process(COMMAND ln -sf ${input} ${CMAKE_BINARY_DIR}/bin/)
 endforeach()

examples/2023-Chen-exo3/CMakeLists.txt

@@ -19,7 +19,7 @@ else()
 endif()
 
 # 2. Copy input file to run directory
-file(GLOB inputs *.inp *.dat)
+file(GLOB inputs *.inp *.dat *.yaml)
 foreach(input ${inputs})
-    file(COPY ${input} DESTINATION ${CMAKE_BINARY_DIR}/bin)
+  execute_process(COMMAND ln -sf ${input} ${CMAKE_BINARY_DIR}/bin/)
 endforeach()

examples/2023-Chen-exo3/earth-dry.yaml

@@ -0,0 +1,14 @@
+
+phases:
+  - name: atm
+    thermo: ideal-moist
+    species: [atm]
+    kinetics: condensation
+    reactions: none
+
+species:
+  - name: atm
+    composition: {O: 0.42, N: 1.56, Ar: 0.01}
+    thermo:
+      model: constant-cp
+      cp0: 29.1 J/mol/K

examples/2023-Chen-exo3/hot-jupiter.yaml

@@ -0,0 +1,14 @@
+
+phases:
+  - name: atm
+    thermo: ideal-moist
+    species: [atm]
+    kinetics: condensation
+    reactions: none
+
+species:
+  - name: atm
+    composition: {H: 1.5, He: 0.15}
+    thermo:
+      model: constant-cp
+      cp0: 29.1 J/mol/K

examples/2023-Chen-exo3/hot_jupiter.cpp

@@ -29,7 +29,6 @@
 #include <application/exceptions.hpp>
 
 // canoe
-#include <air_parcel.hpp>
 #include <configure.hpp>
 #include <impl.hpp>
 
@@ -234,20 +233,23 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
 
   // construct an isothermal atmosphere
   auto pthermo = Thermodynamics::GetInstance();
-  AirParcel air(AirParcel::Type::MoleFrac);
+  auto &w = phydro->w;
+  std::vector<Real> yfrac(1, 1.);
 
   for (int k = ks; k <= ke; ++k)
     for (int j = js; j <= je; ++j) {
-      air.w[IPR] = p0;
-      air.w[IDN] = Ts;
+      pthermo->SetMassFractions<Real>(yfrac.data());
+      pthermo->EquilibrateTP(Ts, p0);
 
       int i = is;
       for (; i <= ie; ++i) {
-        AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
-        pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav, 0.001);
+        pthermo->GetPrimitive(w.at(k, j, i));
+        pthermo->Extrapolate_inplace(pcoord->dx1f(i), "isothermal", grav,
+                                     0.001);
+
         // add noise
-        air.w[IVY] = 10. * distribution(generator);
-        air.w[IVZ] = 10. * distribution(generator);
+        w(IVY, k, j, i) = 10. * distribution(generator);
+        w(IVZ, k, j, i) = 10. * distribution(generator);
       }
     }
 
@@ -277,8 +279,8 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
 
   // transfer to conservative variables
   // bcc is cell-centered magnetic fields, it is only a place holder here
-  // peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is,
-  // ie, js, je, ks, ke);
+  peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie,
+                             js, je, ks, ke);
 }
 
 void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {

examples/2023-Chen-exo3/hot_jupiter.inp

@@ -71,6 +71,8 @@ z_stra      = 2.E6        # Stratosphere height, m
 Gamma_trop  = 2.E-4       # Lapse rate, K/m
 
 <problem>
+thermodynamics_config = hot-jupiter.yaml
+
 Kt         = 1.5E5            # Damping timescale, s
 Ts         = 1600.            # Surface temperature, K
 p0         = 1.E5             # Surface pressure, Pa

examples/2023-Chen-exo3/hs94.cpp

@@ -29,7 +29,6 @@
 #include <application/exceptions.hpp>
 
 // canoe
-#include <air_parcel.hpp>
 #include <configure.hpp>
 #include <impl.hpp>
 
@@ -53,8 +52,6 @@ Real piso = 1E4;
 std::default_random_engine generator;
 std::normal_distribution<double> distribution(0.0, 1.0);
 
-namespace cs = CubedSphereUtility;
-
 // \brief Held-Suarez atmosphere benchmark test. Refernce: Held & Suarez,
 // (1994). Forcing parameters are given in the paper.
 
@@ -265,30 +262,38 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
 
   // construct an adiabatic atmosphere
   auto pthermo = Thermodynamics::GetInstance();
-  AirParcel air(AirParcel::Type::MoleFrac);
+  auto &w = phydro->w;
+  std::vector<Real> yfrac(1, 1.);
 
   for (int k = ks; k <= ke; ++k)
     for (int j = js; j <= je; ++j) {
-      air.w[IPR] = p0;
-      air.w[IDN] = Ts;
+      pthermo->SetMassFractions<Real>(yfrac.data());
+      pthermo->EquilibrateTP(Ts, p0);
+
+      // half a grid to cell center
+      pthermo->Extrapolate_inplace(pcoord->dx1f(is) / 2., "dry", grav);
 
       int i = is;
       for (; i <= ie; ++i) {
         if (pcoord->x1v(i) - Rp > z_iso) break;
 
-        AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
-        pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav, 0.001);
+        pthermo->GetPrimitive(w.at(k, j, i));
+        pthermo->Extrapolate_inplace(pcoord->dx1f(i), "dry", grav);
+
         // add noise
-        air.w[IVY] = 10. * distribution(generator);
-        air.w[IVZ] = 10. * distribution(generator);
+        w(IVY, k, j, i) = 10. * distribution(generator);
+        w(IVZ, k, j, i) = 10. * distribution(generator);
       }
 
       // construct isothermal atmosphere
       for (; i <= ie; ++i) {
-        AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
-        pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav);
+        pthermo->GetPrimitive(w.at(k, j, i));
+        pthermo->Extrapolate_inplace(pcoord->dx1f(i), "isothermal", grav);
       }
     }
+
+  peos->PrimitiveToConserved(w, pfield->bcc, phydro->u, pcoord, is, ie, js, je,
+                             ks, ke);
 }
 
 void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {

examples/2023-Chen-exo3/hs94.inp

@@ -64,14 +64,15 @@ gamma          = 1.4         # gamma = C_p/C_v
 grav_acc1      = -9.81         # gravity accelaration
 implicit_flag  = 1
 
-
 <thermodynamics>
 dT         = 60.       # horizontal temperature difference
 dtheta     = 10.       # vertical potential temperature difference
 Rd         = 287.
 cp         = 1004.
 
 <problem>
+thermodynamics_config = earth-dry.yaml
+
 Omega      = 7.292E-5
 Rp         = 6.371E6
 p0         = 1.E5      # surface pressure

examples/2023-Chen-exo3/polar_dry.cpp

@@ -210,14 +210,16 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   z_iso = pin->GetReal("problem", "z_iso");
   sponge_lat = pin->GetReal("problem", "sponge_lat");
   heat_flux = pin->GetReal("problem", "heat_flux");
+
   // construct an adiabatic atmosphere
   auto pthermo = Thermodynamics::GetInstance();
-  AirParcel air(AirParcel::Type::MoleFrac);
+  auto &w = phydro->w;
+  std::vector<Real> yfrac(1, 1.);
 
   for (int k = ks; k <= ke; ++k)
     for (int j = js; j <= je; ++j) {
-      air.w[IPR] = p0;
-      air.w[IDN] = Ts;
+      pthermo->SetMassFractions<Real>(yfrac.data());
+      pthermo->EquilibrateTP(Ts, p0);
 
       int i = is;
       // for (; i <= ie; ++i) {
@@ -232,11 +234,12 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
 
       // construct isothermal atmosphere
       for (; i <= ie; ++i) {
-        AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
-        pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav, 0.001);
+        pthermo->GetPrimitive(w.at(k, j, i));
+        pthermo->Extrapolate_inplace(pcoord->dx1f(i), "isothermal", grav,
+                                     0.001);
         // add noise
-        air.w[IVY] = 0.00001 * distribution(generator);
-        air.w[IVZ] = 0.00001 * distribution(generator);
+        w(IVY, k, j, i) = 0.00001 * distribution(generator);
+        w(IVZ, k, j, i) = 0.00001 * distribution(generator);
       }
     }
 }

examples/2023-Chen-exo3/test_adiabat.cpp

@@ -13,7 +13,6 @@
 #include <application/exceptions.hpp>
 
 // canoe
-#include <air_parcel.hpp>
 #include <impl.hpp>
 #include <index_map.hpp>
 
@@ -23,6 +22,7 @@
 #include <climath/root.h>
 
 // snap
+#include <snap/thermodynamics/atm_thermodynamics.hpp>
 #include <snap/thermodynamics/thermodynamics.hpp>
 
 Real Ps, Ts, grav;
@@ -36,15 +36,16 @@ void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
 
 void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
   auto pthermo = Thermodynamics::GetInstance();
+  auto &w = phydro->w;
 
   for (int k = ks; k <= ke; ++k)
     for (int j = js; j <= je; ++j)
       for (int i = is; i <= ie; ++i) {
         user_out_var(0, k, j, i) = pthermo->GetTemp(w.at(k, j, i));
-        user_out_var(1, k, j, i) = pthermo->PotentialTemp(w.at(k, j, i), Ps);
+        user_out_var(1, k, j, i) = potential_temp(pthermo, w.at(k, j, i), Ps);
         // msv
         user_out_var(2, k, j, i) =
-            pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i);
+            moist_static_energy(pthermo, w.at(k, j, i), grav * pcoord->x1v(i));
       }
 }
 
@@ -56,17 +57,21 @@ void Mesh::InitUserMeshData(ParameterInput *pin) {
 
 void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   auto pthermo = Thermodynamics::GetInstance();
-
-  AirParcel air(AirParcel::Type::MoleFrac);
+  auto &w = phydro->w;
 
   // construct a reversible adiabat
+  std::vector<Real> yfrac(IVX, 1.);
   for (int k = ks; k <= ke; ++k)
     for (int j = js; j <= je; ++j) {
-      air.w[IPR] = Ps;
-      air.w[IDN] = Ts;
+      pthermo->SetMassFractions<Real>(yfrac.data());
+      pthermo->EquilibrateTP(Ts, Ps);
+
+      // half a grid to cell center
+      pthermo->Extrapolate_inplace(pcoord->dx1f(is) / 2., "reversible", grav);
+
       for (int i = is; i <= ie; ++i) {
-        AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
-        pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav);
+        pthermo->GetPrimitive(w.at(k, j, i));
+        pthermo->Extrapolate_inplace(pcoord->dx1f(i), "reversible", grav);
       }
     }
 }

patches/33.bvals_cc_cpp.patch

@@ -18,7 +18,7 @@ index 223962de..74b3ab79 100644
    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, TypeFlag);
++  cs::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