Add kessler microphysics (#179)

- cloud microphysics are inside equilibrate_uv
- use explicit integration for slow chem
- hard code sedimentation velocity

examples/2019-Li-snap/bryan.inp

@@ -57,7 +57,7 @@ grav_acc1   = -9.81
 gamma       = 1.4           # gamma = C_p/C_v
 
 <problem>
-thermodynamics_config = earth-precip.yaml
+thermodynamics_config = earth-moist.yaml
 
 dT    = 2.0
 xc    = 10.E3

examples/2019-Li-snap/earth-precip.yaml

@@ -51,4 +51,5 @@ reactions:
     rate-constant: {A: 0.01, b: 0}
 
   - equation: H2O(p) => H2O
-    rate-constant: {A: 0.01, b: 0}
+    type: evaporation
+    rate-constant: {A: 0.1, b: 0}

patches/21.new_blockdt.patch

@@ -1,19 +1,27 @@
 diff --git a/src/hydro/new_blockdt.cpp b/src/hydro/new_blockdt.cpp
-index 0dfe3b29..71e4f7aa 100644
+index c664de5d5..9286f0d70 100644
 --- a/src/hydro/new_blockdt.cpp
 +++ b/src/hydro/new_blockdt.cpp
 @@ -29,6 +29,10 @@
  #include "hydro.hpp"
  #include "hydro_diffusion/hydro_diffusion.hpp"
-
+ 
 +// snap injection
 +#include <impl.hpp>
 +#include <snap/implicit/implicit_solver.hpp>
 +
  // MPI/OpenMP header
  #ifdef MPI_PARALLEL
  #include <mpi.h>
-@@ -112,9 +116,21 @@ void Hydro::NewBlockTimeStep() {
+@@ -82,6 +86,7 @@ void Hydro::NewBlockTimeStep() {
+ #pragma ivdep
+         for (int i=is; i<=ie; ++i) {
+           wi[IDN] = w(IDN,k,j,i);
++          for (int n=1; n<IVX; ++n) wi[n] = w(n,k,j,i);
+           wi[IVX] = w(IVX,k,j,i);
+           wi[IVY] = w(IVY,k,j,i);
+           wi[IVZ] = w(IVZ,k,j,i);
+@@ -115,9 +120,21 @@ void Hydro::NewBlockTimeStep() {
                Real speed1 = std::max(cspeed, (std::abs(wi[IVX]) + cs));
                Real speed2 = std::max(cspeed, (std::abs(wi[IVY]) + cs));
                Real speed3 = std::max(cspeed, (std::abs(wi[IVZ]) + cs));

src/snap/riemann/hllc_transform.cpp

@@ -249,6 +249,26 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
       tflxi[n] = sl * tfl[n] + sr * tfr[n];
       pmicro->mass_flux[dir](n, k, j, i) = tflxi[n];
     }*/
+
+    // sedimentation flux
+    if (ivx == IVX) {
+      if (i == iu) return;
+      Real vsed_[NHYDRO];
+      for (int n = 0; n <= NVAPOR + NCLOUD; ++n) vsed_[n] = 0.;
+      for (int n = 1 + NVAPOR + NCLOUD; n <= NVAPOR + NCLOUD + NPRECIP; ++n) {
+        vsed_[n] = -10.;
+      }
+
+      Real hr = (wri[IPR] * (kappar + 1.) + KE_r) / wri[IDN];
+      for (int n = 1 + NVAPOR; n <= NVAPOR + NCLOUD + NPRECIP; ++n) {
+        Real rho = wri[IDN] * wri[n];
+        flx(n, k, j, i) += vsed_[n] * rho;
+        flx(ivx, k, j, i) += vsed_[n] * rho * wri[ivx];
+        flx(ivy, k, j, i) += vsed_[n] * rho * wri[ivy];
+        flx(ivz, k, j, i) += vsed_[n] * rho * wri[ivz];
+        flx(IEN, k, j, i) += vsed_[n] * rho * hr;
+      }
+    }
   }
 
 #if defined(AFFINE) || defined(CUBED_SPHERE)  // need of deprojection

src/snap/riemann/lmars.cpp

@@ -18,9 +18,6 @@
 // snap
 #include <snap/thermodynamics/thermodynamics.hpp>
 
-// microphysics
-#include <microphysics/microphysics.hpp>
-
 void Hydro::RiemannSolver(int const k, int const j, int const il, int const iu,
                           int const ivx, AthenaArray<Real> &wl,
                           AthenaArray<Real> &wr, AthenaArray<Real> &flx,
@@ -30,7 +27,6 @@ void Hydro::RiemannSolver(int const k, int const j, int const il, int const iu,
   int dir = ivx - IVX;
 
   auto pthermo = Thermodynamics::GetInstance();
-  auto pmicro = pmy_block->pimpl->pmicro;
   MeshBlock *pmb = pmy_block;
 
   Real rhobar, pbar, cbar, ubar, hl, hr;
@@ -96,13 +92,22 @@ void Hydro::RiemannSolver(int const k, int const j, int const il, int const iu,
     // the "total" density To get the denisty of the dry species, the "dry
     // mixing ratio" (rdl/rdr) is multiplied
     // FIXME: remove this?
-    /*for (int n = 0; n < NMASS; ++n) {
-      Real vfld = ubar + pmicro->vsedf[dir](n, k, j, i);
-      if (vfld > 0.) {
-        pmicro->mass_flux[dir](n, k, j, i) = vfld * wli[IDN] * rdl;
-      } else {
-        pmicro->mass_flux[dir](n, k, j, i) = vfld * wri[IDN] * rdr;
+    if (ivx == IVX) {
+      if (i == iu) return;
+      Real vsed_[NHYDRO];
+      for (int n = 0; n <= NVAPOR + NCLOUD; ++n) vsed_[n] = 0.;
+      for (int n = 1 + NVAPOR + NCLOUD; n <= NVAPOR + NCLOUD + NPRECIP; ++n) {
+        vsed_[n] = -10.;
+      }
+
+      for (int n = 1 + NVAPOR; n <= NVAPOR + NCLOUD + NPRECIP; ++n) {
+        Real rho = wri[IDN] * wri[n];
+        flx(n, k, j, i) += vsed_[n] * rho;
+        flx(ivx, k, j, i) += vsed_[n] * rho * wri[ivx];
+        flx(ivy, k, j, i) += vsed_[n] * rho * wri[ivy];
+        flx(ivz, k, j, i) += vsed_[n] * rho * wri[ivz];
+        flx(IEN, k, j, i) += vsed_[n] * rho * hr;
       }
-    }*/
+    }
   }
 }

src/snap/thermodynamics/equilibrate_tp.cpp

@@ -30,7 +30,7 @@ void Thermodynamics::EquilibrateTP(Real temp, Real pres) const {
     }
     if (max_abs_rate < 1.E-8) break;
 
-    kinetics_->getActivityConcentrations(xfrac.data());
+    thermo.getMoleFractions(xfrac.data());
 
     /*std::cout << "iter: " << iter << std::endl;
     for (size_t i = 0; i < Size; ++i) {

src/snap/thermodynamics/equilibrate_uv.cpp

@@ -9,8 +9,8 @@
 // snap
 #include "thermodynamics.hpp"
 
-void Thermodynamics::EquilibrateUV() const {
-  kinetics_->setQuantityConcentration();
+void Thermodynamics::EquilibrateUV(Real dt) const {
+  kinetics_->setQuantityConcentration(dt);
 
   Eigen::VectorXd rates(Size);
   Eigen::VectorXd conc(Size);
@@ -26,7 +26,7 @@ void Thermodynamics::EquilibrateUV() const {
     kinetics_->getNetProductionRates(rates.data());
 
     // get concentration
-    kinetics_->getActivityConcentrations(conc.data());
+    thermo.getConcentrations(conc.data());
     thermo.getIntEnergy_RT(intEng.data());
     thermo.getCv_R(cv.data());
 

src/snap/thermodynamics/thermodynamics.hpp

@@ -100,7 +100,7 @@ class Thermodynamics {
   void EquilibrateTP(Real temp, Real pres) const;
 
   //! Thermodnamic equilibrium at current UV
-  void EquilibrateUV() const;
+  void EquilibrateUV(Real dt = 0.) const;
 
   template <typename T>
   void SetMassFractions(StrideIterator<T *> w) const;

src/tasklist/implicit_hydro_tasks.cpp

@@ -309,7 +309,7 @@ TaskStatus ImplicitHydroTasks::UpdateAllConserved(MeshBlock *pmb, int stage) {
 
         pthermo->SetConserved(u.at(k, j, i), m.at(k, j, i));
         // pthermo->Evolve(pmb->pmy_mesh->time, pmb->pmy_mesh->dt);
-        pthermo->EquilibrateUV();
+        pthermo->EquilibrateUV(pmb->pmy_mesh->dt);
         pthermo->GetConserved(u.at(k, j, i), m.at(k, j, i));
 
         /*std::cout << "after: " << std::endl;