Merge minichem 1

README

@@ -4,6 +4,13 @@ when vapor pressure are too low (~1e-10 Pa) at low temperatuers, the saturationa
 
 when vapor pressures are too high (1e5 Pa) eg. for CO2, the saturation adjustment fails immediately (all nan)
 
-exo can use cores 24/96/384/1536
-also permit 216/600/864/1176/ if specify NBLOCKS in the .cmake file
+exo can use cores 24/96/384/1536: 6*4^n
+also permit 216/600/864/1176/1536/1944/2400... (24*n^2) if specify NBLOCKS in the .cmake file
+That is, if meshblock is x, the mesh should be x*4*n, and x*6*n, for example, if x is 10, mesh could be 40/60, or 240/360 etc.
+meshblock = 5 is not very efficient, maybe meshblock is 8 or higher
 
+
+gammal in svp is  delta[n] = (cp_ratio_mass[n] - cp_ratio_mass[i]) * mu_ratio[i] / (1. - 1. / gammad);
+
+cp_ratio_mass is rcp
+mu_ratio is eps

cmake/examples/ice_h2s.cmake

@@ -0,0 +1,20 @@
+# 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)
+
+# canoe configure
+set(NVAPOR 1)
+set(NCLOUD 2)
+set(NPHASE_LEGACY 3)
+set(PNETCDF ON)
+set(MPI ON)
+set(TASKLIST ImplicitHydroTasks)
+#set_if_empty(RSOLVER hllc_transform)
+set(RSOLVER lmars)

cmake/examples/minichem.cmake

@@ -10,10 +10,8 @@ endmacro()
 set_if_empty(NUMBER_GHOST_CELLS 3)
 
 set(MINICHEM ON)
-set(NCHEMISTRY 13)
-#set(NVAPOR 1)
-#set(NCLOUD 2)
-#set(NPHASE_LEGACY 3)
+set(NCHEM 13)
+# set(NVAPOR 1) set(NCLOUD 2) set(NPHASE_LEGACY 3)
 set(NETCDF ON)
 set(PNETCDF ON)
 set(MPI ON)

cmake/examples/rtcldexo.cmake

@@ -9,14 +9,14 @@ endmacro()
 # athena variables
 set_if_empty(NUMBER_GHOST_CELLS 3)
 
-set(CUBED_SPHERE ON) 
-#set(NBLOCKS 600)
+#set(CUBED_SPHERE ON) 
+#set(NBLOCKS 2400)
 #set(COORDINATE_SYSTEM "spherical_polar")
-set(NVAPOR 1)
-set(NCLOUD 2)
+set(NVAPOR 2)
+set(NCLOUD 4)
 set(NPHASE_LEGACY 3)
 set(NETCDF ON)
-#set(PNETCDF OFF)
+set(PNETCDF ON)
 set(MPI ON)
 set(DISORT ON)
 set(PYTHON_BINDINGS ON)

cmake/examples/rtdry.cmake

@@ -3,11 +3,8 @@
 # athena variables
 set_if_empty(NUMBER_GHOST_CELLS 3)
 
-set(CUBED_SPHERE ON)
-set(NBLOCKS 864)
-# set(COORDINATE_SYSTEM "spherical_polar")
-set(NETCDF OFF)
-set(PNETCDF ON)
+set(NETCDF ON)
+set(PNETCDF OFF)
 set(MPI ON)
 set(DISORT ON)
 set(PYTHON_BINDINGS ON)

cmake/minichem.cmake

@@ -7,4 +7,4 @@ set(patch_command
 set(PACKAGE_NAME minichem)
 set(REPO_URL "https://github.com/chengcli/mini_chem")
 set(REPO_TAG "3372a500f038f83228c9f8f944b3fb6b2dedc572")
-add_package(${PACKAGE_NAME} ${REPO_URL} ${REPO_TAG} "${patch_commend}" OFF)
+add_package(${PACKAGE_NAME} ${REPO_URL} ${REPO_TAG} "${patch_command}" OFF)

cmake/modules/FindKintera.cmake

@@ -4,8 +4,6 @@
 #
 # KINTERA_LIBRARY:      Link these libraries when using KINTERA
 #
-# VAPORS_LIBRARY:       Link these libraries when using KINTERA
-#
 # KINTERA_FOUND:        True if kintera found
 #
 # Normal usage would be:
@@ -14,7 +12,7 @@
 #
 # include_directories(${KINTERA_INCLUDE_DIR})
 #
-# target_link_libraries(${KINTERA_LIBRARY} ${VAPORS_LIBRARY})
+# target_link_libraries(${KINTERA_LIBRARY})
 
 include(FindPackageHandleStandardArgs)
 
@@ -102,22 +100,11 @@ find_library(
         $ENV{KINTERA_DIR}/lib
         $ENV{KINTERA_ROOT}/lib)
 
-find_library(
-  VAPORS_LIBRARY vapors_release
-  HINTS ${kintera_lib_dir}
-        /opt/homebrew/lib
-        /usr/lib/x86_64-linux-gnu/
-        KINTERA_DIR/lib
-        KINTERA_LIB
-        $ENV{KINTERA_LIB}
-        $ENV{KINTERA_DIR}/lib
-        $ENV{KINTERA_ROOT}/lib)
-
 if(${CUDAToolKit_FOUND})
-  set(kintera_required_vars VAPORS_LIBRARY KINTERA_LIBRARY KINTERA_CUDA_LIBRARY
+  set(kintera_required_vars KINTERA_LIBRARY KINTERA_CUDA_LIBRARY
                             KINTERA_INCLUDE_DIR)
 else()
-  set(kintera_required_vars VAPORS_LIBRARY KINTERA_LIBRARY KINTERA_INCLUDE_DIR)
+  set(kintera_required_vars KINTERA_LIBRARY KINTERA_INCLUDE_DIR)
 endif()
 
 mark_as_advanced(${kintera_required_vars})

cmake/modules/FindSnap.cmake

@@ -4,8 +4,6 @@
 #
 # SNAP_LIBRARY:      Link these libraries when using SNAP
 #
-# BC_LIBRARY:         Link these libraries when using SNAP
-#
 # SNAP_FOUND:        True if snap found
 #
 # Normal usage would be:
@@ -14,7 +12,7 @@
 #
 # include_directories(${SNAP_INCLUDE_DIR})
 #
-# target_link_libraries(${SNAP_LIBRARY} ${BC_LIBRARY})
+# target_link_libraries(${SNAP_LIBRARY})
 
 include(FindPackageHandleStandardArgs)
 
@@ -100,22 +98,10 @@ find_library(
         $ENV{SNAP_DIR}/lib
         $ENV{SNAP_ROOT}/lib)
 
-find_library(
-  BC_LIBRARY bc_release
-  HINTS ${snap_lib_dir}
-        /opt/homebrew/lib
-        /usr/lib/x86_64-linux-gnu/
-        SNAP_DIR/lib
-        SNAP_LIB
-        $ENV{SNAP_LIB}
-        $ENV{SNAP_DIR}/lib
-        $ENV{SNAP_ROOT}/lib)
-
 if(${CUDAToolKit_FOUND})
-  set(snap_required_vars BC_LIBRARY SNAP_LIBRARY SNAP_CUDA_LIBRARY
-                         SNAP_INCLUDE_DIR)
+  set(snap_required_vars SNAP_LIBRARY SNAP_CUDA_LIBRARY SNAP_INCLUDE_DIR)
 else()
-  set(snap_required_vars BC_LIBRARY SNAP_LIBRARY SNAP_INCLUDE_DIR)
+  set(snap_required_vars SNAP_LIBRARY SNAP_INCLUDE_DIR)
 endif()
 
 mark_as_advanced(${snap_required_vars})

examples/2024-XZhang-cloud-rt/CMakeLists.txt

@@ -1,11 +1,13 @@
 #if (${NVAPOR} EQUAL 1 AND ${NCLOUD} EQUAL 2 AND ${NPHASE_LEGACY} EQUAL 3)
 if (${NVAPOR} EQUAL 1)
   #setup_problem(rt_silicate)
-  setup_problem(beta_rtsi)
+  setup_problem(ice_h2s)
+  #setup_problem(beta_rtsi)
   #setup_problem(rt_h2o)
 elseif (${NVAPOR} EQUAL 2)
   #setup_problem(rt_h2o_nh3)
-  setup_problem(jupiter_crm)
+  #setup_problem(jupiter_crm)
+   setup_problem(rtcldu3d)
 elseif (${NVAPOR} EQUAL 3)
   setup_problem(beta_rtjup3cld)
 else()

examples/2024-XZhang-cloud-rt/ice_h2s.cpp

@@ -0,0 +1,195 @@
+/* -------------------------------------------------------------------------------------
+ * SNAP Example Program
+ *
+ * Contributer:
+ * Cheng Li, University of Michigan
+ *
+ * Year: 2023
+ * Contact: chengcli@umich.edu
+ * Reference: Test Jupiter CRM
+ * -------------------------------------------------------------------------------------
+ */
+
+// 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>
+
+// snap
+#include <snap/thermodynamics/atm_thermodynamics.hpp>
+
+// special includes
+#include "svp_functions.hpp"
+
+Real grav, P0, T0, Tmin, prad, hrate;
+int iH2S;
+
+void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
+  AllocateUserOutputVariables(4 + NVAPOR);
+  SetUserOutputVariableName(0, "temp");
+  SetUserOutputVariableName(1, "theta");
+  SetUserOutputVariableName(2, "thetav");
+  SetUserOutputVariableName(3, "mse");
+  SetUserOutputVariableName(4, "rh_H2S");
+}
+
+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, k, 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);
+        // relative humidity
+        user_out_var(4, k, j, i) = pthermo->RelativeHumidity(this, 1, k, j, i);
+      }
+}
+
+void Forcing(MeshBlock *pmb, Real const time, Real const dt,
+             AthenaArray<Real> const &w, AthenaArray<Real> const &r,
+             AthenaArray<Real> const &bcc, AthenaArray<Real> &du,
+             AthenaArray<Real> &s) {
+  int is = pmb->is, js = pmb->js, ks = pmb->ks;
+  int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
+  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) {
+        auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i);
+
+        air.ToMoleFraction();
+        Real cv =
+            get_cv_mass(air, 0, pthermo->GetRd(), pthermo->GetCvRatioMass());
+
+//        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));
+//        }
+      }
+}
+
+void Mesh::InitUserMeshData(ParameterInput *pin) {
+  grav = -pin->GetReal("hydro", "grav_acc1");
+
+  P0 = pin->GetReal("problem", "P0");
+  T0 = pin->GetReal("problem", "T0");
+
+  Tmin = pin->GetReal("problem", "Tmin");
+//  prad = pin->GetReal("problem", "prad");
+//  hrate = pin->GetReal("problem", "hrate") / 86400.;
+
+  // index
+  auto pindex = IndexMap::GetInstance();
+  iH2S = pindex->GetVaporId("H2S");
+  EnrollUserExplicitSourceFunction(Forcing);
+}
+
+void MeshBlock::ProblemGenerator(ParameterInput *pin) {
+  srand(Globals::my_rank + time(0));
+
+  Application::Logger app("main");
+  app->Log("ProblemGenerator: jupiter_crm");
+
+  auto pthermo = Thermodynamics::GetInstance();
+
+  // mesh limits
+  Real x1min = pmy_mesh->mesh_size.x1min;
+  Real x1max = pmy_mesh->mesh_size.x1max;
+
+  // request temperature and pressure
+  app->Log("request T", T0);
+  app->Log("request P", P0);
+
+  // thermodynamic constants
+  Real gamma = pin->GetReal("hydro", "gamma");
+  Real Rd = pthermo->GetRd();
+  Real cp = gamma / (gamma - 1.) * Rd;
+
+  // set up an adiabatic atmosphere
+  int max_iter = 400, iter = 0;
+  Real Ttol = pin->GetOrAddReal("problem", "init_Ttol", 0.01);
+
+  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 xH2S = pin->GetReal("problem", "qH2S.ppmv") / 1.E6;
+
+  while (iter++ < max_iter) {
+    // read in vapors
+    air.w[iH2S] = xH2S;
+    air.w[IPR] = Ps;
+    air.w[IDN] = Ts;
+
+    // stop at just above P0
+    for (int i = is; i <= ie; ++i) {
+      pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav);
+      if (air.w[IPR] < P0) break;
+    }
+
+    // make up for the difference
+    Ts += T0 - air.w[IDN];
+    if (std::abs(T0 - air.w[IDN]) < Ttol) break;
+
+    app->Log("Iteration #", iter);
+    app->Log("T", air.w[IDN]);
+  }
+
+  if (iter > max_iter) {
+    throw RuntimeError("ProblemGenerator", "maximum iteration reached");
+  }
+
+  // construct atmosphere from bottom up
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j) {
+      air.SetZero();
+      air.w[iH2S] = xH2S;
+      air.w[IPR] = Ps;
+      air.w[IDN] = Ts;
+
+      // half a grid to cell center
+      pthermo->Extrapolate(&air, 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->Extrapolate(&air, 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);
+      }
+    }
+}