Add RRTMG-SW (#135)

- Bind RRTMG-SW to C

CMakeLists.txt

@@ -33,10 +33,10 @@ message(STATUS "2. Set up project parameters ...")
 # load custom task
 message(STATUS "Load custom task = ${TASK}")
 if (NOT ${TASK} STREQUAL "")
-  if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/cmake/${TASK}.cmake)
+  if (EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/cmake/examples/${TASK}.cmake)
     message(STATUS "Load custom setting - ${TASK} -")
     message(STATUS "Include ${CMAKE_CURRENT_SOURCE_DIR}/cmake/${TASK}.cmake")
-    include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/${TASK}.cmake)
+    include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/examples/${TASK}.cmake)
   else()
     message(FATAL_ERROR "Cannot find custom file cmake/${TASK}.cmake")
   endif()
@@ -59,6 +59,7 @@ include(${CMAKE_SOURCE_DIR}/cmake/athena.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/application.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/rfm.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/disort.cmake)
+include(${CMAKE_SOURCE_DIR}/cmake/rrtmg.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/pvfmm.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/pybind11.cmake)
 

cmake/examples/rcemip.cmake

@@ -0,0 +1,9 @@
+# configure RCEMIP case
+set(NUMBER_GHOST_CELLS 3)
+set(NVAPOR 1)
+set(NCLOUD 2)
+set(NETCDF ON)
+set(PNETCDF ON)
+set(MPI ON)
+set(RRTMG_SW ON)
+set(RSOLVER lmars)

cmake/parameters.cmake

@@ -95,6 +95,18 @@ else()
   set(DISORT_OPTION "RT_DISORT")
 endif()
 
+if(NOT RRTMG_SW OR NOT DEFINED RRTMG_SW)
+  set(RRTMG_SW_OPTION "NOT_RT_RRTMG_SW")
+else()
+  set(RRTMG_SW_OPTION "RT_RRTMG_SW")
+endif()
+
+if(NOT RRTMG_LW OR NOT DEFINED RRTMG_LW)
+  set(RRTMG_LW_OPTION "NOT_RT_RRTMG_LW")
+else()
+  set(RRTMG_LW_OPTION "RT_RRTMG_LW")
+endif()
+
 if(NOT PVFMM OR NOT DEFINED PVFMM)
   set(PVFMM_OPTION "DISABLE_PVFMM")
 else()

cmake/rrtmg.cmake

@@ -0,0 +1,29 @@
+include(FetchContent)
+
+set(FETCHCONTENT_QUIET FALSE)
+
+macro(set_if_empty _variable)
+  if("${${_variable}}" STREQUAL "")
+    set(${_variable} ${ARGN})
+  endif()
+endmacro()
+
+option(RRTMG_SW "Build RRTMG_SW" OFF)
+
+if(RRTMG_SW)
+  FetchContent_Declare(
+    RRTMG_SW
+    GIT_REPOSITORY https://github.com/chengcli/RRTMG_SW
+    GIT_TAG master)
+  FetchContent_MakeAvailable(RRTMG_SW)
+endif()
+
+option(RRTMG_LW "Build RRTMG_LW" OFF)
+
+if(RRTMG_LW)
+  FetchContent_Declare(
+    RRTMG_LW
+    GIT_REPOSITORY https://github.com/chengcli/RRTMG_LW
+    GIT_TAG master)
+  FetchContent_MakeAvailable(RRTMG_LW)
+endif()

examples/2024-CLi-RCEMIP/CMakeLists.txt

@@ -0,0 +1,11 @@
+
+# 1. Compile rcemip
+
+# 2. Compile hjupiter
+setup_problem(rcemip)
+
+# 2. Copy input files to run directory
+file(GLOB inputs *.py *.inp *.yaml)
+foreach(input ${inputs})
+  execute_process(COMMAND ln -sf ${input} ${CMAKE_BINARY_DIR}/bin/${inp})
+endforeach()

examples/2024-CLi-RCEMIP/rce_small.inp

@@ -0,0 +1,107 @@
+<comment>
+problem     = RCEMPI-1 Small Domain
+
+<job>
+problem_id  = rce1_small    # problem ID: basename of output filenames
+
+<output0>
+file_type   = rst
+dt          = 10.E5
+
+<output1>
+file_type   = hst           # History data dump
+dt          = 2.E4          # time increment between outputs
+
+<output2>
+file_type   = pnetcdf        # Netcdf format output
+variable    = prim          # variables to be output
+dt          = 1.E4          # time increment between outputs
+
+<output3>
+file_type   = pnetcdf        # Netcdf format output
+variable    = uov           # diagnostic variables
+dt          = 1.E4          # time increment
+
+<output3>
+file_type   = pnetcdf        # Netcdf format output
+variable    = diag           # diagnostic variables
+dt          = 1.E4          # time increment
+
+<time>
+cfl_number  = 0.9
+nlim        = -1            # cycle limit
+tlim        = 100.E5
+xorder      = 5             # horizontal reconstruction order
+integrator  = rk3           # integration method
+
+<mesh>
+nx1         = 80          # Number of zones in X1-direction
+x1min       = 0.E3     # minimum value of X1
+x1max       = 34.E3       # maximum value of X1
+ix1_bc      = reflecting  # Inner-X1 boundary condition flag
+ox1_bc      = reflecting  # Outer-X1 boundary condition flag
+
+nx2         = 100         # Number of zones in X2-direction
+x2min       = 0.          # minimum value of X2
+x2max       = 100.E3      # maximum value of X2
+ix2_bc      = periodic    # Inner-X2 boundary condition flag
+ox2_bc      = periodic    # Outer-X2 boundary condition flag
+
+nx3         = 100         # Number of zones in X3-direction
+x3min       = 0.          # minimum value of X3
+x3max       = 100.E3      # maximum value of X3
+ix3_bc      = periodic    # Inner-X3 boundary condition flag
+ox3_bc      = periodic    # Outer-X3 boundary condition flag
+
+<meshblock>
+nx1         = 80
+nx2         = 50
+nx3         = 50
+
+<hydro>
+gamma       = 1.42
+grav_acc1   = -9.79764
+implicit_flag = 1
+
+<species>
+vapor = H2O
+cloud = H2O(c), H2O(p)
+
+<chemistry>
+microphysics_config = rcemip.yaml
+Nc = 1.e8
+Ni = 1.e5
+
+<thermodynamics>
+Rd          = 287.04   # mu = 2.3175 g/mol
+
+eps1        = 8.18    8.18    8.18
+beta1       = 0.      24.845  24.845
+rcp1        = 0.15    0.33    0.33
+Ttriple1    = 273.16
+Ptriple1    = 611.7
+
+sa.relax    = 1.0
+sa.max_iter = 4
+sa.ftol     = 1.e-2
+
+#<forcing>
+#packages = fix_bot_temp
+
+<problem>
+radius      = 6471.e3
+xCO2        = 348.e-6
+xCH4        = 1650.e-12
+xN2O        = 306.e-12
+
+zt          = 15.e3
+zq1         = 4.e3
+zq2         = 7.5e3
+
+qt          = 1.e-14
+q0          = 12.e-3
+
+P0          = 1.0148e5
+T0          = 295.
+gamma       = 0.0067
+diagnostics = div, curl, b, mean, div_h, tempa, presa, hydroflux

examples/2024-CLi-RCEMIP/rcemip.cpp

@@ -0,0 +1,120 @@
+// athena
+#include <athena/athena.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>
+
+// climath
+#include <climath/core.h>
+
+// snap
+#include <snap/thermodynamics/thermodynamics.hpp>
+
+//! variables used in initial condition
+Real grav, P0, T0, gamma;
+Real xCO2, xCH4, xN2O;
+Real zt, zq1, zq2;
+Real qt, q0;
+int iH2O = 1, iCO2 = 0, iCH4 = 1, iN2O = 2, iO3 = 3;
+
+void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
+  AllocateUserOutputVariables(5);
+  SetUserOutputVariableName(0, "temp");
+  SetUserOutputVariableName(1, "theta");
+  SetUserOutputVariableName(2, "thetav");
+  SetUserOutputVariableName(3, "mse");
+  SetUserOutputVariableName(4, "rh_H2O");
+
+  AllocateRealUserMeshBlockDataField(1);
+
+  // CO2, CH4, N2O, O3 in mole fraction
+  ruser_meshblock_data[0].NewAthenaArray(4, ncells1);
+  for (int i = is; i <= ie; ++i) {
+    ruser_meshblock_data[0](iCO2, i) = xCO2;
+    ruser_meshblock_data[0](iCH4, i) = xCH4;
+    ruser_meshblock_data[0](iN2O, i) = xN2O;
+  }
+}
+
+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);
+        user_out_var(4, k, j, i) = pthermo->RelativeHumidity(this, 1, k, j, i);
+      }
+}
+
+void Mesh::InitUserMeshData(ParameterInput *pin) {
+  grav = -pin->GetReal("hydro", "grav_acc1");
+
+  xCO2 = pin->GetReal("problem", "xCO2");
+  xCH4 = pin->GetReal("problem", "xCH4");
+  xN2O = pin->GetReal("problem", "xN2O");
+
+  zt = pin->GetReal("problem", "zt");
+  zq1 = pin->GetReal("problem", "zq1");
+  zq2 = pin->GetReal("problem", "zq2");
+
+  qt = pin->GetReal("problem", "qt");
+  q0 = pin->GetReal("problem", "q0");
+
+  P0 = pin->GetReal("problem", "P0");
+  T0 = pin->GetReal("problem", "T0");
+  gamma = pin->GetReal("problem", "gamma");
+}
+
+void MeshBlock::ProblemGenerator(ParameterInput *pin) {
+  auto pthermo = Thermodynamics::GetInstance();
+
+  Real Rd = pthermo->GetRd();
+  Real Tv0 = T0 * (1. + 0.608 * q0);
+  Real Tvt = Tv0 - gamma * zt;
+  Real Pt = P0 * pow(Tvt / Tv0, grav / (Rd * gamma));
+
+  // O3 parameters
+  Real g1 = 3.6478;  // ppmv hPa−g2
+  Real g2 = 0.83209;
+  Real g3 = 11.3515;  // hPa
+
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j)
+      for (int i = is; i <= ie; ++i) {
+        Real z = pcoord->x1v(i);
+        Real Tv;
+
+        if (z < zt) {
+          phydro->w(iH2O, k, j, i) = q0 * exp(-z / zq1) * exp(-sqr(z / zq2));
+          phydro->w(IPR, k, j, i) =
+              P0 * pow((Tv0 - gamma * z) / Tv0, grav / (Rd * gamma));
+          Tv = Tv0 - gamma * z;
+        } else {
+          phydro->w(iH2O, k, j, i) = qt;
+          phydro->w(IPR, k, j, i) = Pt * exp(-grav * (z - zt) / (Rd * Tvt));
+          Tv = Tvt;
+        }
+
+        phydro->w(IDN, k, j, i) = phydro->w(IPR, k, j, i) / (Rd * Tv);
+
+        // pa -> hpa
+        Real pre = phydro->w(IPR, k, j, i) / 100.;
+        // ppmv -> mole fraction
+        ruser_meshblock_data[0](iO3, i) =
+            1.e-6 * g1 * pow(pre, g2) * exp(-pre / g3);
+      }
+
+  peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie,
+                             js, je, ks, ke);
+}

examples/2024-CLi-RCEMIP/rcemip.yaml

@@ -0,0 +1,9 @@
+name: giant planet chemistry model
+
+microphysics:
+  - water-system
+
+water-system:
+  scheme: Kessler94
+  dependent-species: [vapor.H2O, cloud.H2O(c), cloud.H2O(p)]
+  parameters: {autoconversion: 1.e-3, accretion: 0.0, evaporation: 1.e-3, sedimentation: -20.}

examples/CMakeLists.txt

@@ -44,3 +44,7 @@ endif()
 if (${TASK} STREQUAL "hjupiter")
   add_subdirectory(2024-XZhang-cloud-rt)
 endif()
+
+if (${TASK} STREQUAL "rcemip")
+  add_subdirectory(2024-CLi-RCEMIP)
+endif()