From 31114ce44b70f2a890a4e51ece6d41cfcbf6e9bb Mon Sep 17 00:00:00 2001 From: mac/cli Date: Mon, 8 Jul 2024 11:47:26 -0400 Subject: [PATCH 1/8] wip --- CMakeLists.txt | 3 +- .../eos_equilibrate_sp.cpp_} | 0 src/modules/eos_equilibrate_tp.cpp_ | 76 +++ .../eos_equilibrate_uv.cpp_} | 0 src/modules/eos_rk4_integrate_z.cpp_ | 110 ++++ src/modules/eos_saturation_adjustment.cpp_ | 27 + src/modules/eos_thermodynamics.cpp_ | 541 ++++++++++++++++++ src/modules/eos_thermodynamics.hpp_ | 463 +++++++++++++++ src/snap/eos/ideal_moist_hydro.cpp | 5 +- src/snap/thermodynamics/equilibrate_tp.cpp | 83 +-- src/snap/thermodynamics/rk4_integrate_z.cpp | 43 +- .../thermodynamics/saturation_adjustment.cpp | 76 ++- src/snap/thermodynamics/thermodynamics.cpp | 181 +++--- src/snap/thermodynamics/thermodynamics.hpp | 43 +- 14 files changed, 1411 insertions(+), 240 deletions(-) rename src/{snap/thermodynamics/equilibrate_sp.cpp => modules/eos_equilibrate_sp.cpp_} (100%) create mode 100644 src/modules/eos_equilibrate_tp.cpp_ rename src/{snap/thermodynamics/equilibrate_uv.cpp => modules/eos_equilibrate_uv.cpp_} (100%) create mode 100644 src/modules/eos_rk4_integrate_z.cpp_ create mode 100644 src/modules/eos_saturation_adjustment.cpp_ create mode 100644 src/modules/eos_thermodynamics.cpp_ create mode 100644 src/modules/eos_thermodynamics.hpp_ diff --git a/CMakeLists.txt b/CMakeLists.txt index a5a85e9e..40b61897 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -58,12 +58,11 @@ endif() message(STATUS "Include ${CMAKE_CURRENT_SOURCE_DIR}/cmake/parameters.cmake") include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/parameters.cmake) find_package(Eigen3 REQUIRED) -find_package(Cantera REQUIRED) find_package(Torch REQUIRED) +find_package(Cantera) ## 3. set up project system libraries ## message(STATUS "3. Set up system libraries") -#message(STATUS "Include ${CMAKE_SOURCE_DIR}/cmake/athena.cmake") if (NOT ${CANTERA_FOUND} STREQUAL "TRUE") include(${CMAKE_SOURCE_DIR}/cmake/yamlpp.cmake) endif() diff --git a/src/snap/thermodynamics/equilibrate_sp.cpp b/src/modules/eos_equilibrate_sp.cpp_ similarity index 100% rename from src/snap/thermodynamics/equilibrate_sp.cpp rename to src/modules/eos_equilibrate_sp.cpp_ diff --git a/src/modules/eos_equilibrate_tp.cpp_ b/src/modules/eos_equilibrate_tp.cpp_ new file mode 100644 index 00000000..b4f508b9 --- /dev/null +++ b/src/modules/eos_equilibrate_tp.cpp_ @@ -0,0 +1,76 @@ +// cantera +#include +#include +#include + +// canoe +#include + +// snap +#include "atm_thermodynamics.hpp" + +void Thermodynamics::EquilibrateTP(AirParcel* air) const { + air->ToMassFraction(); + Real pres = air->w[IPR]; + for (int n = IVX; n < IVX + NCLOUD; ++n) air->w[n] = air->c[n - IVX]; + auto& thermo = kinetics_->thermo(); + thermo.setMassFractionsPartial(&air->w[1]); + thermo.setDensity(air->w[IDN]); + thermo.setPressure(air->w[IPR]); + + EquilibrateTP(); + + for (int n = 0; n < NCLOUD; ++n) air->c[n] = air->w[IVX + n]; + air->w[IPR] = pres; + air->ToMoleFraction(); +} + +void Thermodynamics::EquilibrateTP() const { + std::static_pointer_cast(kinetics_) + ->setQuantityMoleFraction(); + + std::vector rates(kinetics_->nTotalSpecies()); + std::vector mfrac(kinetics_->nTotalSpecies()); + + auto& thermo = kinetics_->thermo(); + Real temp = thermo.temperature(); + Real pres = thermo.pressure(); + + int iter = 0, max_iter = 3; + + while (iter++ < max_iter) { + // std::cout << "Iteration " << iter << std::endl; + + // get mole fraction + kinetics_->getActivityConcentrations(mfrac.data()); + + /* print initial conditions + std::cout << "Mole fractions" << std::endl; + for (size_t i = 0; i < mfrac.size(); ++i) { + std::cout << "Mole fraction " << i << ": " << mfrac[i] << std::endl; + }*/ + + kinetics_->getNetProductionRates(rates.data()); + /*for (size_t i = 0; i < rates.size(); ++i) { + std::cout << "NET Production " << i << ": " << rates[i] << std::endl; + }*/ + + // update mole fraction + for (size_t i = 1; i < rates.size(); ++i) { + mfrac[i] += rates[i]; + } + + thermo.setMoleFractions(mfrac.data()); + thermo.setTemperature(temp); + thermo.setPressure(pres); + + /* print again + std::cout << "After Mole fractions" << std::endl; + for (size_t i = 0; i < mfrac.size(); ++i) { + std::cout << "Mole fraction " << i << ": " << mfrac[i] << std::endl; + } + std::cout << "T = " << thermo.temperature() << std::endl; + std::cout << "P = " << thermo.pressure() << std::endl; + std::cout << "D = " << thermo.density() << std::endl;*/ + } +} diff --git a/src/snap/thermodynamics/equilibrate_uv.cpp b/src/modules/eos_equilibrate_uv.cpp_ similarity index 100% rename from src/snap/thermodynamics/equilibrate_uv.cpp rename to src/modules/eos_equilibrate_uv.cpp_ diff --git a/src/modules/eos_rk4_integrate_z.cpp_ b/src/modules/eos_rk4_integrate_z.cpp_ new file mode 100644 index 00000000..4b398718 --- /dev/null +++ b/src/modules/eos_rk4_integrate_z.cpp_ @@ -0,0 +1,110 @@ +// C/C++ +#include +#include +#include + +// cantera +#include +#include + +// thermodynamics +#include "atm_thermodynamics.hpp" + +void rk4_integrate_z(AirParcel* air, Real dz, std::string method, Real grav, + Real adTdz) { + auto pthermo = Thermodynamics::GetInstance(); + auto& thermo = pthermo->Kinetics()->thermo(); + + auto const& mu_ratio = pthermo->GetMuRatio(); + auto const& cp_ratio_mole = pthermo->GetCpRatioMole(); + + Real step[] = {0.5, 0.5, 1.}; + Real dTdz[4], chi[4]; + Real latent[1 + NVAPOR]; + Real temp = air->w[IDN]; + + std::vector rates(pthermo->Kinetics()->nTotalSpecies()); + std::vector enthalpy(pthermo->Kinetics()->nTotalSpecies()); + + air->ToMassFraction(); + Real pres = air->w[IPR]; + for (int n = IVX; n < IVX + NCLOUD; ++n) air->w[n] = air->c[n - IVX]; + thermo.setMassFractionsPartial(&air->w[1]); + thermo.setDensity(air->w[IDN]); + thermo.setPressure(air->w[IPR]); + + thermo.getEnthalpy_RT(enthalpy.data()); + pthermo->Kinetics()->getNetProductionRates(rates.data()); + + for (int n = 0; n < NCLOUD; ++n) air->c[n] = air->w[IVX + n]; + air->w[IPR] = pres; + air->ToMoleFraction(); + + for (int i = 1; i <= NVAPOR; ++i) { + if (rates[i] > 0.) { + latent[i] = enthalpy[i] - enthalpy[i + NVAPOR]; + } else { + latent[i] = 0.; + } + } + + for (int rk = 0; rk < 4; ++rk) { + pthermo->EquilibrateTP(air); + + if (method != "reversible") { + for (int j = 0; j < NCLOUD; ++j) air->c[j] = 0; + } + + Real q_gas = 1., q_eps = 1.; + for (int i = 1; i <= NVAPOR; ++i) { + q_eps += air->w[i] * (mu_ratio[i] - 1.); + } + + for (int j = 0; j < NCLOUD; ++j) { + q_eps += air->c[j] * (mu_ratio[1 + NVAPOR + j] - 1.); + q_gas += -air->c[j]; + } + + Real g_ov_Rd = grav / pthermo->GetRd(); + Real R_ov_Rd = q_gas / q_eps; + + if (method == "reversible" || method == "pseudo") { + chi[rk] = cal_dlnT_dlnP(*air, cp_ratio_mole, latent); + } else if (method == "dry") { + for (int i = 1; i <= NVAPOR; ++i) latent[i] = 0; + chi[rk] = cal_dlnT_dlnP(*air, cp_ratio_mole, latent); + } else { // isothermal + chi[rk] = 0.; + } + + dTdz[rk] = -chi[rk] * g_ov_Rd / R_ov_Rd + adTdz; + chi[rk] = -R_ov_Rd / g_ov_Rd * dTdz[rk]; + + // integrate over dz + Real chi_avg; + if (rk < 3) { + air->w[IDN] = temp + dTdz[rk] * dz * step[rk]; + chi_avg = chi[rk]; + } else { + air->w[IDN] = + temp + + 1. / 6. * (dTdz[0] + 2. * dTdz[1] + 2. * dTdz[2] + dTdz[3]) * dz; + chi_avg = 1. / 6. * (chi[0] + 2. * chi[1] + 2. * chi[2] + chi[3]); + } + + if (!(air->w[IDN] > 0.)) air->w[IDN] = temp; + + if (fabs(air->w[IDN] - temp) > 0.01) { + air->w[IPR] = pres * pow(air->w[IDN] / temp, 1. / chi_avg); + } else { // isothermal limit + air->w[IPR] = + pres * exp(-2. * g_ov_Rd * dz / (R_ov_Rd * (air->w[IDN] + temp))); + } + } + + // recondensation + pthermo->EquilibrateTP(air); + if (method != "reversible") { + for (int j = 0; j < NCLOUD; ++j) air->c[j] = 0; + } +} diff --git a/src/modules/eos_saturation_adjustment.cpp_ b/src/modules/eos_saturation_adjustment.cpp_ new file mode 100644 index 00000000..8bca5780 --- /dev/null +++ b/src/modules/eos_saturation_adjustment.cpp_ @@ -0,0 +1,27 @@ +// C/C++ +#include +#include +#include + +// athena +#include +#include + +// application +#include + +// canoe +#include +#include + +// snap +#include "atm_thermodynamics.hpp" + +void Thermodynamics::SaturationAdjustment(AirColumn& air_column) const { + // return if there's no vapor + if (NVAPOR == 0) return; + + for (auto& air : air_column) { + EquilibrateUV(&air); + } +} diff --git a/src/modules/eos_thermodynamics.cpp_ b/src/modules/eos_thermodynamics.cpp_ new file mode 100644 index 00000000..a1416e3f --- /dev/null +++ b/src/modules/eos_thermodynamics.cpp_ @@ -0,0 +1,541 @@ +// C/C++ header +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // fill + +// cantera +#include +#include + +// athena +#include +#include +#include +#include +#include + +// application +#include +#include + +// canoe +#include +#include +#include +#include +#include + +// climath +#include + +// snap +#include "atm_thermodynamics.hpp" +#include "molecules.hpp" + +static std::mutex thermo_mutex; + +const std::string Thermodynamics::input_key = "thermodynamics_config"; + +Thermodynamics::~Thermodynamics() { + Application::Logger app("snap"); + app->Log("Destroy Thermodynamics"); +} + +Thermodynamics* Thermodynamics::fromYAMLInput(std::string const& fname) { + mythermo_ = new Thermodynamics(); + + auto thermo = Cantera::newThermo(fname, "gas"); + + if (thermo->nSpecies() != 1 + NVAPOR + NCLOUD) { + throw RuntimeError("Thermodynamics", + "Number of species does not match the input file"); + } + + auto& kinetics = mythermo_->kinetics_; + kinetics = Cantera::newKinetics({thermo}, fname); + + // --------- vapor + cloud thermo --------- + std::vector mu(thermo->nSpecies()); + std::vector cp_mole(thermo->nSpecies()); + + // g/mol + thermo->getMolecularWeights(mu.data()); + + // J/kmol/K + thermo->getPartialMolarCp(cp_mole.data()); + + mythermo_->Rd_ = Cantera::GasConstant / mu[0]; + mythermo_->gammad_ref_ = cp_mole[0] / (cp_mole[0] - Cantera::GasConstant); + + // ---------- dimensionless properties ---------- + + // cp ratios + for (size_t i = 0; i < cp_mole.size(); ++i) { + mythermo_->cp_ratio_mole_[i] = cp_mole[i] / cp_mole[0]; + } + + // molecular weight ratios + for (size_t i = 0; i < mu.size(); ++i) { + mythermo_->mu_ratio_[i] = mu[i] / mu[0]; + } + + // beta parameter (dimensionless internal energy) + for (size_t i = 0; i < thermo->nSpecies(); ++i) { + auto& tp = thermo->species(i)->thermo; + size_t n; + int type; + Real tlow, thigh, pref; + std::vector coeffs(tp->nCoeffs()); + tp->reportParameters(n, type, tlow, thigh, pref, coeffs.data()); + + Real t0 = coeffs[0]; + Real h0 = coeffs[1]; + Real s0 = coeffs[2]; + Real cp0 = coeffs[3]; + mythermo_->beta_[i] = h0 / (Cantera::GasConstant * t0); + } + + return mythermo_; +} + +Thermodynamics* Thermodynamics::fromLegacyInput(ParameterInput* pin) { + if (NCLOUD < (NPHASE - 1) * NVAPOR) { + throw RuntimeError( + "Thermodynamics", + "NCLOUD < (NPHASE-1)*NVAPOR is not supported for legacy input"); + } + + mythermo_ = new Thermodynamics(); + + // Read molecular weight ratios + read_thermo_property(mythermo_->mu_ratio_.data(), "eps", NPHASE, 1., pin); + + // Read cp ratios + read_thermo_property(mythermo_->cp_ratio_mass_.data(), "rcp", NPHASE, 1., + pin); + + // Read beta parameter + read_thermo_property(mythermo_->beta_.data(), "beta", NPHASE, 0., pin); + + // Read triple point temperature + read_thermo_property(mythermo_->t3_.data(), "Ttriple", 1, 0., pin); + + // Read triple point pressure + read_thermo_property(mythermo_->p3_.data(), "Ptriple", 1, 0., pin); + + mythermo_->svp_func1_.resize(1 + NVAPOR); + mythermo_->cloud_index_set_.resize(1 + NVAPOR); + + for (int i = 0; i <= NVAPOR; ++i) { + mythermo_->cloud_index_set_[i].resize(NPHASE - 1); + mythermo_->svp_func1_[i].resize(NPHASE - 1); + for (int j = 1; j < NPHASE; ++j) { + mythermo_->cloud_index_set_[i][j - 1] = (j - 1) * NVAPOR + (i - 1); + } + } + + mythermo_->Rd_ = pin->GetOrAddReal("thermodynamics", "Rd", 1.); + mythermo_->gammad_ref_ = pin->GetReal("hydro", "gamma"); + + return mythermo_; +} + +Thermodynamics const* Thermodynamics::GetInstance() { + // RAII + std::unique_lock lock(thermo_mutex); + + if (mythermo_ == nullptr) { + mythermo_ = new Thermodynamics(); + } + + return mythermo_; +} + +Thermodynamics const* Thermodynamics::InitFromYAMLInput( + std::string const& fname) { + if (mythermo_ != nullptr) { + throw RuntimeError("Thermodynamics", "Thermodynamics has been initialized"); + } + + Application::Logger app("snap"); + app->Log("Initialize Thermodynamics"); + + return fromYAMLInput(fname); +} + +Thermodynamics const* Thermodynamics::InitFromAthenaInput(ParameterInput* pin) { + if (mythermo_ != nullptr) { + throw RuntimeError("Thermodynamics", "Thermodynamics has been initialized"); + } + + Application::Logger app("snap"); + app->Log("Initialize Thermodynamics"); + + if (pin->DoesParameterExist("thermodynamics", input_key)) { + std::string filename = pin->GetString("thermodynamics", input_key); + mythermo_ = fromYAMLInput(filename); + } else { // legacy input + mythermo_ = fromLegacyInput(pin); + } + + // alias + auto& Rd = mythermo_->Rd_; + auto& gammad = mythermo_->gammad_ref_; + auto& mu_ratio = mythermo_->mu_ratio_; + auto& inv_mu_ratio = mythermo_->inv_mu_ratio_; + auto& mu = mythermo_->mu_; + auto& inv_mu = mythermo_->inv_mu_; + + auto& cp_ratio_mole = mythermo_->cp_ratio_mole_; + auto& cp_ratio_mass = mythermo_->cp_ratio_mass_; + auto& cv_ratio_mole = mythermo_->cv_ratio_mole_; + auto& cv_ratio_mass = mythermo_->cv_ratio_mass_; + + auto& latent_energy_mass = mythermo_->latent_energy_mass_; + auto& latent_energy_mole = mythermo_->latent_energy_mole_; + + auto& beta = mythermo_->beta_; + auto& delta = mythermo_->delta_; + auto& t3 = mythermo_->t3_; + auto& p3 = mythermo_->p3_; + auto& cloud_index_set = mythermo_->cloud_index_set_; + + // molecular weight + mu[0] = Constants::Rgas / Rd; + + // inverse mu ratio + for (int n = 0; n < mu_ratio.size(); ++n) { + inv_mu_ratio[n] = 1. / mu_ratio[n]; + mu[n] = mu[0] * mu_ratio[n]; + inv_mu[n] = 1. / mu[n]; + cp_ratio_mass[n] = cp_ratio_mole[n] * inv_mu_ratio[n]; + } + + /* calculate latent energy = $\beta\frac{R_d}{\epsilon}T^r$ + for (int n = 0; n <= NVAPOR; ++n) { + latent_energy_mass[n] = 0.; + latent_energy_mole[n] = 0.; + } + + for (int i = 1; i <= NVAPOR; ++i) { + for (int j = 0; j < NPHASE - 1; ++j) { + int n = cloud_index_set[i][j] + 1 + NVAPOR; + latent_energy_mass[n] = beta[n] * Rd / mu_ratio[n] * t3[i]; + latent_energy_mole[n] = latent_energy_mass[n] * mu[n]; + } + } + + // calculate delta = $(\sigma_j - \sigma_i)*\epsilon_i*\gamma/(\gamma - 1)$ + for (int n = 0; n <= NVAPOR; ++n) delta[n] = 0.; + for (int i = 1; i <= NVAPOR; ++i) { + for (int j = 0; j < cloud_index_set[i].size(); ++j) { + int n = cloud_index_set[i][j] + 1 + NVAPOR; + delta[n] = (cp_ratio_mass[n] - cp_ratio_mass[i]) * mu_ratio[i] / + (1. - 1. / gammad); + } + }*/ + + // finally, set up cv ratio + // calculate cv_ratio = $\sigma_i + (1. - \gamma)/\epsilon_i$ + for (int n = 0; n <= NVAPOR; ++n) { + cv_ratio_mass[n] = + gammad * cp_ratio_mass[n] + (1. - gammad) * inv_mu_ratio[n]; + cv_ratio_mole[n] = cv_ratio_mass[n] * mu_ratio[n]; + } + + for (int n = 1 + NVAPOR; n < Size; ++n) { + cv_ratio_mass[n] = gammad * cp_ratio_mass[n]; + cv_ratio_mole[n] = cv_ratio_mass[n] * mu_ratio[n]; + } + + // saturation adjustment parmaeters + mythermo_->sa_max_iter_ = + pin->GetOrAddReal("thermodynamics", "sa.max_iter", 10); + mythermo_->sa_relax_ = pin->GetOrAddReal("thermodynamics", "sa.relax", 0.8); + mythermo_->sa_ftol_ = pin->GetOrAddReal("thermodynamics", "sa.ftol", 1.e-4); + + return mythermo_; +} + +void Thermodynamics::Destroy() { + std::unique_lock lock(thermo_mutex); + + if (Thermodynamics::mythermo_ != nullptr) { + delete Thermodynamics::mythermo_; + Thermodynamics::mythermo_ = nullptr; + } +} + +Real Thermodynamics::GetPres(MeshBlock const* pmb, int k, int j, int i) const { + Real gm1 = pmb->peos->GetGamma() - 1; + auto& u = pmb->phydro->u; + + Real rho = 0., fsig = 0., feps = 0.; +#pragma omp simd reduction(+ : rho, fsig, feps) + for (int n = 0; n <= NVAPOR; ++n) { + rho += u(n, k, j, i); + fsig += u(n, k, j, i) * cv_ratio_mass_[n]; + feps += u(n, k, j, i) * inv_mu_ratio_[n]; + } + + // TODO(cli): not correct for Cubed sphere + Real KE = + 0.5 * + (sqr(u(IM1, k, j, i)) + sqr(u(IM2, k, j, i)) + sqr(u(IM3, k, j, i))) / + rho; + return gm1 * (u(IEN, k, j, i) - KE) * feps / fsig; +} + +// Eq.71 in Li2019 +Real Thermodynamics::GetChi(MeshBlock const* pmb, int k, int j, int i) const { + Real gammad = pmb->peos->GetGamma(); + auto& w = pmb->phydro->w; + + Real qsig = 1., feps = 1.; +#pragma omp simd reduction(+ : qsig, feps) + for (int n = 1; n <= NVAPOR; ++n) { + feps += w(n, k, j, i) * (inv_mu_ratio_[n] - 1.); + qsig += w(n, k, j, i) * (cp_ratio_mass_[n] - 1.); + } + + return (gammad - 1.) / gammad * feps / qsig; +} + +// Eq.63 in Li2019 +Real Thermodynamics::GetGamma(MeshBlock const* pmb, int k, int j, int i) const { + Real gammad = pmb->peos->GetGamma(); + + auto&& air = AirParcelHelper::gather_from_primitive(pmb, k, j, i); + + Real fsig = 1., feps = 1.; +#pragma omp simd reduction(+ : fsig, feps) + for (int n = 1; n <= NVAPOR; ++n) { + fsig += air.w[n] * (cv_ratio_mass_[n] - 1.); + feps += air.w[n] * (inv_mu_ratio_[n] - 1.); + } + + for (int n = 0; n < NCLOUD; ++n) { + fsig += air.c[n] * (cv_ratio_mass_[n + 1 + NVAPOR] - 1.); + feps += -air.c[n]; + } + + return 1. + (gammad - 1.) * feps / fsig; +} + +Real Thermodynamics::MoistStaticEnergy(MeshBlock const* pmb, Real gz, int k, + int j, int i) const { + auto&& air = AirParcelHelper::gather_from_primitive(pmb, k, j, i); + air.ToMassConcentration(); + + Real temp = GetTemp(pmb, k, j, i); + Real rho = 0., LE = 0., IE = 0.; + +#pragma omp simd reduction(+ : IE, rho) + for (int n = 0; n <= NVAPOR; ++n) { + IE += air.w[n] * GetCpMassRef(n) * temp; + rho += air.w[n]; + } + +#pragma omp simd reduction(+ : IE, LE, rho) + for (int n = 0; n < NCLOUD; ++n) { + IE += air.c[n] * GetCpMassRef(1 + NVAPOR + n) * temp; + LE += -latent_energy_mass_[1 + NVAPOR + n] * air.c[n]; + rho += air.c[n]; + } + + return (IE + LE) / rho + gz; +} + +// TODO(cli): check +Real Thermodynamics::GetCpMass(MeshBlock const* pmb, int k, int j, + int i) const { + Real gammad = pmb->peos->GetGamma(); + auto& w = pmb->phydro->w; + + Real qsig = 1.; +#pragma omp simd reduction(+ : qsig) + for (int n = 1; n <= NVAPOR; ++n) + qsig += w(n, k, j, i) * (cp_ratio_mass_[n] - 1.); + return gammad / (gammad - 1.) * Rd_ * qsig; +} + +// TODO(cli): check +Real Thermodynamics::GetCvMass(MeshBlock const* pmb, int k, int j, + int i) const { + Real gammad = pmb->peos->GetGamma(); + auto& w = pmb->phydro->w; + + Real qsig = 1.; +#pragma omp simd reduction(+ : qsig) + for (int n = 1; n <= NVAPOR; ++n) + qsig += w(n, k, j, i) * (cv_ratio_mass_[n] - 1.); + return 1. / (gammad - 1.) * Rd_ * qsig; +} + +Real Thermodynamics::GetEnthalpyMass(MeshBlock const* pmb, int k, int j, + int i) const { + Real gammad = pmb->peos->GetGamma(); + auto& w = pmb->phydro->w; + + Real qsig = 1.; +#pragma omp simd reduction(+ : qsig) + for (int n = 1; n <= NVAPOR; ++n) + qsig += w(n, k, j, i) * (cp_ratio_mass_[n] - 1.); + return gammad / (gammad - 1.) * Rd_ * qsig * w(IDN, k, j, i); +} + +Real Thermodynamics::GetMu(MeshBlock const* pmb, int k, int j, int i) const { + auto& w = pmb->phydro->w; + + Real feps = 1.; +#pragma omp simd reduction(+ : feps) + for (int n = 1; n <= NVAPOR; ++n) feps += w(n, k, j, i) * (mu_ratio_[n] - 1.); + return mu_[0] * feps; +} + +Real Thermodynamics::RelativeHumidity(MeshBlock const* pmb, int n, int k, int j, + int i) const { + auto&& air = AirParcelHelper::gather_from_primitive(pmb, k, j, i); + air.ToMoleFraction(); + return get_relative_humidity(air, n); +} + +void Thermodynamics::Extrapolate(AirParcel* qfrac, Real dzORdlnp, + std::string method, Real grav, + Real userp) const { + qfrac->ToMassFraction(); + + auto& thermo = kinetics_->thermo(); + + double pres = qfrac->w[IPR]; + for (int n = IVX; n < IVX + NCLOUD; ++n) qfrac->w[n] = qfrac->c[n - IVX]; + thermo.setMassFractionsPartial(&qfrac->w[1]); + thermo.setDensity(qfrac->w[IDN]); + thermo.setPressure(qfrac->w[IPR]); + + EquilibrateTP(); + + thermo.getMassFractions(&qfrac->w[0]); + + for (int n = 0; n < NCLOUD; ++n) { + qfrac->c[n] = qfrac->w[IVX + n]; + } + + qfrac->w[IPR] = pres; + qfrac->w[IDN] = thermo.density(); + qfrac->w[IVX] = 0.; + qfrac->w[IVX] = 0.; + qfrac->w[IVX] = 0.; + + qfrac->ToMoleFraction(); + + // std::cout << *qfrac << std::endl; + + // RK4 integration +#ifdef HYDROSTATIC + rk4_integrate_lnp(qfrac, dzORdlnp, method, userp); +#else + rk4_integrate_z(qfrac, dzORdlnp, method, grav, userp); +#endif +} + +Real Thermodynamics::GetLatentHeatMole(int i, std::vector const& rates, + Real temp) const { + if (std::abs(rates[0]) < 1.E-8) return 0.; + + Real heat = 0.; + for (int j = 1; j < rates.size(); ++j) { + int n = cloud_index_set_[i][j - 1] + 1 + NVAPOR; + heat += rates[j] * GetLatentEnergyMole(n, temp); + } + + return heat / std::abs(rates[0]); +} + +void Thermodynamics::SetStateFromPrimitive(MeshBlock* pmb, int k, int j, + int i) const { + Hydro* phydro = pmb->phydro; + + auto& thermo = kinetics_->thermo(); + size_t total_cells = pmb->ncells1 * pmb->ncells2 * pmb->ncells3; + + thermo.setMassFractionsPartial(&phydro->w(1, k, j, i), total_cells); + thermo.setDensity(phydro->w(IDN, k, j, i)); + thermo.setPressure(phydro->w(IPR, k, j, i)); +} + +void Thermodynamics::SetStateFromConserved(MeshBlock* pmb, int k, int j, + int i) const { + Hydro* phydro = pmb->phydro; + + auto& thermo = kinetics_->thermo(); + size_t total_cells = pmb->ncells1 * pmb->ncells2 * pmb->ncells3; + + thermo.setMassFractions(&phydro->u(0, k, j, i), total_cells); + Real rho = 0.; + for (int n = 0; n < Size; ++n) { + rho += phydro->u(n, k, j, i); + } + thermo.setDensity(rho); + thermo.setPressure(GetPres(pmb, k, j, i)); +} + +// Eq.4.5.11 in Emanuel (1994) +Real Thermodynamics::EquivalentPotentialTemp(MeshBlock* pmb, Real p0, int v, + int k, int j, int i) const { +#if (NVAPOR > 0) + AirParcel&& air_mass = AirParcelHelper::gather_from_primitive(pmb, k, j, i); + + AirParcel air_mole = air_mass; + air_mole.ToMoleFraction(); + + // get dry air mixing ratio + Real xg = 1., qd = 1.; +#pragma omp simd reduction(+ : xg, qd) + for (int n = 0; n < NCLOUD; ++n) { + xg += -air_mole.c[n]; + qd += -air_mass.c[n]; + } + + Real xd = xg; +#pragma omp simd reduction(+ : xd, qd) + for (int n = 1; n <= NVAPOR; ++n) { + xd += -air_mole.w[n]; + qd += -air_mass.w[n]; + } + + Real temp = air_mole.w[IDN]; + Real pres = air_mole.w[IPR]; + + Real qv = air_mass.w[v]; + Real qc = air_mass.c[cloud_index_set_[v][0]]; + Real qt = qv + qc; + + Real Rd = Rd_; + Real Rv = Rd_ / mu_ratio_[v]; + + Real cpd = GetCpMassRef(0); + Real cl = GetCpMassRef(cloud_index_set_[v][0] + 1 + NVAPOR); + + std::vector rates{-0.01, 0.01}; + Real lv = GetLatentHeatMass(v, rates, temp); + + Real rh = RelativeHumidity(pmb, v, k, j, i); + Real pd = xd / xg * pres; + Real cpt = cpd * qd + cl * qt; + + return temp * pow(p0 / pd, Rd * qd / cpt) * pow(rh, -Rv * qv / cpt) * + exp(lv * qv / (cpt * temp)); +#else + return PotentialTemp(pmb, p0, k, j, i); +#endif +} + +Thermodynamics* Thermodynamics::mythermo_ = nullptr; diff --git a/src/modules/eos_thermodynamics.hpp_ b/src/modules/eos_thermodynamics.hpp_ new file mode 100644 index 00000000..33109270 --- /dev/null +++ b/src/modules/eos_thermodynamics.hpp_ @@ -0,0 +1,463 @@ +#ifndef SRC_SNAP_THERMODYNAMICS_THERMODYNAMICS_HPP_ +#define SRC_SNAP_THERMODYNAMICS_THERMODYNAMICS_HPP_ + +// C/C++ +#include +#include +#include +#include +#include +#include +#include +#include + +// athena +#include +#include +#include + +// canoe +#include +#include +#include + +class MeshBlock; +class ParameterInput; + +namespace Cantera { +class ThermoPhase; +class Kinetics; +} // namespace Cantera + +using IndexPair = std::pair; +using IndexSet = std::vector; + +using RealArray3 = std::array; +using RealArrayX = std::vector; + +using SatVaporPresFunc1 = Real (*)(AirParcel const &, int i, int j); +using SatVaporPresFunc2 = Real (*)(AirParcel const &, int i, int j, int k); + +enum { MAX_REACTANT = 3 }; + +using ReactionIndx = std::array; +using ReactionStoi = std::array; +using ReactionInfo = std::pair; + +Real NullSatVaporPres1(AirParcel const &, int, int); +Real NullSatVaporPres2(AirParcel const &, int, int, int); + +void read_thermo_property(Real var[], char const name[], int len, Real v0, + ParameterInput *pin); +Real saha_ionization_electron_density(Real T, Real num, Real ion_ev); + +//! Ideal saturation vapor pressure +//! $p^* = p^r\exp[\beta(1-1/t)-\delta\lnt]$ +//! $p^r$ is the reference pressure, usually choosen to be the triple point +//! pressure \n $t=T/T^r$ is the dimensionless temperature. $T^r$ is the +//! reference temperature. \n Similar to $p^r$, $T^r$ is usually choosen to be +//! the triple point temperature +//! \return $p^*$ [pa] +inline Real SatVaporPresIdeal(Real t, Real p3, Real beta, Real delta) { + return p3 * exp((1. - 1. / t) * beta - delta * log(t)); +} + +// Thermodynamic variables are ordered in the array as the following +// 0: dry air (non-condensible gas) +// 1..NVAPOR: moist air (condensible gas) +// NVAPOR+1..NVAPOR+NCLOUD: clouds + +// 0 is a special buffer place for cloud in equilibrium with vapor at the same +// temperature + +class Thermodynamics { + protected: + enum { NPHASE = NPHASE_LEGACY }; + + //! Constructor for class sets up the initial conditions + //! Protected ctor access thru static member function Instance + Thermodynamics() {} + static Thermodynamics *fromLegacyInput(ParameterInput *pin); + static Thermodynamics *fromYAMLInput(std::string const &fname); + + public: + using SVPFunc1Container = std::vector>; + using SVPFunc2Container = std::map; + + enum { Size = 1 + NVAPOR + NCLOUD }; + + //! thermodynamics input key in the input file [thermodynamics_config] + static const std::string input_key; + + // member functions + ~Thermodynamics(); + + //! Return a pointer to the one and only instance of Thermodynamics + static Thermodynamics const *GetInstance(); + + static Thermodynamics const *InitFromYAMLInput(std::string const &fname); + + static Thermodynamics const *InitFromAthenaInput(ParameterInput *pin); + + //! Destroy the one and only instance of Thermodynamics + static void Destroy(); + + //! Ideal gas constant of dry air in [J/(kg K)] + //! \return $R_d=\hat{R}/\mu_d$ + Real GetRd() const { return Rd_; } + + //! reference adiabatic index of dry air [1] + Real GetGammadRef() const { return gammad_ref_; } + + //! Molecular weight [kg/mol] + //! \param[in] n the index of the thermodynamic species + //! \return $\mu$ + Real GetMu(int n) const { return mu_[n]; } + + //! Inverse molecular weight [mol/kg] + //! \param[in] n the index of the thermodynamic species + //! \return $\mu$ + Real GetInvMu(int n) const { return inv_mu_[n]; } + + //! Ratio of molecular weights [1] + //! \param[in] n the index of the thermodynamic species + //! \return $\epsilon_i=\mu_i/\mu_d$ + Real GetMuRatio(int n) const { return mu_ratio_[n]; } + + //! const pointer to mu_ratio_ + Real const *GetMuRatio() const { return &mu_ratio_[0]; } + + //! Ratio of molecular weights [1] + //! \param[in] n the index of the thermodynamic species + //! \return $\epsilon_i^{-1} =\mu_d/\mu_i$ + Real GetInvMuRatio(int n) const { return inv_mu_ratio_[n]; } + + //! Ratio of specific heat capacity [J/(kg K)] at constant volume [1] + //! \param[in] n the index of the thermodynamic species + //! \return $c_{v,i}/c_{v,d}$ + Real GetCvRatioMass(int n) const { return cv_ratio_mass_[n]; } + + //! const pointer to cv_ratio_mass_ + Real const *GetCvRatioMass() const { return cv_ratio_mass_.data(); } + + //! Ratio of specific heat capacity [J/(mol K)] at constant volume [1] + //! \param[in] n the index of the thermodynamic species + //! \return $\hat{c}_{v,i}/\hat{c}_{v,d}$ + Real GetCvRatioMole(int n) const { return cv_ratio_mole_[n]; } + + //! const pointer to cv_ratio_mole_ + Real const *GetCvRatioMole() const { return cv_ratio_mole_.data(); } + + std::shared_ptr Kinetics() const { return kinetics_; } + + //! \return the index of the cloud + //! \param[in] i the index of the vapor + //! \param[in] j the sequential index of the cloud + int GetCloudIndex(int i, int j) const { return cloud_index_set_[i][j]; } + + //! \return the index set of the cloud + //! \param[in] i the index of the vapor + std::vector GetCloudIndexSet(int i) const { return cloud_index_set_[i]; } + + //! const pointer to cloud_index_set_ + IndexSet const *GetCloudIndexSet() const { return cloud_index_set_.data(); } + + //! Reference specific heat capacity [J/(kg K)] at constant volume + Real GetCvMassRef(int n) const { + Real cvd = Rd_ / (gammad_ref_ - 1.); + return cv_ratio_mass_[n] * cvd; + } + + //! 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]; } + + //! const pointer to cp_ratio_mass_ + Real const *GetCpRatioMass() const { return cp_ratio_mass_.data(); } + + //! Ratio of specific heat capacity [J/(mol K)] at constant pressure + //! \return $\hat{c}_{p,i}/\hat{c}_{p,d}$ + Real GetCpRatioMole(int n) const { return cp_ratio_mole_[n]; } + + //! const pointer to cp_ratio_mole_ + Real const *GetCpRatioMole() const { return cp_ratio_mole_.data(); } + + Real GetCpMassRef(int n) const { + Real cpd = Rd_ * gammad_ref_ / (gammad_ref_ - 1.); + return cp_ratio_mass_[n] * cpd; + } + + //! \brief Temperature dependent specific latent energy [J/kg] of condensates + //! at constant volume $L_{ij}(T) = L_{ij}^r - (c_{ij} - c_{p,i})\times(T - + //! T^r)$ + //! + //! $= L_{ij}^r - \delta_{ij}R_i(T - T^r)$ + //! \return $L_{ij}(T)$ + Real GetLatentEnergyMass(int n, Real temp = 0.) const { + return latent_energy_mass_[n] - delta_[n] * Rd_ * inv_mu_ratio_[n] * temp; + } + + //! \brief Temperature dependent specific latent energy [J/mol] of condensates + //! at constant volume + //! + //! \return $L_{ij}(T)$ + Real GetLatentEnergyMole(int n, Real temp = 0.) const { + return GetLatentEnergyMass(n, temp) * mu_[n]; + } + + //! const pointer to latent_energy_mole_ + Real const *GetLatentEnergyMole() const { return &latent_energy_mole_[0]; } + + Real GetLatentHeatMole(int i, std::vector const &rates, + Real temp) const; + + Real GetLatentHeatMass(int i, std::vector const &rates, + Real temp) const { + return GetLatentHeatMole(i, rates, temp) * inv_mu_[i]; + } + + //! \brief Calculate the equilibrium mole transfer by cloud reaction + //! vapor -> cloud + //! + //! \param[in] qfrac mole fraction representation of air parcel + //! \param[in] ivapor the index of the vapor + //! \param[in] cv_hat $cv_hat$ molar heat capacity + //! \param[in] misty if true, there is an infinite supple of cloud + //! \return molar fraction change of vapor to cloud + RealArrayX TryEquilibriumTP_VaporCloud(AirParcel const &qfrac, int ivapor, + Real cv_hat = 0., + bool misty = false) const; + + //! \brief Calculate the equilibrium mole transfer by cloud reaction + //! vapor1 + vapor2 -> cloud + //! + //! \param[in] air mole fraction representation of air parcel + //! \param[in] ij the index pair of the vapor1 and vapor2 + //! \param[in] cv_hat $\har{cv}$ molar heat capacity + //! \param[in] misty if true, there is an infinite supple of cloud + //! \return molar fraction change of vapor1, vapor2 and cloud + RealArray3 TryEquilibriumTP_VaporVaporCloud(AirParcel const &air, + IndexPair ij, Real cv_hat = 0., + bool misty = false) const; + + //! Construct an 1d atmosphere + //! \param[in,out] qfrac mole fraction representation of air parcel + //! \param[in] dzORdlnp vertical grid spacing + //! \param[in] method choose from [reversible, pseudo, dry, isothermal] + //! \param[in] grav gravitational acceleration + //! \param[in] userp user parameter to adjust the temperature gradient + void Extrapolate(AirParcel *qfrac, Real dzORdlnp, std::string method, + Real grav = 0., Real userp = 0.) const; + + //! Thermodnamic equilibrium at current TP + //! \param[in,out] qfrac mole fraction representation of air parcel + void EquilibrateTP() const; + + void EquilibrateTP(AirParcel *qfrac) const; + + //! Thermodnamic equilibrium at current UV + void EquilibrateUV() const; + + void EquilibrateUV(AirParcel *qfrac) const; + + void EquilibrateSP(double P) const; + + //! Adjust to the maximum saturation state conserving internal energy + //! \param[in,out] ac mole fraction representation of a collection of air + //! parcels + void SaturationAdjustment(AirColumn &ac) const; + + void SetStateFromPrimitive(MeshBlock *pmb, int k, int j, int i) const; + void SetStateFromConserved(MeshBlock *pmb, int k, int j, int i) const; + + //! \brief Calculate potential temperature from primitive variable + //! + //! $\theta = T(\frac{p_0}{p})^{\chi}$ + //! \return $\theta$ + Real PotentialTemp(MeshBlock *pmb, Real p0, int k, int j, int i) const { + auto &w = pmb->phydro->w; + return GetTemp(pmb, k, j, i) * + pow(p0 / w(IPR, k, j, i), GetChi(pmb, k, j, i)); + } + + //! \brief Calculate equivalent potential temperature from primitive variable + //! + //! $\theta_e = T(\frac{p}{p_d})^{Rd/(cpd + cl r_t} \exp(\frac{L_v q_v}{c_p + //! T})$ + Real EquivalentPotentialTemp(MeshBlock *pmb, Real p0, int v, int k, int j, + int i) const; + + //! \brief Calculate temperature from primitive variable + //! + //! $T = p/(\rho R) = p/(\rho \frac{R}{R_d} Rd)$ + //! \return $T$ + Real GetTemp(MeshBlock const *pmb, int k, int j, int i) const { + auto &w = pmb->phydro->w; + return w(IPR, k, j, i) / (w(IDN, k, j, i) * Rd_ * RovRd(pmb, k, j, i)); + } + + //! \brief Mean molecular weight + //! + //! $mu = \mu_d (1 + \sum_i q_i (\epsilon_i - 1))$ + //! \return $mu$ + Real GetMu(MeshBlock const *pmb, int k, int j, int i) const; + + //! \brief Inverse of the mean molecular weight (no cloud) + //! + //! $ \frac{R}{R_d} = \frac{\mu_d}{\mu}$ + //! \return $1/\mu$ + //! Eq.16 in Li2019 + Real RovRd(MeshBlock const *pmb, int k, int j, int i) const { + Real feps = 1.; + auto &w = pmb->phydro->w; +#pragma omp simd reduction(+ : feps) + for (int n = 1; n <= NVAPOR; ++n) + feps += w(n, k, j, i) * (inv_mu_ratio_[n] - 1.); + return feps; + } + + //! \brief Specific heat capacity [J/(kg K)] of the air parcel at constant + //! volume + //! + //! $c_v = c_{v,d}*(1 + \sum_i (q_i*(\hat{c}_{v,i} - 1.))) + //! = \gamma_d/(\gamma_d - 1.)*R_d*T*(1 + \sum_i (q_i*(\hat{c}_{v,i} + //! - 1.)))$ + //! \return $c_v$ + Real GetCvMass(MeshBlock const *pmb, int k, int j, int i) const; + + //! \brief Specific heat capacity [J/(kg K)] of the air parcel at constant + //! pressure + //! + //! $c_p = c_{p,d}*(1 + \sum_i (q_i*(\hat{c}_{p,i} - 1.))) + //! = \gamma_d/(\gamma_d - 1.)*R_d*T*(1 + \sum_i (q_i*(\hat{c}_{p,i} + //! - 1.)))$ + //! \return $c_p$ + Real GetCpMass(MeshBlock const *pmb, int k, int j, int i) const; + + //! \brief Adiabatic index + //! + //! $\chi = \frac{R}{c_p}$ + //! \return $\chi$ + Real GetChi(MeshBlock const *pmb, int k, int j, int i) const; + + //! \brief Polytropic index + //! + //! $\gamma = \frac{c_p}{c_v}$ + //! \return $\gamma$ + Real GetGamma(MeshBlock const *pmb, int k, int j, int i) const; + + //! Pressure from conserved variables + //! \return $p$ + Real GetPres(MeshBlock const *pmb, int k, int j, int i) const; + + //! \brief specific enthalpy [J/kg] of the air parcel + //! + //! $h = c_{p,d}*T*(1 + \sum_i (q_i*(\hat{c}_{p,i} - 1.)))$ + //! $ = \gamma_d/(\gamma_d - 1.)*R_d*T*(1 + \sum_i (q_i*(\hat{c}_{pi} + //! - 1.)))$ + Real GetEnthalpyMass(MeshBlock const *pmb, int k, int j, int i) const; + + //! \brief Moist static energy + //! + //! $h_s = c_{pd}T + gz + L_vq_v + L_s\sum_i q_i$ + //! \return $h_s$ + Real MoistStaticEnergy(MeshBlock const *pmb, Real gz, int k, int j, + int i) const; + + //! \brief Relative humidity + Real RelativeHumidity(MeshBlock const *pmb, int n, int k, int j, int i) const; + + protected: + //! custom functions for vapor (to be overridden by user mods) + void enrollVaporFunctions(); + + protected: + std::shared_ptr kinetics_; + + //! ideal gas constant of dry air in J/kg + Real Rd_; + + //! reference polytropic index of dry air + Real gammad_ref_; + + //! ratio of mean molecular weights + std::array mu_ratio_; + + //! inverse ratio of mean molecular weights + std::array inv_mu_ratio_; + + //! molecular weight [kg/mol] + std::array mu_; + + //! inverse molecular weight [mol/kg] + std::array inv_mu_; + + //! ratio of specific heat capacities [J/mol] at constant pressure + std::array cp_ratio_mole_; + + //! ratio of specific heat capacities [J/mol] at constant pressure + std::array cp_ratio_mass_; + + //! ratio of specific heat capacities [J/mol] at constant volume + std::array cv_ratio_mole_; + + //! ratio of specific heat capacities [J/mol] at constant volume + std::array cv_ratio_mass_; + + //! latent energy at constant volume [J/mol] + std::array latent_energy_mole_; + + //! \brief latent energy at constant volume [J/kg] + //! + //! $L_i^r+\Delta c_{ij}T^r == \beta_{ij}\frac{R_d}{\epsilon_i}T^r$ + std::array latent_energy_mass_; + + //! \brief Dimensionless latent heat + //! + //! $\beta_{ij} = \frac{\Delta U_{ij}}{R_i T^r}$ + //! $\Delta U_{ij}$ is the difference in internal energy between the vapor $i$ + //! and the condensate $j$ $R_i=\hat{R}/m_i=R_d/\epsilon_i$ $T^r$ is the + //! triple point temperature + //! $\beta_{ij} == \frac{L_{ij}^r + \Delta c_{ij}T^r}{R_i T^r}$ + std::array beta_; + + //! \brief Dimensionless differences in specific heat capacity + //! + //! $\delta_{ij} = \frac{c_{ij} - c_{p,i}}{R_i}$ + //! $c_{ij} - c_{p,j}$ is the difference in specific heat + //! capacity at constant pressure. $c_{ij}$ is the heat capacity of + //! condensate $j$ of vapor $i$ + std::array delta_; + + //! triple point temperature [K] + std::array t3_; + + //! triple point pressure [pa] + std::array p3_; + + //! saturation vapor pressure function: Vapor -> Cloud + SVPFunc1Container svp_func1_; + + //! cloud index set + std::vector cloud_index_set_; + + //! reaction information map + std::map cloud_reaction_map_; + + //! saturation vapor pressure function: Vapor + Vapor -> Cloud + SVPFunc2Container svp_func2_; + + //! pointer to the single Thermodynamics instance + static Thermodynamics *mythermo_; + + //! maximum saturation adjustment iterations + int sa_max_iter_ = 5; + + //! saturation adjustment relaxation parameter + Real sa_relax_ = 0.8; + + //! saturation adjustment temperature tolerance + Real sa_ftol_ = 0.01; +}; + +#endif // SRC_SNAP_THERMODYNAMICS_THERMODYNAMICS_HPP_ diff --git a/src/snap/eos/ideal_moist_hydro.cpp b/src/snap/eos/ideal_moist_hydro.cpp index 1b18a7bb..1f9b4f68 100644 --- a/src/snap/eos/ideal_moist_hydro.cpp +++ b/src/snap/eos/ideal_moist_hydro.cpp @@ -16,7 +16,6 @@ #include #include #include -#include // canoe #include @@ -27,7 +26,9 @@ #include // snap -#include "../thermodynamics/thermodynamics.hpp" +#include +#include + #include "eos_helper.hpp" // checks diff --git a/src/snap/thermodynamics/equilibrate_tp.cpp b/src/snap/thermodynamics/equilibrate_tp.cpp index b4f508b9..d13a3af7 100644 --- a/src/snap/thermodynamics/equilibrate_tp.cpp +++ b/src/snap/thermodynamics/equilibrate_tp.cpp @@ -1,76 +1,35 @@ -// cantera -#include -#include -#include - // canoe #include // snap #include "atm_thermodynamics.hpp" -void Thermodynamics::EquilibrateTP(AirParcel* air) const { - air->ToMassFraction(); - Real pres = air->w[IPR]; - for (int n = IVX; n < IVX + NCLOUD; ++n) air->w[n] = air->c[n - IVX]; - auto& thermo = kinetics_->thermo(); - thermo.setMassFractionsPartial(&air->w[1]); - thermo.setDensity(air->w[IDN]); - thermo.setPressure(air->w[IPR]); - - EquilibrateTP(); - - for (int n = 0; n < NCLOUD; ++n) air->c[n] = air->w[IVX + n]; - air->w[IPR] = pres; - air->ToMoleFraction(); -} - -void Thermodynamics::EquilibrateTP() const { - std::static_pointer_cast(kinetics_) - ->setQuantityMoleFraction(); +void Thermodynamics::EquilibrateTP(AirParcel* qfrac) const { + set_total_equivalent_vapor(qfrac, cloud_index_set_.data(), + cloud_reaction_map_); - std::vector rates(kinetics_->nTotalSpecies()); - std::vector mfrac(kinetics_->nTotalSpecies()); + // vapor <=> cloud + for (int i = 1; i <= NVAPOR; ++i) { + auto rates = TryEquilibriumTP_VaporCloud(*qfrac, i); - auto& thermo = kinetics_->thermo(); - Real temp = thermo.temperature(); - Real pres = thermo.pressure(); + // vapor condensation rate + qfrac->w[i] += rates[0]; - int iter = 0, max_iter = 3; - - while (iter++ < max_iter) { - // std::cout << "Iteration " << iter << std::endl; - - // get mole fraction - kinetics_->getActivityConcentrations(mfrac.data()); - - /* print initial conditions - std::cout << "Mole fractions" << std::endl; - for (size_t i = 0; i < mfrac.size(); ++i) { - std::cout << "Mole fraction " << i << ": " << mfrac[i] << std::endl; - }*/ - - kinetics_->getNetProductionRates(rates.data()); - /*for (size_t i = 0; i < rates.size(); ++i) { - std::cout << "NET Production " << i << ": " << rates[i] << std::endl; - }*/ + // cloud concentration rates + for (int n = 1; n < rates.size(); ++n) + qfrac->c[cloud_index_set_[i][n - 1]] += rates[n]; + } - // update mole fraction - for (size_t i = 1; i < rates.size(); ++i) { - mfrac[i] += rates[i]; - } + // vapor + vapor <=> cloud + for (auto const& [ij, info] : cloud_reaction_map_) { + auto rates = TryEquilibriumTP_VaporVaporCloud(*qfrac, ij); + auto& indx = info.first; - thermo.setMoleFractions(mfrac.data()); - thermo.setTemperature(temp); - thermo.setPressure(pres); + // vapor condensation rate + qfrac->w[indx[0]] += rates[0]; + qfrac->w[indx[1]] += rates[1]; - /* print again - std::cout << "After Mole fractions" << std::endl; - for (size_t i = 0; i < mfrac.size(); ++i) { - std::cout << "Mole fraction " << i << ": " << mfrac[i] << std::endl; - } - std::cout << "T = " << thermo.temperature() << std::endl; - std::cout << "P = " << thermo.pressure() << std::endl; - std::cout << "D = " << thermo.density() << std::endl;*/ + // cloud concentration rates + qfrac->c[indx[2]] += rates[3]; } } diff --git a/src/snap/thermodynamics/rk4_integrate_z.cpp b/src/snap/thermodynamics/rk4_integrate_z.cpp index 4b398718..e9ce735e 100644 --- a/src/snap/thermodynamics/rk4_integrate_z.cpp +++ b/src/snap/thermodynamics/rk4_integrate_z.cpp @@ -3,58 +3,37 @@ #include #include -// cantera -#include -#include - // thermodynamics #include "atm_thermodynamics.hpp" void rk4_integrate_z(AirParcel* air, Real dz, std::string method, Real grav, Real adTdz) { auto pthermo = Thermodynamics::GetInstance(); - auto& thermo = pthermo->Kinetics()->thermo(); auto const& mu_ratio = pthermo->GetMuRatio(); auto const& cp_ratio_mole = pthermo->GetCpRatioMole(); Real step[] = {0.5, 0.5, 1.}; - Real dTdz[4], chi[4]; - Real latent[1 + NVAPOR]; Real temp = air->w[IDN]; - - std::vector rates(pthermo->Kinetics()->nTotalSpecies()); - std::vector enthalpy(pthermo->Kinetics()->nTotalSpecies()); - - air->ToMassFraction(); Real pres = air->w[IPR]; - for (int n = IVX; n < IVX + NCLOUD; ++n) air->w[n] = air->c[n - IVX]; - thermo.setMassFractionsPartial(&air->w[1]); - thermo.setDensity(air->w[IDN]); - thermo.setPressure(air->w[IPR]); - - thermo.getEnthalpy_RT(enthalpy.data()); - pthermo->Kinetics()->getNetProductionRates(rates.data()); - - for (int n = 0; n < NCLOUD; ++n) air->c[n] = air->w[IVX + n]; - air->w[IPR] = pres; - air->ToMoleFraction(); - - for (int i = 1; i <= NVAPOR; ++i) { - if (rates[i] > 0.) { - latent[i] = enthalpy[i] - enthalpy[i + NVAPOR]; - } else { - latent[i] = 0.; - } - } + Real dTdz[4], chi[4]; + Real latent[1 + NVAPOR]; for (int rk = 0; rk < 4; ++rk) { pthermo->EquilibrateTP(air); - if (method != "reversible") { for (int j = 0; j < NCLOUD; ++j) air->c[j] = 0; } + for (int i = 1; i <= NVAPOR; ++i) { + // make a trial run to get latent heat + air->w[i] += 1.E-6; + auto rates = pthermo->TryEquilibriumTP_VaporCloud(*air, i); + latent[i] = pthermo->GetLatentHeatMole(i, rates, air->w[IDN]) / + (Constants::Rgas * air->w[IDN]); + air->w[i] -= 1.E-6; + } + Real q_gas = 1., q_eps = 1.; for (int i = 1; i <= NVAPOR; ++i) { q_eps += air->w[i] * (mu_ratio[i] - 1.); diff --git a/src/snap/thermodynamics/saturation_adjustment.cpp b/src/snap/thermodynamics/saturation_adjustment.cpp index 8bca5780..c330e9dc 100644 --- a/src/snap/thermodynamics/saturation_adjustment.cpp +++ b/src/snap/thermodynamics/saturation_adjustment.cpp @@ -22,6 +22,80 @@ void Thermodynamics::SaturationAdjustment(AirColumn& air_column) const { if (NVAPOR == 0) return; for (auto& air : air_column) { - EquilibrateUV(&air); + air.ToMoleFraction(); + AirParcel air0(air); + + // calculate total mole density + Real rmole = get_density_mole(air); + Real umole = get_internal_energy_mole(air, cv_ratio_mole_.data(), + latent_energy_mole_.data(), + cloud_index_set_.data()); + + int iter = 0; + + Real Teq = air.w[IDN]; + std::stringstream msg; + + while (iter++ < sa_max_iter_) { + // msg << "iter = " << iter << std::endl; + // msg << "old var = " << air << std::endl; + + Real cvd = 1. / (get_gammad(air) - 1.); + Real fsig = 1.; + + bool isNeg = false; + + 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; + } + +#pragma omp simd reduction(+ : qv) + for (auto j : cloud_index_set_[i]) qv += air.c[j]; + + fsig += (cv_ratio_mole_[i] - 1.) * qv; + } + + // condensation + for (int i = 1; i <= NVAPOR; ++i) { + auto rates = TryEquilibriumTP_VaporCloud(air, i, cvd * fsig); + + // 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; + } + // cloud concentration rates + for (int j = 1; j < rates.size(); ++j) { + air.c[cloud_index_set_[i][j - 1]] += rates[j]; + } + } + + Real Told = air.w[IDN]; + update_TP_conserving_U(&air, rmole, umole, cv_ratio_mole_.data(), + latent_energy_mole_.data(), + cloud_index_set_.data()); + // msg << "new var = " << air << std::endl; + if (fabs(air.w[IDN] - Teq) < sa_ftol_) break; + + // relax temperature and pressure + Teq = air.w[IDN]; + Real qgas = 1.; +#pragma omp simd reduction(+ : qgas) + for (int n = 0; n < NCLOUD; ++n) qgas += -air.c[n]; + + air.w[IDN] = Told + sa_relax_ * (air.w[IDN] - Told); + air.w[IPR] = rmole * qgas * Constants::Rgas * air.w[IDN]; + } + + if (iter > sa_max_iter_) { + msg << "Variables before iteration q0 = (" << air0 << ")" << std::endl; + msg << "Variables after iteration q = (" << air << ")" << std::endl; + throw RuntimeError("SaturationAdjustment", msg.str()); + } } } diff --git a/src/snap/thermodynamics/thermodynamics.cpp b/src/snap/thermodynamics/thermodynamics.cpp index a1416e3f..c49c28f4 100644 --- a/src/snap/thermodynamics/thermodynamics.cpp +++ b/src/snap/thermodynamics/thermodynamics.cpp @@ -4,15 +4,13 @@ #include #include #include -#include #include #include #include #include // fill -// cantera -#include -#include +// external +#include // athena #include @@ -48,58 +46,23 @@ Thermodynamics::~Thermodynamics() { app->Log("Destroy Thermodynamics"); } -Thermodynamics* Thermodynamics::fromYAMLInput(std::string const& fname) { +Thermodynamics* Thermodynamics::fromYAMLInput(YAML::Node const& node) { mythermo_ = new Thermodynamics(); - auto thermo = Cantera::newThermo(fname, "gas"); - - if (thermo->nSpecies() != 1 + NVAPOR + NCLOUD) { - throw RuntimeError("Thermodynamics", - "Number of species does not match the input file"); - } - - auto& kinetics = mythermo_->kinetics_; - kinetics = Cantera::newKinetics({thermo}, fname); - - // --------- vapor + cloud thermo --------- - std::vector mu(thermo->nSpecies()); - std::vector cp_mole(thermo->nSpecies()); - - // g/mol - thermo->getMolecularWeights(mu.data()); - - // J/kmol/K - thermo->getPartialMolarCp(cp_mole.data()); - - mythermo_->Rd_ = Cantera::GasConstant / mu[0]; - mythermo_->gammad_ref_ = cp_mole[0] / (cp_mole[0] - Cantera::GasConstant); - - // ---------- dimensionless properties ---------- - - // cp ratios - for (size_t i = 0; i < cp_mole.size(); ++i) { - mythermo_->cp_ratio_mole_[i] = cp_mole[i] / cp_mole[0]; - } - - // molecular weight ratios - for (size_t i = 0; i < mu.size(); ++i) { - mythermo_->mu_ratio_[i] = mu[i] / mu[0]; - } - - // beta parameter (dimensionless internal energy) - for (size_t i = 0; i < thermo->nSpecies(); ++i) { - auto& tp = thermo->species(i)->thermo; - size_t n; - int type; - Real tlow, thigh, pref; - std::vector coeffs(tp->nCoeffs()); - tp->reportParameters(n, type, tlow, thigh, pref, coeffs.data()); - - Real t0 = coeffs[0]; - Real h0 = coeffs[1]; - Real s0 = coeffs[2]; - Real cp0 = coeffs[3]; - mythermo_->beta_[i] = h0 / (Cantera::GasConstant * t0); + // non-condensable + mythermo_->Rd_ = node["dry"]["Rd"].as(); + mythermo_->gammad_ref_ = node["dry"]["gammad"].as(); + + // saturation adjustment + if (node["saturation_adjustment"]) { + mythermo_->sa_max_iter_ = + node["saturation_adjustment"]["max_iter"].as(10); + mythermo_->sa_ftol_ = node["saturation_adjustment"]["ftol"].as(1e-4); + mythermo_->sa_relax_ = node["saturation_adjustment"]["relax"].as(0.8); + } else { + mythermo_->sa_max_iter_ = 10; + mythermo_->sa_ftol_ = 1e-4; + mythermo_->sa_relax_ = 0.8; } return mythermo_; @@ -130,8 +93,8 @@ Thermodynamics* Thermodynamics::fromLegacyInput(ParameterInput* pin) { // Read triple point pressure read_thermo_property(mythermo_->p3_.data(), "Ptriple", 1, 0., pin); - mythermo_->svp_func1_.resize(1 + NVAPOR); mythermo_->cloud_index_set_.resize(1 + NVAPOR); + mythermo_->svp_func1_.resize(1 + NVAPOR); for (int i = 0; i <= NVAPOR; ++i) { mythermo_->cloud_index_set_[i].resize(NPHASE - 1); @@ -159,7 +122,7 @@ Thermodynamics const* Thermodynamics::GetInstance() { } Thermodynamics const* Thermodynamics::InitFromYAMLInput( - std::string const& fname) { + YAML::Node const& node) { if (mythermo_ != nullptr) { throw RuntimeError("Thermodynamics", "Thermodynamics has been initialized"); } @@ -167,7 +130,7 @@ Thermodynamics const* Thermodynamics::InitFromYAMLInput( Application::Logger app("snap"); app->Log("Initialize Thermodynamics"); - return fromYAMLInput(fname); + return fromYAMLInput(node); } Thermodynamics const* Thermodynamics::InitFromAthenaInput(ParameterInput* pin) { @@ -180,11 +143,20 @@ Thermodynamics const* Thermodynamics::InitFromAthenaInput(ParameterInput* pin) { if (pin->DoesParameterExist("thermodynamics", input_key)) { std::string filename = pin->GetString("thermodynamics", input_key); - mythermo_ = fromYAMLInput(filename); + std::ifstream stream(filename); + if (stream.good() == false) { + throw RuntimeError("Thermodynamics", + "Cannot open thermodynamic file: " + filename); + } + + mythermo_ = fromYAMLInput(YAML::Load(stream)); } else { // legacy input mythermo_ = fromLegacyInput(pin); } + // enroll vapor functions + mythermo_->enrollVaporFunctions(); + // alias auto& Rd = mythermo_->Rd_; auto& gammad = mythermo_->gammad_ref_; @@ -215,10 +187,10 @@ Thermodynamics const* Thermodynamics::InitFromAthenaInput(ParameterInput* pin) { inv_mu_ratio[n] = 1. / mu_ratio[n]; mu[n] = mu[0] * mu_ratio[n]; inv_mu[n] = 1. / mu[n]; - cp_ratio_mass[n] = cp_ratio_mole[n] * inv_mu_ratio[n]; + cp_ratio_mole[n] = cp_ratio_mass[n] * mu_ratio[n]; } - /* calculate latent energy = $\beta\frac{R_d}{\epsilon}T^r$ + // calculate latent energy = $\beta\frac{R_d}{\epsilon}T^r$ for (int n = 0; n <= NVAPOR; ++n) { latent_energy_mass[n] = 0.; latent_energy_mole[n] = 0.; @@ -240,13 +212,36 @@ Thermodynamics const* Thermodynamics::InitFromAthenaInput(ParameterInput* pin) { delta[n] = (cp_ratio_mass[n] - cp_ratio_mass[i]) * mu_ratio[i] / (1. - 1. / gammad); } - }*/ + } + + // set up extra clouds + auto pindex = IndexMap::GetInstance(); + + for (int n = (NPHASE - 1) * NVAPOR; n < NCLOUD; ++n) { + Molecule mol; + mol.LoadThermodynamicFile(pindex->GetCloudName(n) + ".thermo"); + + int j = 1 + NVAPOR + n; + mu[j] = mol.mu(); + inv_mu[j] = 1. / mu[j]; + + mu_ratio[j] = mol.mu() / mu[0]; + inv_mu_ratio[j] = 1. / mu_ratio[j]; + + cp_ratio_mass[j] = mol.cp() / (mol.mu() * mythermo_->GetCpMassRef(0)); + cp_ratio_mole[j] = cp_ratio_mass[j] * mu_ratio[j]; + + latent_energy_mole[j] = mol.latent(); + latent_energy_mass[j] = mol.latent() / mol.mu(); + + beta[j] = mol.latent() / (Constants::Rgas * Thermodynamics::RefTemp); + delta[j] = 0.; + } // finally, set up cv ratio // calculate cv_ratio = $\sigma_i + (1. - \gamma)/\epsilon_i$ for (int n = 0; n <= NVAPOR; ++n) { - cv_ratio_mass[n] = - gammad * cp_ratio_mass[n] + (1. - gammad) * inv_mu_ratio[n]; + cv_ratio_mass[n] = gammad * cp_ratio_mass[n] + (1. - gammad) / mu_ratio[n]; cv_ratio_mole[n] = cv_ratio_mass[n] * mu_ratio[n]; } @@ -410,34 +405,20 @@ Real Thermodynamics::RelativeHumidity(MeshBlock const* pmb, int n, int k, int j, void Thermodynamics::Extrapolate(AirParcel* qfrac, Real dzORdlnp, std::string method, Real grav, Real userp) const { - qfrac->ToMassFraction(); - - auto& thermo = kinetics_->thermo(); - - double pres = qfrac->w[IPR]; - for (int n = IVX; n < IVX + NCLOUD; ++n) qfrac->w[n] = qfrac->c[n - IVX]; - thermo.setMassFractionsPartial(&qfrac->w[1]); - thermo.setDensity(qfrac->w[IDN]); - thermo.setPressure(qfrac->w[IPR]); + qfrac->ToMoleFraction(); - EquilibrateTP(); + // equilibrate vapor with clouds + for (int i = 1; i <= NVAPOR; ++i) { + auto rates = TryEquilibriumTP_VaporCloud(*qfrac, i); - thermo.getMassFractions(&qfrac->w[0]); + // vapor condensation rate + qfrac->w[i] += rates[0]; - for (int n = 0; n < NCLOUD; ++n) { - qfrac->c[n] = qfrac->w[IVX + n]; + // cloud concentration rates + for (int j = 1; j < rates.size(); ++j) + qfrac->c[cloud_index_set_[i][j - 1]] += rates[j]; } - qfrac->w[IPR] = pres; - qfrac->w[IDN] = thermo.density(); - qfrac->w[IVX] = 0.; - qfrac->w[IVX] = 0.; - qfrac->w[IVX] = 0.; - - qfrac->ToMoleFraction(); - - // std::cout << *qfrac << std::endl; - // RK4 integration #ifdef HYDROSTATIC rk4_integrate_lnp(qfrac, dzORdlnp, method, userp); @@ -459,34 +440,6 @@ Real Thermodynamics::GetLatentHeatMole(int i, std::vector const& rates, return heat / std::abs(rates[0]); } -void Thermodynamics::SetStateFromPrimitive(MeshBlock* pmb, int k, int j, - int i) const { - Hydro* phydro = pmb->phydro; - - auto& thermo = kinetics_->thermo(); - size_t total_cells = pmb->ncells1 * pmb->ncells2 * pmb->ncells3; - - thermo.setMassFractionsPartial(&phydro->w(1, k, j, i), total_cells); - thermo.setDensity(phydro->w(IDN, k, j, i)); - thermo.setPressure(phydro->w(IPR, k, j, i)); -} - -void Thermodynamics::SetStateFromConserved(MeshBlock* pmb, int k, int j, - int i) const { - Hydro* phydro = pmb->phydro; - - auto& thermo = kinetics_->thermo(); - size_t total_cells = pmb->ncells1 * pmb->ncells2 * pmb->ncells3; - - thermo.setMassFractions(&phydro->u(0, k, j, i), total_cells); - Real rho = 0.; - for (int n = 0; n < Size; ++n) { - rho += phydro->u(n, k, j, i); - } - thermo.setDensity(rho); - thermo.setPressure(GetPres(pmb, k, j, i)); -} - // Eq.4.5.11 in Emanuel (1994) Real Thermodynamics::EquivalentPotentialTemp(MeshBlock* pmb, Real p0, int v, int k, int j, int i) const { diff --git a/src/snap/thermodynamics/thermodynamics.hpp b/src/snap/thermodynamics/thermodynamics.hpp index 33109270..8f53f329 100644 --- a/src/snap/thermodynamics/thermodynamics.hpp +++ b/src/snap/thermodynamics/thermodynamics.hpp @@ -11,6 +11,9 @@ #include #include +// external +#include + // athena #include #include @@ -24,11 +27,6 @@ class MeshBlock; class ParameterInput; -namespace Cantera { -class ThermoPhase; -class Kinetics; -} // namespace Cantera - using IndexPair = std::pair; using IndexSet = std::vector; @@ -38,8 +36,8 @@ using RealArrayX = std::vector; using SatVaporPresFunc1 = Real (*)(AirParcel const &, int i, int j); using SatVaporPresFunc2 = Real (*)(AirParcel const &, int i, int j, int k); +//! \todo(CLI): move to configure.hpp enum { MAX_REACTANT = 3 }; - using ReactionIndx = std::array; using ReactionStoi = std::array; using ReactionInfo = std::pair; @@ -78,7 +76,7 @@ class Thermodynamics { //! Protected ctor access thru static member function Instance Thermodynamics() {} static Thermodynamics *fromLegacyInput(ParameterInput *pin); - static Thermodynamics *fromYAMLInput(std::string const &fname); + static Thermodynamics *fromYAMLInput(YAML::Node const &node); public: using SVPFunc1Container = std::vector>; @@ -86,6 +84,9 @@ class Thermodynamics { enum { Size = 1 + NVAPOR + NCLOUD }; + static constexpr Real RefTemp = 300.; + static constexpr Real RefPres = 1.e5; + //! thermodynamics input key in the input file [thermodynamics_config] static const std::string input_key; @@ -95,7 +96,7 @@ class Thermodynamics { //! Return a pointer to the one and only instance of Thermodynamics static Thermodynamics const *GetInstance(); - static Thermodynamics const *InitFromYAMLInput(std::string const &fname); + static Thermodynamics const *InitFromYAMLInput(YAML::Node const &node); static Thermodynamics const *InitFromAthenaInput(ParameterInput *pin); @@ -148,8 +149,6 @@ class Thermodynamics { //! const pointer to cv_ratio_mole_ Real const *GetCvRatioMole() const { return cv_ratio_mole_.data(); } - std::shared_ptr Kinetics() const { return kinetics_; } - //! \return the index of the cloud //! \param[in] i the index of the vapor //! \param[in] j the sequential index of the cloud @@ -216,6 +215,10 @@ class Thermodynamics { return GetLatentHeatMole(i, rates, temp) * inv_mu_[i]; } + RealArrayX CalcSurfEvapRates(AirParcel const &qfrac, int i, Real &amd, + Real btemp, Real dTs, Real cSurf, Real dt, + Real Cde, Real Mbar) const; + //! \brief Calculate the equilibrium mole transfer by cloud reaction //! vapor -> cloud //! @@ -251,25 +254,13 @@ class Thermodynamics { //! Thermodnamic equilibrium at current TP //! \param[in,out] qfrac mole fraction representation of air parcel - void EquilibrateTP() const; - void EquilibrateTP(AirParcel *qfrac) const; - //! Thermodnamic equilibrium at current UV - void EquilibrateUV() const; - - void EquilibrateUV(AirParcel *qfrac) const; - - void EquilibrateSP(double P) const; - //! Adjust to the maximum saturation state conserving internal energy //! \param[in,out] ac mole fraction representation of a collection of air //! parcels void SaturationAdjustment(AirColumn &ac) const; - void SetStateFromPrimitive(MeshBlock *pmb, int k, int j, int i) const; - void SetStateFromConserved(MeshBlock *pmb, int k, int j, int i) const; - //! \brief Calculate potential temperature from primitive variable //! //! $\theta = T(\frac{p_0}{p})^{\chi}$ @@ -372,8 +363,6 @@ class Thermodynamics { void enrollVaporFunctions(); protected: - std::shared_ptr kinetics_; - //! ideal gas constant of dry air in J/kg Real Rd_; @@ -441,12 +430,12 @@ class Thermodynamics { //! cloud index set std::vector cloud_index_set_; - //! reaction information map - std::map cloud_reaction_map_; - //! saturation vapor pressure function: Vapor + Vapor -> Cloud SVPFunc2Container svp_func2_; + //! reaction information map + std::map cloud_reaction_map_; + //! pointer to the single Thermodynamics instance static Thermodynamics *mythermo_; From eb71b2fc41a2d0143df55af968beeb91e4b4dfc2 Mon Sep 17 00:00:00 2001 From: mac/cli Date: Mon, 8 Jul 2024 12:08:07 -0400 Subject: [PATCH 2/8] wip --- .github/workflows/mac.yml | 2 +- .github/workflows/main.yml | 23 ++++++++++++++++++++++- CMakeLists.txt | 6 ++++++ 3 files changed, 29 insertions(+), 2 deletions(-) diff --git a/.github/workflows/mac.yml b/.github/workflows/mac.yml index ec9ef614..7dd8e7d8 100644 --- a/.github/workflows/mac.yml +++ b/.github/workflows/mac.yml @@ -338,7 +338,7 @@ jobs: run: brew bundle - name: install python modules - run: sudo pip3 install -r requirements.txt + run: pip3 install -r requirements.txt - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index e8fe6f0b..8966327d 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -49,6 +49,9 @@ jobs: with: lfs: false + - name: install python modules + run: sudo pip install -r requirements.txt + - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type @@ -73,6 +76,9 @@ jobs: with: lfs: false + - name: install python modules + run: sudo pip install -r requirements.txt + - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type @@ -100,6 +106,9 @@ jobs: with: lfs: false + - name: install python modules + run: sudo pip install -r requirements.txt + - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type @@ -194,6 +203,9 @@ jobs: with: lfs: false + - name: install python modules + run: sudo pip install -r requirements.txt + - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type @@ -241,6 +253,9 @@ jobs: with: lfs: false + - name: install python modules + run: sudo pip install -r requirements.txt + - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type @@ -288,6 +303,9 @@ jobs: with: lfs: false + - name: install python modules + run: sudo pip install -r requirements.txt + - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. # See https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html?highlight=cmake_build_type @@ -685,7 +703,7 @@ jobs: lfs: false - name: install python modules - run: sudo pip3 install -r requirements.txt + run: pip3 install -r requirements.txt - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. @@ -717,6 +735,9 @@ jobs: with: lfs: false + - name: install python modules + run: pip3 install -r requirements.txt + - name: check if artifact exists run: echo "ARTIFACT_EXISTS=$(./tools/check_artifact.sh test_blast_output.nc)" >> $GITHUB_ENV diff --git a/CMakeLists.txt b/CMakeLists.txt index 40b61897..4e349b80 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -61,6 +61,12 @@ find_package(Eigen3 REQUIRED) find_package(Torch REQUIRED) find_package(Cantera) +if (NOT ${CANTERA_FOUND} STREQUAL "TRUE") + message(STATUS "Not using Cantera library") + set(CANTERA_INCLUDE_DIR "") + set(CANTERA_LIBRARY "") +endif() + ## 3. set up project system libraries ## message(STATUS "3. Set up system libraries") if (NOT ${CANTERA_FOUND} STREQUAL "TRUE") From 575e9c4a0495b52d427a7a42e004277521c695f6 Mon Sep 17 00:00:00 2001 From: mac/cli Date: Mon, 8 Jul 2024 12:20:45 -0400 Subject: [PATCH 3/8] wip --- .github/workflows/mac.yml | 6 ++++++ .github/workflows/main.yml | 24 ++++++++++++------------ CMakeLists.txt | 6 ------ cmake/modules/FindCantera.cmake | 4 ++++ 4 files changed, 22 insertions(+), 18 deletions(-) diff --git a/.github/workflows/mac.yml b/.github/workflows/mac.yml index 7dd8e7d8..fc84d1f0 100644 --- a/.github/workflows/mac.yml +++ b/.github/workflows/mac.yml @@ -22,6 +22,9 @@ jobs: with: lfs: false + - name: set up python libraries + run: pip3 install -r requirements.txt + - name: check if artifact exists run: echo "ARTIFACT_EXISTS=$(./tools/check_artifact.sh straka_output.nc)" >> $GITHUB_ENV @@ -222,6 +225,9 @@ jobs: with: lfs: false + - name: set up python libraries + run: pip3 install -r requirements.txt + - name: check if artifact exists run: echo "ARTIFACT_EXISTS=$(./tools/check_artifact.sh swxy_output.nc)" >> $GITHUB_ENV diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 8966327d..f0dfac21 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -25,7 +25,7 @@ jobs: - run: echo "🖥️ The workspace, ${{ github.workspace }}, is now ready to test your code on the runner." - name: Install cpplint - run: sudo pip install cpplint + run: pip install cpplint - uses: actions/checkout@v3 with: @@ -50,7 +50,7 @@ jobs: lfs: false - name: install python modules - run: sudo pip install -r requirements.txt + run: pip install -r requirements.txt - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. @@ -77,7 +77,7 @@ jobs: lfs: false - name: install python modules - run: sudo pip install -r requirements.txt + run: pip install -r requirements.txt - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. @@ -107,7 +107,7 @@ jobs: lfs: false - name: install python modules - run: sudo pip install -r requirements.txt + run: pip install -r requirements.txt - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. @@ -155,7 +155,7 @@ jobs: lfs: false - name: install python modules - run: sudo pip install -r requirements.txt + run: pip install -r requirements.txt - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. @@ -204,7 +204,7 @@ jobs: lfs: false - name: install python modules - run: sudo pip install -r requirements.txt + run: pip install -r requirements.txt - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. @@ -254,7 +254,7 @@ jobs: lfs: false - name: install python modules - run: sudo pip install -r requirements.txt + run: pip install -r requirements.txt - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. @@ -304,7 +304,7 @@ jobs: lfs: false - name: install python modules - run: sudo pip install -r requirements.txt + run: pip install -r requirements.txt - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. @@ -353,7 +353,7 @@ jobs: lfs: false - name: install python modules - run: sudo pip install -r requirements.txt + run: pip install -r requirements.txt - name: create build directory # Configure CMake in a 'build' subdirectory. `CMAKE_BUILD_TYPE` is only required if you are using a single-configuration generator such as make. @@ -385,6 +385,9 @@ jobs: with: lfs: false + - name: install python modules + run: pip install -r requirements.txt + - name: check if artifact exists run: echo "ARTIFACT_EXISTS=$(./tools/check_artifact.sh robert_output.nc)" >> $GITHUB_ENV @@ -398,9 +401,6 @@ jobs: if: env.ARTIFACT_EXISTS == 'false' run: git lfs pull -I examples/2019-Li-snap/straka_output.nc - - name: set up python libraries - run: pip3 install -r requirements.txt - - name: create build directory run: cmake -B ${{github.workspace}}/build -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} -DTASK=straka diff --git a/CMakeLists.txt b/CMakeLists.txt index 4e349b80..40b61897 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -61,12 +61,6 @@ find_package(Eigen3 REQUIRED) find_package(Torch REQUIRED) find_package(Cantera) -if (NOT ${CANTERA_FOUND} STREQUAL "TRUE") - message(STATUS "Not using Cantera library") - set(CANTERA_INCLUDE_DIR "") - set(CANTERA_LIBRARY "") -endif() - ## 3. set up project system libraries ## message(STATUS "3. Set up system libraries") if (NOT ${CANTERA_FOUND} STREQUAL "TRUE") diff --git a/cmake/modules/FindCantera.cmake b/cmake/modules/FindCantera.cmake index ef338a53..fc412cd3 100644 --- a/cmake/modules/FindCantera.cmake +++ b/cmake/modules/FindCantera.cmake @@ -49,4 +49,8 @@ if(CANTERA_INCLUDE_DIR) message(STATUS "Symbolic link created: ${LINK_PATH} -> ${SOURCE_PATH}") endif() include_directories(${CMAKE_CURRENT_BINARY_DIR}/ext1) +else() + message(STATUS "Not using Cantera library") + set(CANTERA_INCLUDE_DIR "") + set(CANTERA_LIBRARY "") endif() From b4ce97321551740623cb3bf11b506d74e6f2aa0e Mon Sep 17 00:00:00 2001 From: mac/cli Date: Mon, 8 Jul 2024 14:40:19 -0400 Subject: [PATCH 4/8] fix lib --- .github/workflows/mac.yml | 3 +++ CMakeLists.txt | 3 ++- cmake/modules/FindCantera.cmake | 1 + cmake/yamlpp.cmake | 2 ++ src/CMakeLists.txt | 7 ++++++- 5 files changed, 14 insertions(+), 2 deletions(-) diff --git a/.github/workflows/mac.yml b/.github/workflows/mac.yml index fc84d1f0..c2afc453 100644 --- a/.github/workflows/mac.yml +++ b/.github/workflows/mac.yml @@ -375,6 +375,9 @@ jobs: with: lfs: false + - name: install python modules + run: pip3 install -r requirements.txt + - name: check if artifact exists run: echo "ARTIFACT_EXISTS=$(./tools/check_artifact.sh test_blast_output.nc)" >> $GITHUB_ENV diff --git a/CMakeLists.txt b/CMakeLists.txt index 40b61897..c169162e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -63,7 +63,8 @@ find_package(Cantera) ## 3. set up project system libraries ## message(STATUS "3. Set up system libraries") -if (NOT ${CANTERA_FOUND} STREQUAL "TRUE") +message(STATUS ${CANTERA_FOUND}) +if (NOT ${CANTERA_FOUND}) include(${CMAKE_SOURCE_DIR}/cmake/yamlpp.cmake) endif() include(${CMAKE_SOURCE_DIR}/cmake/gtest.cmake) diff --git a/cmake/modules/FindCantera.cmake b/cmake/modules/FindCantera.cmake index fc412cd3..fa791dd6 100644 --- a/cmake/modules/FindCantera.cmake +++ b/cmake/modules/FindCantera.cmake @@ -53,4 +53,5 @@ else() message(STATUS "Not using Cantera library") set(CANTERA_INCLUDE_DIR "") set(CANTERA_LIBRARY "") + set(CANTERA_FOUND FALSE) endif() diff --git a/cmake/yamlpp.cmake b/cmake/yamlpp.cmake index 4397b473..a1855740 100644 --- a/cmake/yamlpp.cmake +++ b/cmake/yamlpp.cmake @@ -19,3 +19,5 @@ endif() # Where yaml-cpp's .h files can be found. include_directories(${yaml-cpp_SOURCE_DIR}/include) + +set(YAML_CPP_LIBRARIES yaml-cpp) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 31a190d5..07092247 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -33,6 +33,11 @@ add_library(${namel}_${buildl} ${src_files} ) +target_include_directories(${namel}_${buildl} + PRIVATE + ${EIGEN3_INCLUDE_DIR} + ) + set_target_properties(${namel}_${buildl} PROPERTIES COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}} @@ -82,7 +87,7 @@ set(CANOE_LIBRARY_${buildu} "diagnostics_${buildl}" "forcing_${buildl}" "modules_${buildl}" - #"yaml-cpp" + ${YAML_CPP_LIBRARIES} ${CPPDISORT_LIBRARY_${buildu}} ${PYTHON_LIBRARY_RELEASE} ${NETCDF_LIBRARIES} From 0c267d1e7d4f4cf1880075ba5efda46ca01932b6 Mon Sep 17 00:00:00 2001 From: stormy/cli Date: Mon, 8 Jul 2024 15:56:47 -0400 Subject: [PATCH 5/8] move torch --- CMakeLists.txt | 2 +- src/CMakeLists.txt | 8 ++++-- src/torch/CMakeLists.txt | 28 +++++++++++++++++++ .../decomposition => torch/column}/decom.cpp | 0 .../decomposition => torch/column}/decom.hpp | 0 src/{snap => torch}/eos/eos.cpp | 0 src/{snap => torch}/eos/eos.hpp | 0 src/{snap => torch}/eos/eos_hydro_ideal.cpp | 0 src/{ => torch}/exo3/exo3.cpp | 0 src/{ => torch}/exo3/exo3.hpp | 0 src/{snap => torch}/hydro/flux_divergence.cpp | 0 src/{snap => torch}/hydro/hydro.hpp | 0 .../recon}/interpolation.cpp | 0 .../recon}/interpolation.hpp | 0 .../reconstruct => torch/recon}/recon.hpp | 0 .../reconstruct => torch/recon}/weno5.cpp | 0 src/{snap => torch}/riemann/riemann.cpp | 0 src/{snap => torch}/riemann/riemann.hpp | 0 .../riemann/rs_hydro_lmars.cpp | 0 tests/CMakeLists.txt | 12 ++++---- 20 files changed, 42 insertions(+), 8 deletions(-) create mode 100644 src/torch/CMakeLists.txt rename src/{snap/decomposition => torch/column}/decom.cpp (100%) rename src/{snap/decomposition => torch/column}/decom.hpp (100%) rename src/{snap => torch}/eos/eos.cpp (100%) rename src/{snap => torch}/eos/eos.hpp (100%) rename src/{snap => torch}/eos/eos_hydro_ideal.cpp (100%) rename src/{ => torch}/exo3/exo3.cpp (100%) rename src/{ => torch}/exo3/exo3.hpp (100%) rename src/{snap => torch}/hydro/flux_divergence.cpp (100%) rename src/{snap => torch}/hydro/hydro.hpp (100%) rename src/{snap/reconstruct => torch/recon}/interpolation.cpp (100%) rename src/{snap/reconstruct => torch/recon}/interpolation.hpp (100%) rename src/{snap/reconstruct => torch/recon}/recon.hpp (100%) rename src/{snap/reconstruct => torch/recon}/weno5.cpp (100%) rename src/{snap => torch}/riemann/riemann.cpp (100%) rename src/{snap => torch}/riemann/riemann.hpp (100%) rename src/{snap => torch}/riemann/rs_hydro_lmars.cpp (100%) diff --git a/CMakeLists.txt b/CMakeLists.txt index c169162e..10ede10f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -58,8 +58,8 @@ endif() message(STATUS "Include ${CMAKE_CURRENT_SOURCE_DIR}/cmake/parameters.cmake") include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/parameters.cmake) find_package(Eigen3 REQUIRED) -find_package(Torch REQUIRED) find_package(Cantera) +#find_package(Torch) ## 3. set up project system libraries ## message(STATUS "3. Set up system libraries") diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 07092247..53cdf7b5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -61,7 +61,11 @@ add_subdirectory(special) add_subdirectory(single_column) add_subdirectory(diagnostics) add_subdirectory(forcing) -add_subdirectory(modules) + +if (${Torch_FOUND}) + add_subdirectory(torch) + add_subdirectory(modules) +endif() #add_subdirectory(nbody) set(CANOE_LIBRARY_${buildu} @@ -86,7 +90,7 @@ set(CANOE_LIBRARY_${buildu} "scm_${buildl}" "diagnostics_${buildl}" "forcing_${buildl}" - "modules_${buildl}" + #"modules_${buildl}" ${YAML_CPP_LIBRARIES} ${CPPDISORT_LIBRARY_${buildu}} ${PYTHON_LIBRARY_RELEASE} diff --git a/src/torch/CMakeLists.txt b/src/torch/CMakeLists.txt new file mode 100644 index 00000000..3d425eaf --- /dev/null +++ b/src/torch/CMakeLists.txt @@ -0,0 +1,28 @@ +set(namel torch) +string(TOUPPER ${namel} nameu) + +file(GLOB src_files + eos/*.cpp + hydro/*.cpp + riemann/*.cpp + recon/*.cpp + column/*.cpp + ) + +string(TOLOWER ${CMAKE_BUILD_TYPE} buildl) +string(TOUPPER ${CMAKE_BUILD_TYPE} buildu) + +add_library(${namel}_${buildl} + OBJECT + ${src_files} + ) + +set_target_properties(${namel}_${buildl} + PROPERTIES + COMPILE_FLAGS ${CMAKE_CXX_FLAGS_${buildu}} + ) + +target_include_directories(${namel}_${buildl} + PRIVATE + ${EIGEN3_INCLUDE_DIR} + ) diff --git a/src/snap/decomposition/decom.cpp b/src/torch/column/decom.cpp similarity index 100% rename from src/snap/decomposition/decom.cpp rename to src/torch/column/decom.cpp diff --git a/src/snap/decomposition/decom.hpp b/src/torch/column/decom.hpp similarity index 100% rename from src/snap/decomposition/decom.hpp rename to src/torch/column/decom.hpp diff --git a/src/snap/eos/eos.cpp b/src/torch/eos/eos.cpp similarity index 100% rename from src/snap/eos/eos.cpp rename to src/torch/eos/eos.cpp diff --git a/src/snap/eos/eos.hpp b/src/torch/eos/eos.hpp similarity index 100% rename from src/snap/eos/eos.hpp rename to src/torch/eos/eos.hpp diff --git a/src/snap/eos/eos_hydro_ideal.cpp b/src/torch/eos/eos_hydro_ideal.cpp similarity index 100% rename from src/snap/eos/eos_hydro_ideal.cpp rename to src/torch/eos/eos_hydro_ideal.cpp diff --git a/src/exo3/exo3.cpp b/src/torch/exo3/exo3.cpp similarity index 100% rename from src/exo3/exo3.cpp rename to src/torch/exo3/exo3.cpp diff --git a/src/exo3/exo3.hpp b/src/torch/exo3/exo3.hpp similarity index 100% rename from src/exo3/exo3.hpp rename to src/torch/exo3/exo3.hpp diff --git a/src/snap/hydro/flux_divergence.cpp b/src/torch/hydro/flux_divergence.cpp similarity index 100% rename from src/snap/hydro/flux_divergence.cpp rename to src/torch/hydro/flux_divergence.cpp diff --git a/src/snap/hydro/hydro.hpp b/src/torch/hydro/hydro.hpp similarity index 100% rename from src/snap/hydro/hydro.hpp rename to src/torch/hydro/hydro.hpp diff --git a/src/snap/reconstruct/interpolation.cpp b/src/torch/recon/interpolation.cpp similarity index 100% rename from src/snap/reconstruct/interpolation.cpp rename to src/torch/recon/interpolation.cpp diff --git a/src/snap/reconstruct/interpolation.hpp b/src/torch/recon/interpolation.hpp similarity index 100% rename from src/snap/reconstruct/interpolation.hpp rename to src/torch/recon/interpolation.hpp diff --git a/src/snap/reconstruct/recon.hpp b/src/torch/recon/recon.hpp similarity index 100% rename from src/snap/reconstruct/recon.hpp rename to src/torch/recon/recon.hpp diff --git a/src/snap/reconstruct/weno5.cpp b/src/torch/recon/weno5.cpp similarity index 100% rename from src/snap/reconstruct/weno5.cpp rename to src/torch/recon/weno5.cpp diff --git a/src/snap/riemann/riemann.cpp b/src/torch/riemann/riemann.cpp similarity index 100% rename from src/snap/riemann/riemann.cpp rename to src/torch/riemann/riemann.cpp diff --git a/src/snap/riemann/riemann.hpp b/src/torch/riemann/riemann.hpp similarity index 100% rename from src/snap/riemann/riemann.hpp rename to src/torch/riemann/riemann.hpp diff --git a/src/snap/riemann/rs_hydro_lmars.cpp b/src/torch/riemann/rs_hydro_lmars.cpp similarity index 100% rename from src/snap/riemann/rs_hydro_lmars.cpp rename to src/torch/riemann/rs_hydro_lmars.cpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 44f86f2e..847365b3 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -8,11 +8,13 @@ configure_file(globals.cpp.in globals.cpp @ONLY) enable_testing() # torch tests -setup_test(test_eos) -setup_test(test_weno) -setup_test(test_reconstruct) -setup_test(test_riemann) -setup_test(test_hydro) +if (${Torch_FOUND}) + setup_test(test_eos) + setup_test(test_weno) + setup_test(test_reconstruct) + setup_test(test_riemann) + setup_test(test_hydro) +endif() # athena tests setup_test(test_mesh) From 9c90b53bd1942fcd6e65c16c740b0d0a073f909b Mon Sep 17 00:00:00 2001 From: stormy/cli Date: Mon, 8 Jul 2024 16:55:58 -0400 Subject: [PATCH 6/8] wip --- examples/2023-Li-saturn-vla/saturn_radio.cpp | 2 +- examples/2023-Li-uranus/uranus_mwr.cpp | 2 +- examples/2024-JHu-juno-mwr/juno_mwr.cpp | 2 +- src/CMakeLists.txt | 5 +++++ src/diagnostics/todo/convective_heat_flx.cpp | 2 +- src/snap/eos/shallow_yz_hydro.cpp | 2 +- tests/CMakeLists.txt | 2 +- 7 files changed, 11 insertions(+), 6 deletions(-) diff --git a/examples/2023-Li-saturn-vla/saturn_radio.cpp b/examples/2023-Li-saturn-vla/saturn_radio.cpp index ebe61af9..3c1c4d43 100644 --- a/examples/2023-Li-saturn-vla/saturn_radio.cpp +++ b/examples/2023-Li-saturn-vla/saturn_radio.cpp @@ -14,7 +14,6 @@ #include #include #include -#include // application #include @@ -32,6 +31,7 @@ #include // snap +#include #include #include #include diff --git a/examples/2023-Li-uranus/uranus_mwr.cpp b/examples/2023-Li-uranus/uranus_mwr.cpp index 071778f4..cb8cdb00 100644 --- a/examples/2023-Li-uranus/uranus_mwr.cpp +++ b/examples/2023-Li-uranus/uranus_mwr.cpp @@ -14,7 +14,6 @@ #include #include #include -#include // application #include @@ -32,6 +31,7 @@ #include // snap +#include #include #include #include diff --git a/examples/2024-JHu-juno-mwr/juno_mwr.cpp b/examples/2024-JHu-juno-mwr/juno_mwr.cpp index 5758524a..224495b8 100644 --- a/examples/2024-JHu-juno-mwr/juno_mwr.cpp +++ b/examples/2024-JHu-juno-mwr/juno_mwr.cpp @@ -14,7 +14,6 @@ #include #include #include -#include // application #include @@ -32,6 +31,7 @@ #include // snap +#include #include #include diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 53cdf7b5..a79e96b9 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -68,6 +68,10 @@ if (${Torch_FOUND}) endif() #add_subdirectory(nbody) +if(UNIX AND NOT APPLE) + set(EXTRA_LIBS dl stdc++fs) +endif() + set(CANOE_LIBRARY_${buildu} "canoe_${buildl}" "athenapp_${buildl}" @@ -103,6 +107,7 @@ set(CANOE_LIBRARY_${buildu} ${TORCH_LIBRARY} ${TORCH_CPU_LIBRARY} ${C10_LIBRARY} + ${EXTRA_LIBS} CACHE INTERNAL "canoe library") set(CANOE_INCLUDE_DIR diff --git a/src/diagnostics/todo/convective_heat_flx.cpp b/src/diagnostics/todo/convective_heat_flx.cpp index 05ccd79e..b6ba7412 100644 --- a/src/diagnostics/todo/convective_heat_flx.cpp +++ b/src/diagnostics/todo/convective_heat_flx.cpp @@ -8,9 +8,9 @@ #include #include #include -#include // snap +#include #include // canoe diff --git a/src/snap/eos/shallow_yz_hydro.cpp b/src/snap/eos/shallow_yz_hydro.cpp index 1ce21ac3..f6286030 100644 --- a/src/snap/eos/shallow_yz_hydro.cpp +++ b/src/snap/eos/shallow_yz_hydro.cpp @@ -12,10 +12,10 @@ #include #include #include -#include // canoe #include +#include // exo3 #include diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 847365b3..b9747480 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -14,6 +14,7 @@ if (${Torch_FOUND}) setup_test(test_reconstruct) setup_test(test_riemann) setup_test(test_hydro) + setup_test(test_thermodynamics) endif() # athena tests @@ -49,7 +50,6 @@ if (${NVAPOR} EQUAL 2) setup_test(test_air_parcel) setup_test(test_radiation) setup_test(test_microwave_opacity) - setup_test(test_thermodynamics) setup_test(test_microphysics) elseif (${NCLOUD} EQUAL 5) setup_test(test_ammonium_hydrosulfide) From 4245eeaffd3e81d015e9647fe717b2770cb88058 Mon Sep 17 00:00:00 2001 From: stormy/cli Date: Mon, 8 Jul 2024 17:07:36 -0400 Subject: [PATCH 7/8] wip --- src/exchanger/exchanger.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/exchanger/exchanger.hpp b/src/exchanger/exchanger.hpp index 89210d6f..75b79beb 100644 --- a/src/exchanger/exchanger.hpp +++ b/src/exchanger/exchanger.hpp @@ -1,6 +1,10 @@ #ifndef SRC_EXCHANGER_EXCHANGER_HPP_ #define SRC_EXCHANGER_EXCHANGER_HPP_ +// C/C++ +#include +#include + // Athena++ #include #include From 4453d261ec82179115c24495094e6ec9cc792122 Mon Sep 17 00:00:00 2001 From: stormy/cli Date: Mon, 8 Jul 2024 17:36:49 -0400 Subject: [PATCH 8/8] wip --- CMakeLists.txt | 1 - tests/test_thermodynamics.inp | 2 -- 2 files changed, 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 10ede10f..3356e19b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -63,7 +63,6 @@ find_package(Cantera) ## 3. set up project system libraries ## message(STATUS "3. Set up system libraries") -message(STATUS ${CANTERA_FOUND}) if (NOT ${CANTERA_FOUND}) include(${CMAKE_SOURCE_DIR}/cmake/yamlpp.cmake) endif() diff --git a/tests/test_thermodynamics.inp b/tests/test_thermodynamics.inp index 1e981244..167f0089 100644 --- a/tests/test_thermodynamics.inp +++ b/tests/test_thermodynamics.inp @@ -32,8 +32,6 @@ vapor = H2O, NH3 cloud = H2O(c), NH3(c), H2O(p), NH3(p) -thermodynamics_config = jup-thermo.yaml - Rd = 3777. # mu = 2.3175 g/mol eps1 = 8.18 8.18 8.18 beta1 = 0. 24.845 24.845