From 9a48fb136bb7c7df47ede175a77fd5531e86839b Mon Sep 17 00:00:00 2001 From: Xi Zhang Date: Thu, 11 Jul 2024 15:11:24 -0400 Subject: [PATCH] fix boiling --- src/microphysics/kessler94.cpp | 65 +++++++++++++++++++ .../thermodynamics/saturation_adjustment.cpp | 19 +++--- src/snap/thermodynamics/thermodynamics.hpp | 11 +++- 3 files changed, 84 insertions(+), 11 deletions(-) diff --git a/src/microphysics/kessler94.cpp b/src/microphysics/kessler94.cpp index 20c16141..93f26581 100644 --- a/src/microphysics/kessler94.cpp +++ b/src/microphysics/kessler94.cpp @@ -8,6 +8,8 @@ #include // microphysics +#include + #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: " <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= "<svp_func1_[iv][0](*air, iv, j) << std::endl; +//std::cout << "iv=:" << iv << " j= "<svp_func1_[iv][0](*air, iv, j) << std::endl; + +if (xs > 0.) { // test + //std::cout << "svp nonzero: " <getSVPFunc1(*air, iv, j, n)<<" +P:"<< air->w[IPR]<<" T:"<< tem << " Bsvp:" <svp_func1_[iv][n](*air, iv, j)<<" +P:"<< air->w[IPR]<<" T:"<< tem << " Bsvp:" <svp_func1_[iv][0](*air, iv, j) << std::endl; + //std::cout << "svp nonzero: " <svp_func1_[iv][n](*air, iv, n)<<" +P:"<< air->w[IPR]<<" T:"<< tem << " Bsvp:" <svp_func1_[iv][0](airmole, 0, +0)<<" P:"<< air->w[IPR]<<" T:"<< tem << " Bsvp:" <svp_func1_[1][1](*air, iv, n)<<" +P:"<< air->w[IPR]<<" T:"<< tem << " Bsvp:" < 1.) { // boiling + //std::cout << "boiling: " <getSVPFunc1(*air, iv, j, n)<< +std::endl; +// std::cout << "boiling: " <svp_func1_[iv][n](*air, iv, j)<< +std::endl; + + air->w[iv] += air->w[ip]; + air->w[ip] = 0.; + } +} } +*/ void Kessler94::SetVsedFromConserved(AthenaArray 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 index c330e9dc..a8cbed68 100644 --- 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 index 2251e4a9..158edd47 100644 --- 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 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 p3_; //! saturation vapor pressure function: Vapor -> Cloud - SVPFunc1Container svp_func1_; + // SVPFunc1Container svp_func1_; //! cloud index set - std::vector cloud_index_set_; + // std::vector cloud_index_set_; //! saturation vapor pressure function: Vapor + Vapor -> Cloud SVPFunc2Container svp_func2_;