merge amars convective model to canoe (#2)

- add four condensable species
- works for ~10 days without CO2
- CO2 breaks

cmake/amars.cmake

@@ -0,0 +1,21 @@
+# configure file for rad unit test case
+
+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 4)
+set(NCLOUD 8)
+# set(NPHASE_LEGACY 2)
+set(NETCDF ON)
+set(PNETCDF ON)
+set(MPI ON)
+set(RSOLVER lmars)
+set(DISORT ON)
+set(PYTHON_BINDINGS ON)

examples/2024-CMetz-amars/CMakeLists.txt

@@ -0,0 +1,11 @@
+
+# 1. Compile amars
+if (${NVAPOR} EQUAL 4 AND ${NCLOUD} EQUAL 8 AND ${NPHASE_LEGACY} EQUAL 3)
+  setup_problem(amars)
+endif()
+
+# 2. Copy input files to run directory
+file(GLOB inputs *.inp *.yaml)
+foreach(input ${inputs})
+    file(COPY ${input} DESTINATION ${CMAKE_BINARY_DIR}/bin)
+endforeach()

examples/2024-CMetz-amars/amars.cpp

@@ -0,0 +1,213 @@
+// 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/thermodynamics.hpp>
+
+// harp
+#include <harp/radiation.hpp>
+
+// special includes
+#include <special/amars_enroll_vapor_functions_v1.hpp>
+
+Real grav, P0, T0, Tmin;
+int iH2O, iCO2, iH2S, iSO2;
+
+void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
+  AllocateUserOutputVariables(4 + NVAPOR);
+  SetUserOutputVariableName(0, "temp");
+  SetUserOutputVariableName(1, "theta");
+  SetUserOutputVariableName(2, "thetav");
+  SetUserOutputVariableName(3, "mse");
+  for (int n = 1; n <= NVAPOR; ++n) {
+    std::string name = "rh" + std::to_string(n);
+    SetUserOutputVariableName(3 + n, name.c_str());
+  }
+}
+
+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);
+        for (int n = 1; n <= NVAPOR; ++n)
+          user_out_var(3 + n, k, j, i) =
+              pthermo->RelativeHumidity(this, n, 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);
+
+        Real cv = pthermo->GetCvMass(air, 0);
+
+        // 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));
+        // }
+
+        // if (air.w[IDN] < Tmin) {
+        //   u(IEN,k,j,i) += w(IDN,k,j,i)*cv*(Tmin - temp)/sponge_tau*dt;
+        // }
+      }
+}
+
+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");
+
+  // index
+  auto pindex = IndexMap::GetInstance();
+  iH2O = pindex->GetVaporId("H2O");
+  iCO2 = pindex->GetVaporId("CO2");
+  iH2S = pindex->GetVaporId("H2S");
+  iSO2 = pindex->GetVaporId("SO2");
+  EnrollUserExplicitSourceFunction(Forcing);
+}
+
+void MeshBlock::ProblemGenerator(ParameterInput *pin) {
+  srand(Globals::my_rank + time(0));
+
+  Application::Logger app("main");
+  app->Log("ProblemGenerator: amars_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 xH2O = pin->GetReal("problem", "qH2O.ppmv") / 1.E6;
+  Real xCO2 = pin->GetReal("problem", "qCO2.ppmv") / 1.E6;
+  Real xH2S = pin->GetReal("problem", "qH2S.ppmv") / 1.E6;
+  Real xSO2 = pin->GetReal("problem", "qSO2.ppmv") / 1.E6;
+
+  while (iter++ < max_iter) {
+    // read in vapors
+    air.w[iH2O] = xH2O;
+    air.w[iCO2] = xCO2;
+    air.w[iH2S] = xH2S;
+    air.w[iSO2] = xSO2;
+    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),
+                           Thermodynamics::Method::PseudoAdiabat, 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[iH2O] = xH2O;
+      air.w[iCO2] = xCO2;
+      air.w[iH2S] = xH2S;
+      air.w[iSO2] = xSO2;
+      air.w[IPR] = Ps;
+      air.w[IDN] = Ts;
+
+      // half a grid to cell center
+      pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2.,
+                           Thermodynamics::Method::PseudoAdiabat, grav);
+
+      int i = is;
+      for (; i <= ie; ++i) {
+        if (air.w[IDN] < Tmin) break;
+        air.w[IVX] = 0.1 * sin(2. * M_PI * rand() / RAND_MAX);
+        AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
+        pthermo->Extrapolate(&air, pcoord->dx1f(i),
+                             Thermodynamics::Method::PseudoAdiabat, grav,
+                             1.e-3);
+      }
+
+      // 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);
+      }
+
+      peos->ConservedToPrimitive(phydro->u, phydro->w, pfield->b, phydro->w,
+                                 pfield->bcc, pcoord, is, ie, j, j, k, k);
+
+      pimpl->prad->CalFlux(this, k, j, is, ie + 1);
+    }
+}