Refactor Thermodynamics and Inversion (#137)

- refactor thermodynamics
- refactor inversion
- may contain unfixed errors

.github/workflows/mac.yml

@@ -358,5 +358,5 @@ jobs:
         run: ./combine.py -o test
 
       - name: Test
-        working-directory: ${{github.workspace}}/examples/2023-jupiter-mwr-eq
+        working-directory: ${{github.workspace}}/examples/2024-CLi-juno-mwr
         run: echo "cheers"

.github/workflows/main.yml

@@ -705,5 +705,5 @@ jobs:
         run: python3 combine.py -o test
 
       - name: Test
-        working-directory: ${{github.workspace}}/examples/2023-jupiter-mwr-eq
+        working-directory: ${{github.workspace}}/examples/2024-CLi-juno-mwr
         run: echo "cheers"

CMakeLists.txt

@@ -20,8 +20,16 @@ else()
   message(STATUS "")
 endif()
 
+# load all modules
 LIST(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/modules/)
 
+# load all macros
+FILE(GLOB _macro_files "${CMAKE_SOURCE_DIR}/cmake/macros/*.cmake")
+FOREACH(_file ${_macro_files})
+  MESSAGE(STATUS "Include ${_file}")
+  INCLUDE(${_file})
+ENDFOREACH()
+
 ## 1. set up compiler flags ##
 message(STATUS "1. Set up project compiler flags ...")
 message(STATUS "Include ${CMAKE_SOURCE_DIR}/cmake/compilers.cmake")
@@ -35,10 +43,10 @@ message(STATUS "Load custom task = ${TASK}")
 if (NOT ${TASK} STREQUAL "")
   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")
+    message(STATUS "Include ${CMAKE_CURRENT_SOURCE_DIR}/cmake/examples/${TASK}.cmake")
     include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/examples/${TASK}.cmake)
   else()
-    message(FATAL_ERROR "Cannot find custom file cmake/${TASK}.cmake")
+    message(FATAL_ERROR "Cannot find custom file cmake/examples/${TASK}.cmake")
   endif()
   set(TASK ${TASK} CACHE STRING "Custom task to be performed")
 else()
@@ -70,24 +78,17 @@ add_subdirectory(src)
 add_subdirectory(tools)
 add_subdirectory(data)
 
-## 5. set up python binding ##
-if (PYTHON_BINDINGS)
-  add_subdirectory(python)
-endif()
-
-## 6. set up unit tests and examples
+## 5. set up examples and tests
 message(STATUS "6. Set up unit tests")
 
-# load all macros
-FILE(GLOB _macro_files "${CMAKE_SOURCE_DIR}/cmake/macros/*.cmake")
-FOREACH(_file ${_macro_files})
-  MESSAGE(STATUS "Include ${_file}")
-  INCLUDE(${_file})
-ENDFOREACH()
-
 add_subdirectory(examples)
 add_subdirectory(tests)
 
+## 6. set up python binding ##
+if (PYTHON_BINDINGS)
+  add_subdirectory(python)
+endif()
+
 ## 7. add uninstall targets ##
 configure_file(
     "${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_uninstall.cmake.in"

cmake/compilers.cmake

@@ -35,10 +35,8 @@ endif()
 if(CMAKE_CXX_COMPILER_ID MATCHES "Intel")
   add_link_options("-fuse-ld=lld")
   set(CMAKE_INTERPROCEDURAL_OPTIMIZATION ON)
-  set(CMAKE_CXX_FLAGS_RELEASE
-      "-O3")
-  set(CMAKE_C_FLAGS_RELEASE
-      "-O3")
+  set(CMAKE_CXX_FLAGS_RELEASE "-O3")
+  set(CMAKE_C_FLAGS_RELEASE "-O3")
 
   set(CMAKE_CXX_FLAGS_DEBUG "-g3")
   set(CMAKE_C_FLAGS_DEBUG "-g3")

cmake/examples/juno.cmake

@@ -19,5 +19,6 @@ set_if_empty(NTRACER 2)
 # canoe configure
 set(HYDROSTATIC ON)
 set(NETCDF ON)
-set(FITS ON)
+# set(FITS ON)
 set(RadianceSolver lambert)
+set(PYTHON_BINDINGS ON)

cmake/examples/juno_mcmc.cmake

@@ -0,0 +1,24 @@
+# configuration for Jupiter Juno forward and inversion calculation
+
+macro(SET_IF_EMPTY _variable)
+  if("${${_variable}}" STREQUAL "")
+    set(${_variable} ${ARGN})
+  endif()
+endmacro()
+
+# athena variables
+set_if_empty(NUMBER_GHOST_CELLS 2)
+set_if_empty(NVAPOR 2)
+
+# canoe variables
+set_if_empty(NCLOUD 4)
+set_if_empty(NTRACER 2)
+
+# canoe task set(TASKLIST InversionTasks)
+
+# canoe configure
+set(HYDROSTATIC ON)
+set(NETCDF ON)
+set(FITS ON)
+set(RadianceSolver lambert)
+set(PYTHON_BINDINGS ON)

cmake/macros/macro_set_if_empty.cmake

@@ -0,0 +1,9 @@
+# A small macro used for set a variable if it is empty
+#
+# Usage: set_if_empty(name)
+
+macro(set_if_empty _variable)
+  if("${${_variable}}" STREQUAL "")
+    set(${_variable} ${ARGN})
+  endif()
+endmacro()

cmake/parameters.cmake

@@ -1,11 +1,5 @@
 # define default parameters
 
-macro(SET_IF_EMPTY _variable)
-  if("${${_variable}}" STREQUAL "")
-    set(${_variable} ${ARGN})
-  endif()
-endmacro()
-
 list(APPEND CMAKE_PREFIX_PATH "/usr/include/mpich-x84_64" "/usr/lib64/mpich/")
 
 # populate the default values

cmake/yamlpp.cmake

@@ -5,8 +5,7 @@ set(FETCHCONTENT_QUIET TRUE)
 FetchContent_Declare(
   yaml-cpp
   DOWNLOAD_EXTRACT_TIMESTAMP TRUE
-  URL https://github.com/jbeder/yaml-cpp/archive/refs/tags/0.8.0.tar.gz
-)
+  URL https://github.com/jbeder/yaml-cpp/archive/refs/tags/0.8.0.tar.gz)
 
 FetchContent_GetProperties(yaml-cpp)
 
@@ -16,7 +15,7 @@ if(NOT yaml-cpp_POPULATED)
   add_subdirectory(${yaml-cpp_SOURCE_DIR} ${yaml-cpp_BINARY_DIR})
 endif()
 
-#FetchContent_MakeAvailable(yamlcpp)
+# FetchContent_MakeAvailable(yamlcpp)
 
 # Where yaml-cpp's .h files can be found.
 include_directories(${yaml-cpp_SOURCE_DIR}/include)

examples/2019-Li-snap/bryan.cpp

@@ -36,7 +36,7 @@
 #include <climath/root.hpp>
 
 // snap
-#include <snap/thermodynamics/thermodynamics.hpp>
+#include <snap/thermodynamics/atm_thermodynamics.hpp>
 
 // special includes
 #include "bryan_vapor_functions.hpp"
@@ -120,13 +120,11 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
       air.c[iH2Oc] = 0.;
 
       // half a grid to cell center
-      pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2.,
-                           Thermodynamics::Method::ReversibleAdiabat, grav);
+      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),
-                             Thermodynamics::Method::ReversibleAdiabat, grav);
+        pthermo->Extrapolate(&air, pcoord->dx1f(i), "reversible", grav);
       }
     }
 
@@ -141,7 +139,8 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
         if (L < 1.) {
           auto &&air = AirParcelHelper::gather_from_conserved(this, k, j, i);
           air.ToMoleFraction();
-          Real temp_v = air.w[IDN] * pthermo->RovRd(air);
+          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;
@@ -153,8 +152,9 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
                                pthermo->TryEquilibriumTP_VaporCloud(air, iH2O);
                            air.w[iH2O] += rates[0];
                            air.c[iH2Oc] += rates[1];
+                           Real rovrd = get_rovrd(air, pthermo->GetMuRatio());
 
-                           return temp * pthermo->RovRd(air) - temp_v;
+                           return temp * rovrd - temp_v;
                          });
 
           if (err) throw RuntimeError("pgen", "TVSolver doesn't converge");

examples/2019-Li-snap/jupiter_crm.cpp

@@ -34,7 +34,7 @@
 #include <climath/interpolation.h>
 
 // snap
-#include <snap/thermodynamics/thermodynamics.hpp>
+#include <snap/thermodynamics/atm_thermodynamics.hpp>
 
 // special includes
 #include <special/giants_enroll_vapor_functions_v1.hpp>
@@ -87,7 +87,9 @@ void Forcing(MeshBlock *pmb, Real const time, Real const dt,
       for (int i = is; i <= ie; ++i) {
         auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i);
 
-        Real cv = pthermo->GetCvMass(air, 0);
+        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 *
@@ -159,8 +161,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
 
     // stop at just above P0
     for (int i = is; i <= ie; ++i) {
-      pthermo->Extrapolate(&air, pcoord->dx1f(i),
-                           Thermodynamics::Method::PseudoAdiabat, grav);
+      pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav);
       if (air.w[IPR] < P0) break;
     }
 
@@ -186,24 +187,20 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
       air.w[IDN] = Ts;
 
       // half a grid to cell center
-      pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2.,
-                           Thermodynamics::Method::ReversibleAdiabat, grav);
+      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),
-                             Thermodynamics::Method::PseudoAdiabat, grav,
-                             1.e-5);
+        pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5);
       }
 
       // 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),
-                             Thermodynamics::Method::Isothermal, grav);
+        pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav);
       }
     }
 }

examples/2024-CLi-earth-rt/run_ktable_earth.py

@@ -4,8 +4,8 @@
 sys.path.append("../python")
 sys.path.append(".")
 
-from pyharp import radiation_band
-from utilities import load_configure, find_resource
+from canoe.harp import radiation_band
+from canoe import load_configure, find_resource
 from numpy import *
 from rfmlib import *
 
@@ -44,7 +44,7 @@ def create_rfm_atmosphere(nlyr: int) -> dict:
 
     band = radiation_band(band_name, info)
 
-    nspec = band.get_num_spec_grids()
+    nspec = band.get_num_specgrids()
     wmin, wmax = band.get_range()
     wres = (wmax - wmin) / (nspec - 1)
 

examples/2023-jupiter-mwr-eq/CMakeLists.txt → examples/2024-CLi-juno-mwr/CMakeLists.txt

@@ -5,8 +5,13 @@
 # 1. Compile juno_mwr.cpp
 setup_problem(juno_mwr)
 
+# PGEN variable will be used in python bindings
+set(PGEN ${CMAKE_CURRENT_SOURCE_DIR}/juno_mwr
+  CACHE INTERNAL "Path to the source file for problem generator")
+
 # 4. Copy input file to run directory
-file(GLOB inputs *.inp *.yaml)
+file(GLOB inputs *.inp *.yaml *.py)
 foreach(input ${inputs})
-    file(COPY ${input} DESTINATION ${CMAKE_BINARY_DIR}/bin)
+  # softlink inp files
+  execute_process(COMMAND ln -sf ${input} ${CMAKE_BINARY_DIR}/bin/${inp})
 endforeach()