Add Jupiter CRM case (#183)

- Add Jupiter CRM case
- Multi-species particle sedimentation
- Size-dependent sedimentation
- Side-independent sedimentation

CMakeLists.txt

@@ -64,7 +64,7 @@ include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/parameters.cmake)
 message(STATUS "${PROJECT_NAME}-3. Setting up system libraries")
 find_package(Eigen3 REQUIRED)
 find_package(Cantera)
-#find_package(Torch)
+find_package(Torch)
 
 include(${CMAKE_SOURCE_DIR}/cmake/yamlpp.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/gtest.cmake)

cmake/examples/jcrm.cmake

@@ -4,11 +4,12 @@
 set(NUMBER_GHOST_CELLS 3)
 
 # canoe configure
-set(NVAPOR 2)
-set(NCLOUD 3)
+set(NVAPOR 1)
+set(NCLOUD 1)
+set(NPRECIP 1)
 set(NETCDF ON)
 set(PNETCDF ON)
 set(MPI ON)
 set(EQUATION_OF_STATE ideal_moist)
-# set(RSOLVER lmars)
-set(RSOLVER hllc_transform)
+set(RSOLVER lmars)
+# set(RSOLVER hllc_transform)

cmake/examples/jupiter_h2o.cmake

@@ -1,20 +1,15 @@
 # configure file for test jupiter crm
 
-macro(SET_IF_EMPTY _variable)
-  if("${${_variable}}" STREQUAL "")
-    set(${_variable} ${ARGN})
-  endif()
-endmacro()
-
 # athena variables
-set_if_empty(NUMBER_GHOST_CELLS 3)
+set(NUMBER_GHOST_CELLS 3)
 
 # canoe configure
 set(NVAPOR 1)
 set(NCLOUD 2)
 set(NPHASE_LEGACY 3)
 set(PNETCDF ON)
 set(MPI ON)
+set(EQUATION_OF_STATE ideal_moist)
 set(TASKLIST ImplicitHydroTasks)
-#set_if_empty(RSOLVER hllc_transform)
+# set_if_empty(RSOLVER hllc_transform)
 set(RSOLVER lmars)

cmake/parameters.cmake

@@ -14,6 +14,9 @@ set_if_empty(NSTATIC 0)
 set_if_empty(NINT_PARTICLE_DATA 0)
 set_if_empty(NREAL_PARTICLE_DATA 0)
 
+# dirty
+set_if_empty(NBLOCKS 0)
+
 if(NOT AFFINE OR NOT DEFINED AFFINE)
   set(AFFINE_OPTION "NOT_AFFINE")
 else()

configure.hpp.in

@@ -1,8 +1,7 @@
-#ifndef CONFIGURE_HPP_
-#define CONFIGURE_HPP_
+#pragma once
 
 enum {
-  //cubesphere nblocks data
+  // dirty cubesphere nblocks data
   NBLOCKS = @NBLOCKS@,
 
   // thermodynamics data
@@ -57,6 +56,19 @@ enum {
 // nblocks for exo3 option (USE_NBLOCKS or NOT_USE_NBLOCKS)
 #define @USE_NBLOCKS@
 
-#define ENABLE_FIX
-
-#endif  // CONFIGURE_HPP_
+//#define ENABLE_FIX
+
+// logger
+#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_@LOGGER_LEVEL@
+#define LOG_TRACE(...) SPDLOG_LOGGER_TRACE(__VA_ARGS__)
+#define LOG_DEBUG(...) SPDLOG_LOGGER_DEBUG(__VA_ARGS__)
+#define LOG_INFO(...) SPDLOG_LOGGER_INFO(__VA_ARGS__)
+#define LOG_WARN(...) SPDLOG_LOGGER_WARN(__VA_ARGS__)
+#define LOG_ERROR(...) SPDLOG_LOGGER_ERROR(__VA_ARGS__)
+#define LOG_CRITICAL(...) SPDLOG_LOGGER_CRITICAL(__VA_ARGS__)
+
+#ifdef __CUDACC__
+  #define DISPATCH_MACRO __host__ __device__
+#else
+  #define DISPATCH_MACRO
+#endif

examples/2019-Li-snap/jcrm.cpp

@@ -118,7 +118,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   Real Ps = P0 * pow(Ts / T0, cp / Rd);
 
   // read deep composition
-  std::vector<Real> yfrac(1 + NVAPOR, 0.);
+  std::vector<Real> yfrac(IVX, 0.);
   Real qdry = 1.;
   for (int n = 1; n <= NVAPOR; ++n) {
     std::string name = "q" + pthermo->SpeciesName(n);
@@ -132,6 +132,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   while (iter++ < max_iter) {
     pthermo->SetMassFractions<Real>(yfrac.data());
     pthermo->EquilibrateTP(Ts, Ps);
+    std::cout << "Ts = " << Ts << ", Ps = " << Ps << std::endl;
 
     // stop at just above Z = 0
     int i = is;

examples/2019-Li-snap/jcrm.inp

@@ -64,14 +64,17 @@ grav_acc1       = -24.79
 gamma           = 1.4    # override by thermodynamics
 implicit_flag   = 1
 
+<microphysics>
+particle_radius = 0., 1.e-5
+particle_density = 1.e3, 1.e3
+
 <forcing>
 packages = bot_heating, top_cooling
 bot_heating.flux = 50.
 top_cooling.flux = -50.
 
 <problem>
-thermodynamics_config = jup-thermo-H2O-NH3.yaml
-#microphysics_config = water_ammonia.yaml
+thermodynamics_config = jup-precip.yaml
 
 qH2O        = 0.0500
 qNH3        = 0.0028

examples/2019-Li-snap/jup-precip.yaml

@@ -0,0 +1,55 @@
+
+phases:
+  - name: atm
+    thermo: ideal-moist
+    species: [dry, H2O, H2O(l), H2O(p)]
+    kinetics: condensation
+
+species:
+  - name: dry
+    composition: {H: 1.5, He: 0.15}
+    thermo:
+      model: constant-cp
+      T0: 273.16
+      cp0: 29.1 J/mol/K
+
+  - name: H2O
+    composition: {H: 2, O: 1}
+    thermo:
+      model: constant-cp
+      T0: 273.16
+      cp0: 21.1 J/mol/K
+
+  - name: H2O(l)
+    composition: {H: 2, O: 1}
+    thermo:
+      model: constant-cp
+      T0: 273.16
+      cp0: 62.556 J/mol/K
+      h0: -45.103 kJ/mol  # h0 = - (beta - delta) * R * T0
+    equation-of-state:
+      model: constant-volume
+      molar-volume: 0.0 cm^3/mol
+
+  - name: H2O(p)
+    composition: {H: 2, O: 1}
+    thermo:
+      model: constant-cp
+      T0: 273.16
+      cp0: 62.556 J/mol/K
+      h0: -45.103 kJ/mol  # h0 = - (beta - delta) * R * T0
+    equation-of-state:
+      model: constant-volume
+      molar-volume: 0.0 cm^3/mol
+
+reactions:
+  - equation: H2O <=> H2O(l)
+    type: nucleation
+    rate-constant: {formula: ideal, T3: 273.16, P3: 611.7, beta: 24.845, delta: 4.986, minT: 100.}
+
+  - equation: H2O(l) => H2O(p)
+    rate-constant: {A: 0.001, b: 0}
+
+  - equation: H2O(p) => H2O
+    type: evaporation
+    rate-constant: {A: 0.01, b: 0}

examples/2024-YZhou-jupiter/jupiter_crm.cpp

@@ -150,7 +150,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   Real Ts = T0 - grav / cp * x1min;
   Real Ps = P0 * pow(Ts / T0, cp / Rd);
   Real xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6;
-  Real xNH3 = pin->GetReal("problem", "qNH3.ppmv") / 1.E6;
+  // Real xNH3 = pin->GetReal("problem", "qNH3.ppmv") / 1.E6;
 
   while (iter++ < max_iter) {
     // read in vapors

examples/2024-YZhou-jupiter/jupiter_h2o.cpp

@@ -37,7 +37,7 @@
 #include <snap/thermodynamics/atm_thermodynamics.hpp>
 
 // special includes
-#include <special/giants_enroll_vapor_functions_v1.hpp>
+// #include <special/giants_enroll_vapor_functions_v1.hpp>
 
 Real grav, P0, T0, Tmin, prad, hrate;
 int iH2O;
@@ -53,20 +53,21 @@ 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), P0);
+        user_out_var(1, k, j, i) = potential_temp(pthermo, w.at(k, j, i), P0);
         // theta_v
         user_out_var(2, k, j, i) =
-            user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i);
+            user_out_var(1, k, j, i) * pthermo->RovRd(w.at(k, j, i));
         // mse
         user_out_var(3, 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));
         // relative humidity
-        user_out_var(4, k, j, i) = pthermo->RelativeHumidity(this, 1, k, j, i);
+        user_out_var(4, k, j, i) = relative_humidity(pthermo, w.at(k, j, i))[1];
       }
 }
 
@@ -81,12 +82,7 @@ void Forcing(MeshBlock *pmb, Real const time, Real const dt,
   for (int k = ks; k <= ke; ++k)
     for (int j = js; j <= je; ++j)
       for (int i = is; i <= ie; ++i) {
-        auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i);
-
-        air.ToMoleFraction();
-        Real cv =
-            get_cv_mass(air, 0, pthermo->GetRd(), pthermo->GetCvRatioMass());
-
+        Real cv = pthermo->GetCv(w.at(k, j, i));
         if (w(IPR, k, j, i) < prad) {
           du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv *
                               (1. + 1.E-4 * sin(2. * M_PI * rand() / RAND_MAX));
@@ -95,6 +91,7 @@ void Forcing(MeshBlock *pmb, Real const time, Real const dt,
 }
 
 void Mesh::InitUserMeshData(ParameterInput *pin) {
+  auto pthermo = Thermodynamics::GetInstance();
   grav = -pin->GetReal("hydro", "grav_acc1");
 
   P0 = pin->GetReal("problem", "P0");
@@ -105,8 +102,7 @@ void Mesh::InitUserMeshData(ParameterInput *pin) {
   hrate = pin->GetReal("problem", "hrate") / 86400.;
 
   // index
-  auto pindex = IndexMap::GetInstance();
-  iH2O = pindex->GetVaporId("H2O");
+  iH2O = pthermo->SpeciesIndex("H2O");
   EnrollUserExplicitSourceFunction(Forcing);
 }
 
@@ -117,6 +113,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   app->Log("ProblemGenerator: jupiter_crm");
 
   auto pthermo = Thermodynamics::GetInstance();
+  auto &w = phydro->w;
 
   // mesh limits
   Real x1min = pmy_mesh->mesh_size.x1min;
@@ -137,20 +134,21 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
 
   AirParcel air(AirParcel::Type::MoleFrac);
 
-  // estimate surface temperature and pressure
   Real Ts = T0 - grav / cp * x1min;
   Real Ps = P0 * pow(Ts / T0, cp / Rd);
-  Real xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6;
+  Real yH2O = pin->GetReal("problem", "qH2O.gkg") / 1.E3;
+
+  std::vector<Real> yfrac(IVX, 0.);
+  yfrac[iH2O] = yH2O;
+  yfrac[0] = 1. - yH2O;
 
   while (iter++ < max_iter) {
-    // read in vapors
-    air.w[iH2O] = xH2O;
-    air.w[IPR] = Ps;
-    air.w[IDN] = Ts;
+    pthermo->SetMassFractions<Real>(yfrac.data());
+    pthermo->EquilibrateTP(Ts, Ps);
 
     // stop at just above P0
     for (int i = is; i <= ie; ++i) {
-      pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav);
+      pthermo->Extrapolate_inplace(pcoord->dx1f(i), "pseudo", grav);
       if (air.w[IPR] < P0) break;
     }
 
@@ -169,26 +167,24 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   // construct atmosphere from bottom up
   for (int k = ks; k <= ke; ++k)
     for (int j = js; j <= je; ++j) {
-      air.SetZero();
-      air.w[iH2O] = xH2O;
-      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(&air, pcoord->dx1f(is) / 2., "reversible", grav);
+      pthermo->Extrapolate_inplace(pcoord->dx1f(is) / 2., "reversible", grav);
 
       int i = is;
       for (; i <= ie; ++i) {
-        if (air.w[IDN] < Tmin) break;
-        AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
-        pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5);
+        pthermo->GetPrimitive(w.at(k, j, i));
+        if (pthermo->GetTemp() < Tmin) break;
+        pthermo->Extrapolate_inplace(pcoord->dx1f(i), "pseudo", grav);
       }
 
       // Replace adiabatic atmosphere with isothermal atmosphere if temperature
       // is too low
       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);
       }
     }
 }

src/CMakeLists.txt

@@ -62,10 +62,10 @@ add_subdirectory(special)
 add_subdirectory(diagnostics)
 add_subdirectory(forcing)
 
-if (${Torch_FOUND})
-  add_subdirectory(torch)
-  add_subdirectory(modules)
-endif()
+#if (${Torch_FOUND})
+#  add_subdirectory(torch)
+#  add_subdirectory(modules)
+#endif()
 #add_subdirectory(nbody)
 
 if(UNIX AND NOT APPLE)

src/add_arg.h

@@ -0,0 +1,21 @@
+#define ADD_ARG(T, name)                                                  \
+ public:                                                                  \
+  inline auto name(const T& new_##name) -> decltype(*this) { /* NOLINT */ \
+    this->name##_ = new_##name;                                           \
+    return *this;                                                         \
+  }                                                                       \
+  inline auto name(T&& new_##name) -> decltype(*this) { /* NOLINT */      \
+    this->name##_ = std::move(new_##name);                                \
+    return *this;                                                         \
+  }                                                                       \
+  DISPATCH_MACRO                                                          \
+  inline const T& name() const noexcept { /* NOLINT */                    \
+    return this->name##_;                                                 \
+  }                                                                       \
+  DISPATCH_MACRO                                                          \
+  inline T& name() noexcept { /* NOLINT */                                \
+    return this->name##_;                                                 \
+  }                                                                       \
+                                                                          \
+ private:                                                                 \
+  T name##_ /* NOLINT */

src/checks.cpp

@@ -195,7 +195,7 @@ void fix_eos_cons2prim(MeshBlock* pmb, AthenaArray<Real>& prim,
     }
   }
 
-  Real temp = pthermo->GetTemp(pmb, k, j, ifix - 1);
+  Real temp = pthermo->GetTemp(pmb->phydro->w.at(k, j, ifix - 1));
   for (int i = ifix; i <= iu; ++i) {
     Real z = pcoord->x1v(i) - pcoord->x1v(ifix - 1);
     prim(IDN, k, j, i) =

src/microphysics/kessler94.cpp → src/microphysics/kessler94.cpp.old

@@ -11,7 +11,7 @@
 #include <snap/thermodynamics/thermodynamics.hpp>
 
 // microphysics
-#include <snap/thermodynamics/vapors/water_vapors.hpp>
+//#include <snap/thermodynamics/vapors/water_vapors.hpp>
 
 #include "microphysical_schemes.hpp"