Finish Jupiter equator model (#36)

Brewfile

@@ -14,3 +14,4 @@ brew "netcdf-cxx"
 brew "gfortran"
 brew "wget"
 brew "node"
+brew "cfitsio"

CMakeLists.txt

@@ -61,6 +61,7 @@ include(${CMAKE_SOURCE_DIR}/cmake/application.cmake)
 ## 4. set up project configure file and library ##
 message(STATUS "4. Set up project libraries")
 configure_file(${CMAKE_SOURCE_DIR}/configure.hpp.in configure.hpp @ONLY)
+configure_file(${CMAKE_SOURCE_DIR}/main.cpp.in main_${TASK}.cpp @ONLY)
 add_subdirectory(src)
 add_subdirectory(tools)
 add_subdirectory(data)

cmake/application.cmake

@@ -6,7 +6,7 @@ FetchContent_Declare(
   application
   # GIT_REPOSITORY https://github.com/chengcli/application/ GIT_TAG cli/flush)
   DOWNLOAD_EXTRACT_TIMESTAMP TRUE
-  URL https://github.com/chengcli/application/archive/refs/tags/v0.4.2.tar.gz)
+  URL https://github.com/chengcli/application/archive/refs/tags/v0.5.1.tar.gz)
 
 FetchContent_MakeAvailable(application)
 

cmake/athena.cmake

@@ -4,9 +4,10 @@ set(FETCHCONTENT_QUIET FALSE)
 
 FetchContent_Declare(
   athenapp
-  # GIT_REPOSITORY https://github.com/chengcli/athenapp/ GIT_TAG snap-mods)
-  DOWNLOAD_EXTRACT_TIMESTAMP TRUE
-  URL https://github.com/chengcli/athenapp/archive/refs/tags/v0.5.tar.gz)
+  GIT_REPOSITORY https://github.com/chengcli/athenapp/
+  GIT_TAG snap-mods)
+# DOWNLOAD_EXTRACT_TIMESTAMP TRUE URL
+# https://github.com/chengcli/athenapp/archive/refs/tags/v0.5.tar.gz)
 
 FetchContent_MakeAvailable(athenapp)
 

cmake/juno.cmake

@@ -1,19 +1,22 @@
 # configuration for Jupiter Juno forward and inversion calculation
 
-macro(set_if_empty _variable)
+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 3)
+set_if_empty(NVAPOR 2)
 
 # canoe variables
-set_if_empty(NCLOUD 5)
-set_if_empty(NTRACER 3)
-set_if_empty(NCHEMISTRY 0)
+set_if_empty(WATER_VAPOR_ID 1)
+set_if_empty(AMMONIA_VAPOR_ID 2)
+set_if_empty(NCLOUD 4)
+set_if_empty(NTRACER 2)
+
+# canoe task set(TASKLIST InversionTasks)
 
 # canoe configure
 set(HYDROSTATIC ON)

cmake/modules/Findcfitsio.cmake

@@ -0,0 +1,31 @@
+# * Try to find CFITSIO. Once executed, this module will define: Variables
+#   defined by this module: CFITSIO_FOUND        - system has CFITSIO
+#   CFITSIO_INCLUDE_DIR  - the CFITSIO include directory (cached)
+#   CFITSIO_INCLUDE_DIRS - the CFITSIO include directories (identical to
+#   CFITSIO_INCLUDE_DIR) CFITSIO_LIBRARY      - the CFITSIO library (cached)
+#   CFITSIO_LIBRARIES    - the CFITSIO libraries (identical to CFITSIO_LIBRARY)
+#
+# This module will use the following enviornmental variable when searching for
+# CFITSIO: CFITSIO_ROOT_DIR     - CFITSIO root directory
+#
+
+if(NOT cfitsio_FOUND)
+  find_path(
+    CFITSIO_INCLUDE_DIR fitsio.h
+    HINTS $ENV{CFITSIO_ROOT_DIR} /usr /usr/local/
+    PATH_SUFFIXES include include/cfitsio)
+  find_library(
+    CFITSIO_LIBRARY cfitsio /usr/lib /usr/lib64
+    HINTS $ENV{CFITSIO_ROOT_DIR}
+    PATH_SUFFIXES lib)
+
+  mark_as_advanced(CFITSIO_INCLUDE_DIR CFITSIO_LIBRARY)
+
+  include(FindPackageHandleStandardArgs)
+  find_package_handle_standard_args(cfitsio DEFAULT_MSG CFITSIO_LIBRARY
+                                    CFITSIO_INCLUDE_DIR)
+
+  set(CFITSIO_INCLUDE_DIRS ${CFITSIO_INCLUDE_DIR})
+  set(CFITSIO_LIBRARIES ${CFITSIO_LIBRARY})
+
+endif(NOT cfitsio_FOUND)

cmake/parameters.cmake

@@ -1,6 +1,6 @@
 # define default parameters
 
-macro(set_if_empty _variable)
+macro(SET_IF_EMPTY _variable)
   if("${${_variable}}" STREQUAL "")
     set(${_variable} ${ARGN})
   endif()
@@ -23,6 +23,8 @@ set_if_empty(WATER_VAPOR_ID -1)
 
 set_if_empty(EQUATION_OF_STATE "ideal_moist")
 
+set_if_empty(TASKLIST TimeIntegratorTaskList)
+
 set_if_empty(PLANET "Jupiter")
 
 if(NOT NETCDF OR NOT DEFINED NETCDF)
@@ -39,6 +41,13 @@ else()
   find_package(PNetCDF REQUIRED)
 endif()
 
+if(NOT FITS OR NOT DEFINED FITS)
+  set(FITS_OPTION "NO_FITSOUTPUT")
+else()
+  set(FITS_OPTION "FITSOUTPUT")
+  find_package(cfitsio REQUIRED)
+endif()
+
 if(NOT HYDROSTATIC OR NOT DEFINED HYDROSTATIC)
   set(HYDROSTATIC_OPTION "NOT_HYDROSTATIC")
 else()

configure.hpp.in

@@ -1,8 +1,5 @@
-#ifndef CONFIGURE_HPP
-#define CONFIGURE_HPP
-
-// C/C++
-#include <string>
+#ifndef CONFIGURE_HPP_
+#define CONFIGURE_HPP_
 
 enum {
   NCLOUD = @NCLOUD@,
@@ -26,7 +23,10 @@ enum class VariableType {prim, cons, chem};
 // PNetCDF output (PNETCDFOUTPUT or NO_PNETCDFOUTPUT)
 #define @PNETCDF_OPTION@
 
+// FITS option (FITSOUTPUT or NOT_FITSOUTPUT)
+#define @FITS_OPTION@
+
 // HYDROSTATIC option (HYDROSTATIC or NOT_HYDROSTATIC)
 #define @HYDROSTATIC_OPTION@
 
-#endif
+#endif  // CONFIGURE_HPP_

examples/2023-jupiter-mwr-eq/CMakeLists.txt

@@ -5,7 +5,8 @@ string(TOUPPER ${CMAKE_BUILD_TYPE} buildu)
 
 set(namel juno)
 
-add_executable(${namel}.${buildl} ${namel}.cpp)
+add_executable(${namel}.${buildl} ${namel}.cpp
+  ${CMAKE_BINARY_DIR}/main_${TASK}.cpp)
 
 set_target_properties(
     ${namel}.${buildl}

examples/2023-jupiter-mwr-eq/juno.cpp

@@ -15,17 +15,17 @@
 #include <athena/outputs/outputs.hpp>
 #include <athena/stride_iterator.hpp>
 
+// application
+#include <application/application.hpp>
+#include <application/exceptions.hpp>
+
 // canoe
 #include <configure.hpp>
 #include <impl.hpp>
 #include <constants.hpp>
 #include <variable.hpp>
 #include <index_map.hpp>
 
-// application
-#include <application/application.hpp>
-#include <application/exceptions.hpp>
-
 // climath
 #include <climath/interpolation.h>
 
@@ -53,9 +53,8 @@
 // inversion
 #include <inversion/profile_inversion.hpp>
 
-
-Real grav, P0, T0, xHe, xCH4, Tmin, clat;
-Real xH2S, xNa, xKCl, metallicity;
+Real grav, P0, T0, Tmin, clat;
+Real xHe, xCH4, xH2S, xNa, xKCl, metallicity;
 
 void MeshBlock::InitUserMeshBlockData(ParameterInput *pin)
 {
@@ -97,21 +96,29 @@ void Mesh::InitUserMeshData(ParameterInput *pin)
   Tmin = pin->GetReal("problem", "Tmin");
   clat = pin->GetOrAddReal("problem", "clat", 0.);
 
-  //xH2S = pin->GetReal("problem", "xH2S");
-  //metallicity = pin->GetOrAddReal("problem", "metallicity", 0.);
-  //xNa = pin->GetReal("problem", "xNa");
-  //xNa *= pow(10., metallicity);
-  //xKCl = pin->GetReal("problem", "xKCl");
-  //xKCl *= pow(10., metallicity);
+  xH2S = pin->GetReal("problem", "xH2S");
+
+  metallicity = pin->GetOrAddReal("problem", "metallicity", 0.);
+
+  xNa = pin->GetReal("problem", "xNa");
+  xNa *= pow(10., metallicity);
+
+  xKCl = pin->GetReal("problem", "xKCl");
+  xKCl *= pow(10., metallicity);
+
+  xHe = pin->GetReal("problem", "xKCl");
+
+  xCH4 = pin->GetReal("problem", "xKCl");
 }
 
 void update_gamma(Real* gamma, Real const q[]) {
   //std::cout << "I'm here" << std::endl;
   Real T = q[IDN], cp_h2, cp_he, cp_ch4;
-  if (T < 300.)
+  if (T < 300.) {
     cp_h2 = Hydrogen::cp_norm(T);
-  else
+  } else {
     cp_h2 = Hydrogen::cp_nist(T);
+  }
   cp_he = Helium::cp_nist(T);
   cp_ch4 = Methane::cp_nist(T);
 
@@ -124,6 +131,8 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
   Application::Logger app("main");
   app->Log("ProblemGenerator: juno");
 
+  auto pthermo = pimpl->pthermo;
+
   // mesh limits
   Real x1min = pmy_mesh->mesh_size.x1min;
   Real x1max = pmy_mesh->mesh_size.x1max;
@@ -135,7 +144,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
 
   // thermodynamic constants
   Real gamma = pin->GetReal("hydro", "gamma");
-  Real Rd = pin->GetReal("thermodynamics", "Rd");
+  Real Rd = pthermo->GetRd();
   Real cp = gamma/(gamma - 1.)*Rd;
 
   // estimate surface temperature and pressure
@@ -162,8 +171,6 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
   for (int i = 1; i < nx1; ++i)
     z1[i] = z1[i-1] + H0*dlnp;
 
-  auto pthermo = pimpl->pthermo;
-
   // set up an adiabatic atmosphere
   int max_iter = 200, iter = 0;
   Real t0;
@@ -190,12 +197,12 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
     Ts += T0 - t0;
     if (fabs(T0 - t0) < 0.01) break;
 
-    app->Log("iteration #", iter);
+    app->Log("Iteration #", iter);
     app->Log("T", t0);
   }
 
   if (iter > max_iter) {
-    throw InvalidValueError("ProblemGenerator", "maximum iteration reached");
+    throw RuntimeError("ProblemGenerator", "maximum iteration reached");
   }
 
   // change to log quantity
@@ -223,10 +230,9 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
           }
         }
 
-    // set tracers, electron, Na and K
+    // set tracers, electron and Na
     int ielec = pindex->GetTracerId("e-");
     int iNa = pindex->GetTracerId("Na");
-    int iK = pindex->GetTracerId("K");
     auto ptracer = pimpl->ptracer;
 
     for (int k = ks; k <= ke; ++k)
@@ -263,82 +269,26 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin)
   peos->PrimitiveToConserved(phydro->w, pfield->bcc,
     phydro->u, pcoord, is, ie, js, je, ks, ke);
 
-  FreeCArray(w1);
-  delete[] z1;
-  delete[] t1;
-}
-
-int main(int argc, char **argv)  {
-  Application::Logger app("main");
-
-  // no MPI
-  Globals::my_rank = 0;
-  Globals::nranks  = 1;
-
-  IOWrapper infile;
-  infile.Open("juno.inp", IOWrapper::FileMode::read);
-
-  ParameterInput *pinput = new ParameterInput;
-  pinput->LoadFromFile(infile);
-  infile.Close();
-
-  // set up mesh
-  int restart = false;
-  int mesh_only = false;
-
-  Mesh *pmesh = new Mesh(pinput, mesh_only);
-
-  // set up components
-  for (int b = 0; b < pmesh->nblocal; ++b) {
-    MeshBlock *pmb = pmesh->my_blocks(b);
-    pmb->pindex = std::make_shared<MeshBlock::IndexMap>(pmb, pinput);
-    pmb->pimpl = std::make_shared<MeshBlock::Impl>(pmb, pinput);
-  }
-
-  // initialize mesh
-  pmesh->Initialize(restart, pinput);
+  // Microwave radiative transfer needs temperatures at cell interfaces, which are
+  // interpolated from cell centered hydrodynamic variables.
+  // Normally, the boundary conditions are taken care of internally.
+  // But, since we call radiative tranfer directly in pgen, we would need to update the
+  // boundary conditions manually. The following lines of code updates the boundary
+  // conditions.
+  phydro->hbvar.SwapHydroQuantity(phydro->w, HydroBoundaryQuantity::prim);
+  pbval->ApplyPhysicalBoundaries(0., 0., pbval->bvars_main_int);
 
   // calculate radiative transfer
-  for (int b = 0; b < pmesh->nblocal; ++b) {
-    auto pmb = pmesh->my_blocks(b);
-    auto phydro = pmb->phydro;
-    auto prad = pmb->pimpl->prad;
-
-    int is = pmb->is, js = pmb->js, ks = pmb->ks;
-    int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
-
-    AthenaArray<Real> tempa, qNH3a, qH2Oa;
-    // read profile updates from input
-    std::string file_tempa = pinput->GetOrAddString("problem", "file.tempa", "");
-    if (file_tempa != "") {
-      ReadDataTable(&tempa, file_tempa);
-    }
-    std::string file_nh3a = pinput->GetOrAddString("problem", "file.nh3a", "");
-    if (file_nh3a != "") {
-      ReadDataTable(&qNH3a, file_nh3a);
-    }
-    std::string file_h2oa = pinput->GetOrAddString("problem", "file.h2oa", "");
-    if (file_h2oa != "") {
-      ReadDataTable(&qH2Oa, file_h2oa);
-    }
-
-    for (int k = ks; k <= ke; ++k) {
-      // run RT models
-      app->Log("run microwave radiative transfer");
-      for (int j = js-1; j <= je; ++j)
-        prad->CalRadiance(0., k, j, is, ie+1);
-    }
+  auto prad = pimpl->prad;
 
+  for (int k = ks; k <= ke; ++k) {
+    // run RT models
+    app->Log("run microwave radiative transfer");
+    for (int j = js; j <= je; ++j)
+      prad->CalRadiance(0., k, j, is, ie+1);
   }
 
-  // make output
-  Outputs *pouts;
-  pouts = new Outputs(pmesh, pinput);
-  pouts->MakeOutputs(pmesh, pinput);
-
-  delete pinput;
-  delete pmesh;
-  delete pouts;
-
-  Application::Destroy();
+  FreeCArray(w1);
+  delete[] z1;
+  delete[] t1;
 }