Revise thermodynamics step-1 (#154)

- This PR completes the first stage of revising the thermodynamics. - It uses Cantera to do I/O and thermodynamic calculations - It compiles and has been tested in test_thermodynamcis - But it does not work for any real case at this moment - To be completed in step-2
chengcli · Jun 27, 2024 · b428c62 · b428c62
1 parent fb616f5
commit b428c62
Show file tree

Hide file tree

Showing 37 changed files with 917 additions and 458 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -58,11 +58,14 @@ endif()
 message(STATUS "Include ${CMAKE_CURRENT_SOURCE_DIR}/cmake/parameters.cmake")
 include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/parameters.cmake)
 find_package(Eigen3 REQUIRED)
+find_package(Cantera REQUIRED)
 
 ## 3. set up project system libraries ##
 message(STATUS "3. Set up system libraries")
 #message(STATUS "Include ${CMAKE_SOURCE_DIR}/cmake/athena.cmake")
-include(${CMAKE_SOURCE_DIR}/cmake/yamlpp.cmake)
+if (NOT ${CANTERA_FOUND} STREQUAL "TRUE")
+  include(${CMAKE_SOURCE_DIR}/cmake/yamlpp.cmake)
+endif()
 include(${CMAKE_SOURCE_DIR}/cmake/gtest.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/athena.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/application.cmake)

diff --git a/cmake/macros/macro_setup_test.cmake b/cmake/macros/macro_setup_test.cmake
@@ -22,7 +22,9 @@ macro(setup_test namel)
             ${NETCDF_INCLUDES}
             ${PNETCDF_INCLUDE_DIR}
             ${OpenMP_CXX_INCLUDE_DIR}
-            ${FFTW_INCLUDE_DIRS})
+            ${FFTW_INCLUDE_DIRS}
+            ${CANTERA_INCLUDE_DIR}
+            ${CANTERA_INCLUDE_DIR}/cantera/ext) # for yaml-cpp
 
   target_link_libraries(
     ${namel}.${buildl} gtest_main $<$<BOOL:${PVFMM}>:pvfmmStatic>

diff --git a/cmake/modules/FindCantera.cmake b/cmake/modules/FindCantera.cmake
@@ -1,12 +1,15 @@
 # * Find Cantera Find the Cantera include and library
 #
+# Define the following variables
+#
 # CANTERA_INCLUDE_DIR  - where to find cantera/kinetics.h CANTERA_LIBRARY      -
 # link library CANTERA_FOUND        - True if Cantera found
 #
-# Normal usage would be find_package (Cantera REQUIRED)
-# target_include_directories (${CANTERA_INCLUDE_DIR}) target_link_libraries
-# (${CANTERA_LIBRARY})
+# Normal usage would be
 #
+# find_package (Cantera REQUIRED) target_include_directories
+# (${CANTERA_INCLUDE_DIR}) target_link_libraries (${CANTERA_LIBRARY})
+
 find_path(CANTERA_INCLUDE_DIR cantera/kinetics.h HINTS $ENV{HOME}/opt/include
                                                        $ENV{CANTERA_DIR})
 
@@ -22,4 +25,28 @@ if(CANTERA_INCLUDE_DIR)
   include(FindPackageHandleStandardArgs)
   find_package_handle_standard_args(Cantera DEFAULT_MSG CANTERA_LIBRARY
                                     CANTERA_INCLUDE_DIR)
+
+  set(SOURCE_PATH "${CANTERA_INCLUDE_DIR}/cantera/ext/")
+  set(LINK_PATH "ext1")
+  if(UNIX)
+    execute_process(
+      COMMAND ln -sf ${SOURCE_PATH} ${LINK_PATH}
+      WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+      RESULT_VARIABLE result)
+  elseif(WIN32)
+    execute_process(
+      COMMAND cmd.exe /c mklink /D ${LINK_PATH} ${SOURCE_PATH}
+      WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+      RESULT_VARIABLE result)
+  else()
+    message(
+      FATAL_ERROR "Symbolic link creation is not supported on this platform.")
+  endif()
+
+  if(result)
+    message(FATAL_ERROR "Failed to create symbolic link: ${LINK_PATH}")
+  else()
+    message(STATUS "Symbolic link created: ${LINK_PATH} -> ${SOURCE_PATH}")
+  endif()
+  include_directories(${CMAKE_CURRENT_BINARY_DIR}/ext1)
 endif()
diff --git a/examples/2019-Li-snap/bryan-thermo.yaml b/examples/2019-Li-snap/bryan-thermo.yaml
@@ -0,0 +1,34 @@
+
+phases:
+  - name: gas
+    thermo: ideal-moist
+    species: [dry, H2O, H2O(l)]
+    kinetics: condensation
+
+species:
+  - name: dry
+    composition: {O: 0.42, N: 1.56, Ar: 0.01}
+    thermo:
+      model: constant-cp
+      cp0: 29.1 J/mol/K
+
+  - name: H2O
+    composition: {H: 2, O: 1}
+    thermo:
+      model: constant-cp
+      cp0: 37.4 J/mol/K
+
+  - name: H2O(l)
+    composition: {H: 2, O: 1}
+    thermo:
+      model: constant-cp
+      cp0: 75.3 J/mol/K
+      h0: -28.5 kJ/mol
+    equation-of-state:
+      model: constant-volume
+      density: 1.0 g/cm^3
+
+reactions:
+  - equation: H2O <=> H2O(l)
+    type: nucleation
+    rate-constant: {formula: ideal, T3: 273.16, P3: 611.7, beta: 24.845, delta: 4.986, minT: 273.16}
diff --git a/examples/2019-Li-snap/bryan.cpp.new b/examples/2019-Li-snap/bryan.cpp.new
@@ -0,0 +1,164 @@
+/* -------------------------------------------------------------------------------------
+ * SNAP Example Program
+ *
+ * Contributer:
+ * Cheng Li, University of Michigan
+ *
+ * Year: 2023
+ * Contact: [email protected]
+ * Reference: Bryan and Fritsch, 2002
+ * -------------------------------------------------------------------------------------
+ */
+
+// athena
+#include <athena/athena.hpp>
+#include <athena/athena_arrays.hpp>
+#include <athena/bvals/bvals.hpp>
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/eos/eos.hpp>
+#include <athena/field/field.hpp>
+#include <athena/hydro/hydro.hpp>
+#include <athena/mesh/mesh.hpp>
+#include <athena/parameter_input.hpp>
+
+// application
+#include <application/application.hpp>
+#include <application/exceptions.hpp>
+
+// canoe
+#include <impl.hpp>
+#include <index_map.hpp>
+
+// climath
+#include <climath/core.h>
+#include <climath/interpolation.h>
+
+#include <climath/root.hpp>
+
+// snap
+#include <snap/thermodynamics/atm_thermodynamics.hpp>
+
+// special includes
+#include "bryan_vapor_functions.hpp"
+
+int iH2O, iH2Oc;
+Real p0, grav;
+
+void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
+  AllocateUserOutputVariables(7);
+  SetUserOutputVariableName(0, "temp");
+  SetUserOutputVariableName(1, "theta");
+  SetUserOutputVariableName(2, "theta_v");
+  SetUserOutputVariableName(3, "mse");
+  SetUserOutputVariableName(4, "theta_e");
+  SetUserOutputVariableName(5, "rh");
+  SetUserOutputVariableName(6, "qtol");
+}
+
+void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
+  auto pthermo = Thermodynamics::GetInstance();
+
+  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(this, k, j, i);
+        user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, p0, k, j, i);
+        // theta_v
+        user_out_var(2, k, j, i) =
+            user_out_var(1, j, i) * pthermo->RovRd(this, k, j, i);
+        // mse
+        user_out_var(3, k, j, i) =
+            pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i);
+        /* theta_e
+        user_out_var(4, k, j, i) =
+            pthermo->EquivalentPotentialTemp(this, p0, iH2O, k, j, i);
+        user_out_var(5, k, j, i) =
+            pthermo->RelativeHumidity(this, iH2O, k, j, i);*/
+        // total mixing ratio
+        auto &&air = AirParcelHelper::gather_from_primitive(this, k, j, i);
+        user_out_var(6, k, j, i) = air.w[iH2O] + air.c[iH2Oc];
+      }
+}
+
+void Mesh::InitUserMeshData(ParameterInput *pin) {
+  auto pindex = IndexMap::GetInstance();
+
+  grav = -pin->GetReal("hydro", "grav_acc1");
+  p0 = pin->GetReal("problem", "p0");
+
+  // index
+  iH2O = pindex->GetVaporId("H2O");
+  iH2Oc = pindex->GetCloudId("H2O(c)");
+}
+
+void MeshBlock::ProblemGenerator(ParameterInput *pin) {
+  auto pthermo = Thermodynamics::GetInstance();
+
+  Real Ps = p0;
+  Real Ts = pin->GetReal("problem", "Ts");
+
+  Real xc = pin->GetReal("problem", "xc");
+  Real zc = pin->GetReal("problem", "zc");
+  Real xr = pin->GetReal("problem", "xr");
+  Real zr = pin->GetReal("problem", "zr");
+  Real dT = pin->GetReal("problem", "dT");
+  Real qt = pin->GetReal("problem", "qt");
+
+  AirParcel air(AirParcel::Type::MassFrac);
+  air.w[iH2O] = qt;
+  air.c[iH2Oc] = 0.;
+
+  air.ToMoleFraction();
+  qt = air.w[iH2O];
+
+  // construct a reversible adiabat
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j) {
+      air.w[iH2O] = qt;
+      air.w[IPR] = Ps;
+      air.w[IDN] = Ts;
+      air.c[iH2Oc] = 0.;
+
+      // half a grid to cell center
+      pthermo->Extrapolate(&air, 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), "reversible", grav);
+      }
+    }
+
+  // add temperature anomaly
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j)
+      for (int i = is; i <= ie; ++i) {
+        Real x1 = pcoord->x1v(i);
+        Real x2 = pcoord->x2v(j);
+        Real L = sqrt(sqr((x2 - xc) / xr) + sqr((x1 - zc) / zr));
+
+        if (L < 1.) {
+          auto &&air = AirParcelHelper::gather_from_conserved(this, k, j, i);
+          air.ToMoleFraction();
+          Real rovrd = get_rovrd(air, pthermo->GetMuRatio());
+          Real temp_v = air.w[IDN] * rovrd;
+          temp_v *= 1. + dT * sqr(cos(M_PI * L / 2.)) / 300.;
+
+          Real temp;
+          int err = root(air.w[IDN], air.w[IDN] + dT, 1.E-8, &temp,
+                         [&pthermo, &air, temp_v](Real temp) {
+                           air.w[IDN] = temp;
+
+                           pthermo->EquilibrateTP(&air);
+                           Real rovrd = get_rovrd(air, pthermo->GetMuRatio());
+
+                           return temp * rovrd - temp_v;
+                         });
+
+          if (err) throw RuntimeError("pgen", "TVSolver doesn't converge");
+
+          air.w[IDN] = temp;
+          pthermo->EquilibrateTP(&air);
+          AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
+        }
+      }
+}
diff --git a/examples/2019-Li-snap/jup-thermo.yaml b/examples/2019-Li-snap/jup-thermo.yaml
@@ -0,0 +1,120 @@
+# Thermodynamic data for Jupiter's atmosphere
+#
+# constant-cp model assumes a constant specific heat, which by default is zero.
+# Parameters that can be specified are `cp0`, `t0`, `h0`, and `s0`.
+# If omitted, t0 = 298.15 K, cp0 = 0, h0 = 0, and s0 = 0 at 101325 Pa (OneAtm).
+#
+# The thermo properties are computed as follows:
+#
+#   h = h0 + cp0 * (t - t0)
+#   s = s0 + cp0 * ln(t/t0)
+
+phases:
+  - name: gas
+    thermo: ideal-moist
+    species: [dry, H2O, NH3, H2S, H2O(l), NH3(l), H2O(s), NH3(s), NH4SH(s)]
+    kinetics: condensation
+
+species:
+  - name: dry
+    composition: {H: 1.686, He: 0.157}
+    thermo:
+      model: constant-cp
+      cp0: 29.1 J/mol/K
+
+  - name: H2O
+    composition: {H: 2, O: 1}
+    thermo:
+      model: constant-cp
+      cp0: 37.4 J/mol/K
+
+  - name: NH3
+    composition: {N: 1, H: 3}
+    thermo:
+      model: constant-cp
+      cp0: 37.4 J/mol/K
+
+  - name: H2S
+    composition: {H: 2, S: 1}
+    thermo:
+      model: constant-cp
+      cp0: 37.4 J/mol/K
+
+  - name: H2O(l)
+    composition: {H: 2, O: 1}
+    thermo:
+      model: constant-cp
+      cp0: 75.3 J/mol/K
+      h0: -28.5 kJ/mol
+    equation-of-state:
+      model: constant-volume
+      density: 1.0 g/cm^3
+
+  - name: NH3(l)
+    composition: {N: 1, H: 3}
+    thermo:
+      model: constant-cp
+      cp0: 80.0 J/mol/K
+      h0: -45.9 kJ/mol
+    equation-of-state:
+      model: constant-volume
+      density: 1.0 g/cm^3
+
+  - name: H2O(s)
+    composition: {H: 2, O: 1}
+    thermo:
+      model: constant-cp
+      cp0: 75.3 J/mol/K
+      h0: -28.5 kJ/mol
+    equation-of-state:
+      model: constant-volume
+      density: 0.917 g/cm^3
+
+  - name: NH3(s)
+    composition: {N: 1, H: 3}
+    thermo:
+      model: constant-cp
+      cp0: 80.0 J/mol/K
+      h0: -45.9 kJ/mol
+    equation-of-state:
+      model: constant-volume
+      density: 0.917 g/cm^3
+
+  - name: NH4SH(s)
+    composition: {N: 1, H: 5, S: 1}
+    thermo:
+      model: constant-cp
+      cp0: 80.0 J/mol/K
+      h0: -105.9 kJ/mol
+    equation-of-state:
+      model: constant-volume
+      density: 1.5 g/cm^3
+
+reactions:
+  - equation: H2O <=> H2O(l)
+    type: nucleation
+    rate-constant: {formula: ideal, T3: 273.16, P3: 611.7, beta: 24.845, delta: 4.986, minT: 273.16}
+
+  - equation: H2O <=> H2O(s)
+    type: nucleation
+    rate-constant: {formula: ideal, T3: 273.16, P3: 611.7, beta: 22.98, delta: 0.52, maxT: 273.16}
+
+  - equation: H2O(l) <=> H2O(s)
+    type: freezing
+    rate-constant: {T3: 273.16}
+
+  - equation: NH3 <=> NH3(l)
+    type: nucleation
+    rate-constant: {formula: ideal, T3: 195.4, P3: 6060., beta: 20.08, delta: 5.62, minT: 195.4}
+
+  - equation: NH3 <=> NH3(s)
+    type: nucleation
+    rate-constant: {formula: ideal, T3: 195.4, P3: 6060., beta: 20.64, delta: 1.43, maxT: 195.4}
+
+  - equation: NH3(l) <=> NH3(s)
+    type: freezing
+    rate-constant: {T3: 195.4}
+
+  - equation: NH3 + H2S <=> NH4SH(s)
+    type: nucleation
+    rate-constant: {formula: lewis}