Merge exo3 and add harp

cmake/examples/rtdry.cmake

@@ -1,11 +1,5 @@
 # configuration for euler equations on cubed sphere
 
-macro(SET_IF_EMPTY _variable)
-  if("${${_variable}}" STREQUAL "")
-    set(${_variable} ${ARGN})
-  endif()
-endmacro()
-
 # athena variables
 set_if_empty(NUMBER_GHOST_CELLS 3)
 

cmake/parameters.cmake

@@ -57,9 +57,3 @@ else()
     COMPONENTS C
     REQUIRED)
 endif()
-
-if(NOT NBLOCKS OR NOT DEFINED NBLOCKS)
-  set(USE_NBLOCKS "NOT_USE_NBLOCKS")
-else()
-  set(USE_NBLOCKS "USE_NBLOCKS")
-endif()

configure.h.in

@@ -26,6 +26,3 @@ enum {
 
 // MPI parallelization (MPI_PARALLEL or NOT_MPI_PARALLEL)
 #define @MPI_OPTION@
-
-// nblocks for exo3 option (USE_NBLOCKS or NOT_USE_NBLOCKS)
-#define @USE_NBLOCKS@

examples/2023-Chen-exo3/CMakeLists.txt

@@ -15,7 +15,7 @@ else()
   setup_problem(test_adiabat)
   setup_problem(hs94)
   setup_problem(hot_jupiter)
-  setup_problem(polar_dry)
+  #setup_problem(polar_dry)
 endif()
 
 # 2. Copy input file to run directory

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

@@ -1,14 +1,15 @@
-
-phases:
-  - name: atm
-    thermo: ideal-moist
-    species: [atm]
-    kinetics: condensation
-    reactions: none
+reference-state:
+  Tref: 0.
+  Pref: 1.e5
 
 species:
-  - name: atm
+  - name: dry
     composition: {O: 0.42, N: 1.56, Ar: 0.01}
-    thermo:
-      model: constant-cp
-      cp0: 29.1 J/mol/K
+    cv_R: 2.5
+
+geometry:
+  type: cartesian
+
+dynamics:
+  equation-of-state:
+    type: ideal-gas

examples/2023-Chen-exo3/hot_jupiter.cpp

@@ -29,16 +29,13 @@
 #include <application/exceptions.hpp>
 
 // canoe
-#include <configure.hpp>
 #include <impl.hpp>
+#include <interface/eos.hpp>
 
 // exo3
 #include <exo3/cubed_sphere.hpp>
 #include <exo3/cubed_sphere_utility.hpp>
 
-// snap
-#include <snap/thermodynamics/thermodynamics.hpp>
-
 #define _sqr(x) ((x) * (x))
 #define _qur(x) ((x) * (x) * (x) * (x))
 
@@ -232,50 +229,24 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   Rp = pin->GetReal("problem", "Rp");
 
   // construct an isothermal atmosphere
-  auto pthermo = Thermodynamics::GetInstance();
-  auto &w = phydro->w;
-  std::vector<Real> yfrac(1, 1.);
+  kintera::ThermoX thermo(pimpl->peos->pthermo->options);
 
-  for (int k = ks; k <= ke; ++k)
-    for (int j = js; j <= je; ++j) {
-      pthermo->SetMassFractions<Real>(yfrac.data());
-      pthermo->EquilibrateTP(Ts, p0);
-
-      int i = is;
-      for (; i <= ie; ++i) {
-        pthermo->GetPrimitive(w.at(k, j, i));
-        pthermo->Extrapolate_inplace(pcoord->dx1f(i), "isothermal", grav,
-                                     0.001);
-
-        // add noise
-        w(IVY, k, j, i) = 10. * distribution(generator);
-        w(IVZ, k, j, i) = 10. * distribution(generator);
-      }
-    }
+  auto temp = Ts * torch::ones({ncells3, ncells2}, torch::kDouble);
+  auto pres = p0 * torch::ones({ncells3, ncells2}, torch::kDouble);
+
+  // half a grid to cell center
+  pres *= exp(-grav * pcoord->dx1f(is) / (2 * Rd * temp));
+
+  auto w = get_all(phydro->w);
+  for (int i = is; i <= ie; ++i) {
+    w[IPR].select(2, i) = pres;
+    w[IDN].select(2, i) = pres / (Rd * temp);
+    pres *= exp(-grav * pcoord->dx1f(i) / (Rd * temp));
+  }
 
-  // Set up perturbation as the initial condition to break the symmetry of the
-  // original initial condition.
-  // Real k3 = 5.;
-
-  // for (int k = ks; k <= ke; ++k) {
-  //   for (int j = js; j <= je; ++j) {
-  //     for (int i = is; i <= ie; ++i) {
-  //       Real theta, phi;
-  //       pexo3->GetLatLon(&theta, &phi, k, j, i);
-  //       Real z = pcoord->x1v(i)-Rp;
-  //       // Set-up isothermal atmosphere
-  //       Real temp = Ts;
-  //       phydro->w(IPR,k,j,i) = p0*exp(-z*grav/(temp*Rd));
-  //       phydro->w(IDN,k,j,i) = phydro->w(IPR,k,j,i)/
-  //            (Rd* ( temp + 20.*(distribution(generator) - 0.5) * (1. +
-  //            cos(k3*phi)) ) );
-  //       //phydro->w(IDN,k,j,i) = phydro->w(IPR,k,j,i)/(Rd*temp);
-  //       phydro->w(IM1,k,j,i) = 0.;
-  //       phydro->w(IM2,k,j,i) = 0.;
-  //       phydro->w(IM3,k,j,i) = 0.;
-  //     }
-  //   }
-  // }
+  // add noise
+  w[IVY] += 10. * torch::randn_like(w[IVY]);
+  w[IVZ] += 10. * torch::randn_like(w[IVZ]);
 
   // transfer to conservative variables
   // bcc is cell-centered magnetic fields, it is only a place holder here

examples/2023-Chen-exo3/hot_jupiter.inp

@@ -55,13 +55,13 @@ PressureScaleHeight = 7.56E5
 
 <meshblock>
 nx1        = 32
-nx2        = 20
-nx3        = 20
+nx2        = 40
+nx3        = 40
 
 <hydro>
 grav_acc1   = -8.0
 gamma       = 1.4         # gamma = C_p/C_v
-implicit_flag = 1
+implicit_flag = 9
 
 <thermodynamics>
 Rd          = 3779.       # Gas constant
@@ -71,7 +71,7 @@ z_stra      = 2.E6        # Stratosphere height, m
 Gamma_trop  = 2.E-4       # Lapse rate, K/m
 
 <problem>
-thermodynamics_config = hot-jupiter.yaml
+config_file = earth-dry.yaml
 
 Kt         = 1.5E5            # Damping timescale, s
 Ts         = 1600.            # Surface temperature, K

examples/2023-Chen-exo3/hs94.cpp

@@ -29,16 +29,13 @@
 #include <application/exceptions.hpp>
 
 // canoe
-#include <configure.hpp>
 #include <impl.hpp>
+#include <interface/eos.hpp>
 
 // exo3
 #include <exo3/cubed_sphere.hpp>
 #include <exo3/cubed_sphere_utility.hpp>
 
-// snap
-#include <snap/thermodynamics/thermodynamics.hpp>
-
 #define _sqr(x) ((x) * (x))
 #define _qur(x) ((x) * (x) * (x) * (x))
 
@@ -261,39 +258,42 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   z_iso = pin->GetReal("problem", "z_iso");
 
   // construct an adiabatic atmosphere
-  auto pthermo = Thermodynamics::GetInstance();
-  auto &w = phydro->w;
-  std::vector<Real> yfrac(1, 1.);
+  kintera::ThermoX thermo(pimpl->peos->pthermo->options);
 
-  for (int k = ks; k <= ke; ++k)
-    for (int j = js; j <= je; ++j) {
-      pthermo->SetMassFractions<Real>(yfrac.data());
-      pthermo->EquilibrateTP(Ts, p0);
+  auto temp = Ts * torch::ones({ncells3, ncells2}, torch::kDouble);
+  auto pres = p0 * torch::ones({ncells3, ncells2}, torch::kDouble);
+  auto xfrac = torch::ones({ncells3, ncells2, 1}, torch::kDouble);
 
-      // half a grid to cell center
-      pthermo->Extrapolate_inplace(pcoord->dx1f(is) / 2., "dry", grav);
+  // half a grid to cell center
+  thermo->extrapolate_ad(temp, pres, xfrac, grav, pcoord->dx1f(is) / 2.);
 
-      int i = is;
-      for (; i <= ie; ++i) {
-        if (pcoord->x1v(i) - Rp > z_iso) break;
+  auto w = get_all(phydro->w);
+  int i = is;
+  for (; i <= ie; ++i) {
+    auto conc = thermo->compute("TPX->V", {temp, pres, xfrac});
+    w[IPR].select(2, i) = pres;
+    w[IDN].select(2, i) = thermo->compute("V->D", {conc});
+    w.slice(0, 1, IVX).select(3, i) = thermo->compute("X->Y", {xfrac});
 
-        pthermo->GetPrimitive(w.at(k, j, i));
-        pthermo->Extrapolate_inplace(pcoord->dx1f(i), "dry", grav);
+    if (pcoord->x1v(i) - Rp > z_iso) break;
+    thermo->extrapolate_ad(temp, pres, xfrac, grav, pcoord->dx1f(i));
+  }
 
-        // add noise
-        w(IVY, k, j, i) = 10. * distribution(generator);
-        w(IVZ, k, j, i) = 10. * distribution(generator);
-      }
+  // isothermal extrapolation
+  for (; i <= ie; ++i) {
+    pres *= exp(-grav * pcoord->dx1f(i) / (Rd * temp));
+    auto conc = thermo->compute("TPX->V", {temp, pres, xfrac});
+    w[IPR].select(2, i) = pres;
+    w[IDN].select(2, i) = thermo->compute("V->D", {conc});
+    w.slice(0, 1, IVX).select(3, i) = thermo->compute("X->Y", {xfrac});
+  }
 
-      // construct isothermal atmosphere
-      for (; i <= ie; ++i) {
-        pthermo->GetPrimitive(w.at(k, j, i));
-        pthermo->Extrapolate_inplace(pcoord->dx1f(i), "isothermal", grav);
-      }
-    }
+  // add noise
+  w[IVY] += 1. * torch::randn_like(w[IVY]);
+  w[IVZ] += 1. * torch::randn_like(w[IVZ]);
 
-  peos->PrimitiveToConserved(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/hs94.inp

@@ -71,7 +71,7 @@ Rd         = 287.
 cp         = 1004.
 
 <problem>
-thermodynamics_config = earth-dry.yaml
+config_file = earth-dry.yaml
 
 Omega      = 7.292E-5
 Rp         = 6.371E6

examples/2023-Chen-exo3/test_adiabat.cpp

@@ -14,38 +14,32 @@
 
 // canoe
 #include <impl.hpp>
-#include <index_map.hpp>
-
-// climath
-#include <climath/core.h>
-#include <climath/interpolation.h>
-#include <climath/root.h>
-
-// snap
-#include <snap/thermodynamics/atm_thermodynamics.hpp>
-#include <snap/thermodynamics/thermodynamics.hpp>
+#include <interface/eos.hpp>
 
 Real Ps, Ts, grav;
 
 void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
-  AllocateUserOutputVariables(3);
+  AllocateUserOutputVariables(2);
   SetUserOutputVariableName(0, "temp");
   SetUserOutputVariableName(1, "theta");
-  SetUserOutputVariableName(2, "mse");
 }
 
 void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
-  auto pthermo = Thermodynamics::GetInstance();
   auto &w = phydro->w;
 
+  auto temp = get_temp(pimpl->peos, w);
+  auto pres = get_pres(w);
+  auto chi = (get_gammad() - 1.) / get_gammad();
+  auto theta = temp * pow(Ps / pres, chi);
+
+  auto temp_a = temp.accessor<Real, 3>();
+  auto theta_a = theta.accessor<Real, 3>();
+
   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) = potential_temp(pthermo, w.at(k, j, i), Ps);
-        // msv
-        user_out_var(2, k, j, i) =
-            moist_static_energy(pthermo, w.at(k, j, i), grav * pcoord->x1v(i));
+        user_out_var(0, k, j, i) = temp_a[k][j][i];
+        user_out_var(1, k, j, i) = theta_a[k][j][i];
       }
 }
 
@@ -56,22 +50,25 @@ void Mesh::InitUserMeshData(ParameterInput *pin) {
 }
 
 void MeshBlock::ProblemGenerator(ParameterInput *pin) {
-  auto pthermo = Thermodynamics::GetInstance();
-  auto &w = phydro->w;
+  kintera::ThermoX thermo(pimpl->peos->pthermo->options);
 
-  // construct a reversible adiabat
-  std::vector<Real> yfrac(IVX, 1.);
-  for (int k = ks; k <= ke; ++k)
-    for (int j = js; j <= je; ++j) {
-      pthermo->SetMassFractions<Real>(yfrac.data());
-      pthermo->EquilibrateTP(Ts, Ps);
+  auto temp = Ts * torch::ones({ncells3, ncells2}, torch::kDouble);
+  auto pres = Ps * torch::ones({ncells3, ncells2}, torch::kDouble);
+  auto xfrac = torch::ones({ncells3, ncells2, 1}, torch::kDouble);
 
-      // half a grid to cell center
-      pthermo->Extrapolate_inplace(pcoord->dx1f(is) / 2., "reversible", grav);
+  // half a grid to cell center
+  thermo->extrapolate_ad(temp, pres, xfrac, grav, pcoord->dx1f(is) / 2.);
 
-      for (int i = is; i <= ie; ++i) {
-        pthermo->GetPrimitive(w.at(k, j, i));
-        pthermo->Extrapolate_inplace(pcoord->dx1f(i), "reversible", grav);
-      }
-    }
+  auto w = get_all(phydro->w);
+  for (int i = is; i <= ie; ++i) {
+    auto conc = thermo->compute("TPX->V", {temp, pres, xfrac});
+    w[IPR].select(2, i) = pres;
+    w[IDN].select(2, i) = thermo->compute("V->D", {conc});
+    w.slice(0, 1, IVX).select(3, i) = thermo->compute("X->Y", {xfrac});
+    thermo->extrapolate_ad(temp, pres, xfrac, grav, pcoord->dx1f(i));
+  }
+
+  // Change primitive variables to conserved variables
+  peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie,
+                             js, je, ks, ke);
 }

examples/2023-Chen-exo3/test_adiabat.inp

@@ -53,10 +53,10 @@ nx3        = 48
 grav_acc1   = -9.8
 gamma       = 1.4         # gamma = C_p/C_v
 implicit_flag   = 1
-
-<thermodynamics>
 Rd          = 287.
 
 <problem>
+config_file = earth-dry.yaml
+
 Ps    = 1.E5
 Ts    = 300.

src/exo3/cs_find_block_id.cpp

@@ -80,7 +80,7 @@ int CubedSphere::FindBlockID(LogicalLocation const &loc) {
 void CubedSphere::GetLocalIndex(int *lv2_lx2, int *lv2_lx3, int *local_lx2,
                                 int *local_lx3, int *bound_lim,
                                 LogicalLocation const &loc) {
-#ifdef USE_NBLOCKS
+#if NBLOCKS > 0
   // Updated method, need to manually setup in configure.h, allow 6*n^2 blocks
   *bound_lim = (int)(sqrt(NBLOCKS / 6) - 0.5);
   // Find relative location within block