Fix boiling (#163)

Co-authored-by: Xi Zhang <[email protected]>
chengcli · Jul 12, 2024 · 5c04146 · 5c04146
1 parent e01dd77
commit 5c04146
Show file tree

Hide file tree

Showing 3 changed files with 84 additions and 11 deletions.
diff --git a/src/microphysics/kessler94.cpp b/src/microphysics/kessler94.cpp
@@ -8,6 +8,8 @@
 #include <snap/thermodynamics/thermodynamics.hpp>
 
 // microphysics
+#include <snap/thermodynamics/vapors/water_vapors.hpp>
+
 #include "microphysical_schemes.hpp"
 
 Kessler94::Kessler94(std::string name, YAML::Node const &node)
@@ -81,7 +83,70 @@ void Kessler94::EvolveOneStep(AirParcel *air, Real time, Real dt) {
 
   air->c[jbuf] += sol(0);
   for (int n = 1; n < Size; ++n) air->w[mySpeciesId(n)] += sol(n);
+
+  // boiling condition xiz 2024
+  //  get indices
+  int iv = mySpeciesId(0);
+  int ic = mySpeciesId(1);
+  int ip = mySpeciesId(2);
+
+  AirParcel airmole = *air;   // Copy the air object
+  airmole.ToMoleFraction();   // Modify the copy
+  Real tem = airmole.w[IDN];  // Extract the temperature or other property
+
+  Real xs = pthermo->svp_func1_[iv][0](airmole, iv, 0) / airmole.w[IPR];
+  if (xs > 1.) {  // boiling
+    //     std::cout << "boiling: " <<pthermo->svp_func1_[iv][0](airmole, iv, 0)
+    //     <<" P:"<< airmole.w[IPR] << std::endl;
+    air->w[iv] += air->w[ip];
+    air->w[ip] = 0.;
+  }
+}
+/*
+for (int n = 0; n < pthermo->cloud_index_set_[iv].size(); ++n) {
+//int j = pthermo->GetCloudIndex(iv,n);
+//  Real xs = pthermo->getSVPFunc1(*air, iv, j, n) / air->w[IPR];
+int j = pthermo->cloud_index_set_[iv][n];
+Real xs = pthermo->svp_func1_[iv][0](*air, iv, j) / air->w[IPR];
+
+//std::cout << "iv=:" << iv << " bpres=:" << pthermo->getSVPFunc1(*air, iv, 0)
+<< std::endl;
+//std::cout << "iv=:" << iv << "ic=:" << ic<< "ip=:" << ip<< " j= "<<j<<"
+bpres=:" << pthermo->svp_func1_[iv][0](*air, iv, j) << std::endl;
+//std::cout << "iv=:" << iv << " j= "<<j<<" bpres=:" <<
+pthermo->svp_func1_[iv][0](*air, iv, j) << std::endl;
+
+if (xs > 0.) { // test
+   //std::cout << "svp nonzero: " <<pthermo->getSVPFunc1(*air, iv, j, n)<<"
+P:"<< air->w[IPR]<<" T:"<< tem << " Bsvp:" <<sat_vapor_p_H2O_BriggsS(tem) <<
+std::endl;
+   //std::cout << "svp nonzero: " <<pthermo->svp_func1_[iv][n](*air, iv, j)<<"
+P:"<< air->w[IPR]<<" T:"<< tem << " Bsvp:" <<sat_vapor_p_H2O_BriggsS(tem) <<
+std::endl;
+//     std::cout << "iv=:" << iv << " j= "<<j<< " n= "<<n<<" bpres=:" <<
+pthermo->svp_func1_[iv][0](*air, iv, j) << std::endl;
+   //std::cout << "svp nonzero: " <<pthermo->svp_func1_[iv][n](*air, iv, n)<<"
+P:"<< air->w[IPR]<<" T:"<< tem << " Bsvp:" <<sat_vapor_p_H2O_BriggsS(tem) <<
+std::endl; std::cout << "svp nonzero: " <<pthermo->svp_func1_[iv][0](airmole, 0,
+0)<<" P:"<< air->w[IPR]<<" T:"<< tem << " Bsvp:" <<sat_vapor_p_H2O_BriggsS(tem)
+<< std::endl;
+   //std::cout << "svp nonzero: " <<pthermo->svp_func1_[1][1](*air, iv, n)<<"
+P:"<< air->w[IPR]<<" T:"<< tem << " Bsvp:" <<sat_vapor_p_H2O_BriggsS(tem) <<
+std::endl;
+}
+
+if (xs > 1.) { // boiling
+   //std::cout << "boiling: " <<pthermo->getSVPFunc1(*air, iv, j, n)<<
+std::endl;
+//     std::cout << "boiling: " <<pthermo->svp_func1_[iv][n](*air, iv, j)<<
+std::endl;
+
+   air->w[iv] += air->w[ip];
+   air->w[ip] = 0.;
+  }
+}
 }
+*/
 
 void Kessler94::SetVsedFromConserved(AthenaArray<Real> vsed[3],
                                      Hydro const *phydro, int kl, int ku,

diff --git a/src/snap/thermodynamics/saturation_adjustment.cpp b/src/snap/thermodynamics/saturation_adjustment.cpp
@@ -47,11 +47,12 @@ void Thermodynamics::SaturationAdjustment(AirColumn& air_column) const {
 
       for (int i = 1; i <= NVAPOR; ++i) {
         Real qv = air.w[i];
-        if (qv < 0) {
-          isNeg = true;
-          std::cout << "gas " << i << " is negative: " << qv << std::endl;
-          std::cout << "Temp is " << air.w[IDN] << std::endl;
-        }
+        //        if (qv < 0) {
+        //          isNeg = true;
+        //          std::cout << "gas " << i << " is negative: " << qv <<
+        //          std::endl; std::cout << "Temp is " << air.w[IDN] <<
+        //          std::endl;
+        //        }
 
 #pragma omp simd reduction(+ : qv)
         for (auto j : cloud_index_set_[i]) qv += air.c[j];
@@ -65,10 +66,10 @@ void Thermodynamics::SaturationAdjustment(AirColumn& air_column) const {
 
         // vapor condensation rate
         air.w[i] += rates[0];
-        if (isNeg && i == 1) {
-          std::cout << "new qv: " << air.w[i] << std::endl;
-          std::cout << "rate: " << rates[0] << std::endl;
-        }
+        //        if (isNeg && i == 1) {
+        //          std::cout << "new qv: " << air.w[i] << std::endl;
+        //          std::cout << "rate: " << rates[0] << std::endl;
+        //        }
         // cloud concentration rates
         for (int j = 1; j < rates.size(); ++j) {
           air.c[cloud_index_set_[i][j - 1]] += rates[j];

diff --git a/src/snap/thermodynamics/thermodynamics.hpp b/src/snap/thermodynamics/thermodynamics.hpp
@@ -167,6 +167,13 @@ class Thermodynamics {
     return cv_ratio_mass_[n] * cvd;
   }
 
+  // xiz add svp_func1_1 and cloud_index_set_ to public 2024
+  //! saturation vapor pressure function: Vapor -> Cloud
+  SVPFunc1Container svp_func1_;
+
+  //! cloud index set
+  std::vector<IndexSet> cloud_index_set_;
+
   //! Ratio of specific heat capacity [J/(kg K)] at constant pressure
   //! \return $c_{p,i}/c_{p,d}$
   Real GetCpRatioMass(int n) const { return cp_ratio_mass_[n]; }
@@ -438,10 +445,10 @@ class Thermodynamics {
   std::array<Real, 1 + NVAPOR> p3_;
 
   //! saturation vapor pressure function: Vapor -> Cloud
-  SVPFunc1Container svp_func1_;
+  //  SVPFunc1Container svp_func1_;
 
   //! cloud index set
-  std::vector<IndexSet> cloud_index_set_;
+  //  std::vector<IndexSet> cloud_index_set_;
 
   //! saturation vapor pressure function: Vapor + Vapor -> Cloud
   SVPFunc2Container svp_func2_;