Skip to content

Commit

Permalink
wip
Browse files Browse the repository at this point in the history
  • Loading branch information
@chengcli
chengcli committed Mar 29, 2024
1 parent eaeca49 commit 94f852e
Show file tree
Hide file tree
Showing 5 changed files with 26 additions and 29 deletions.
14 changes: 7 additions & 7 deletions examples/2019-Li-snap/bryan.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@
#include <climath/root.hpp>

// snap
#include <snap/thermodynamics/thermodynamics.hpp>
#include <snap/thermodynamics/atm_thermodynamics.hpp>

// special includes
#include "bryan_vapor_functions.hpp"
Expand Down Expand Up @@ -120,13 +120,11 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
air.c[iH2Oc] = 0.;

// half a grid to cell center
pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2.,
Thermodynamics::Method::ReversibleAdiabat, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav);

for (int i = is; i <= ie; ++i) {
AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::ReversibleAdiabat, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "reversible", grav);
}
}

Expand All @@ -141,7 +139,8 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
if (L < 1.) {
auto &&air = AirParcelHelper::gather_from_conserved(this, k, j, i);
air.ToMoleFraction();
Real temp_v = air.w[IDN] * pthermo->RovRd(air);
Real rovrd = get_rovrd(air, pthermo->GetMuRatio());
Real temp_v = air.w[IDN] * rovrd;
temp_v *= 1. + dT * sqr(cos(M_PI * L / 2.)) / 300.;

Real temp;
Expand All @@ -153,8 +152,9 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
pthermo->TryEquilibriumTP_VaporCloud(air, iH2O);
air.w[iH2O] += rates[0];
air.c[iH2Oc] += rates[1];
Real rovrd = get_rovrd(air, pthermo->GetMuRatio());

return temp * pthermo->RovRd(air) - temp_v;
return temp * rovrd - temp_v;
});

if (err) throw RuntimeError("pgen", "TVSolver doesn't converge");
Expand Down
19 changes: 8 additions & 11 deletions examples/2019-Li-snap/jupiter_crm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@
#include <climath/interpolation.h>

// snap
#include <snap/thermodynamics/thermodynamics.hpp>
#include <snap/thermodynamics/atm_thermodynamics.hpp>

// special includes
#include <special/giants_enroll_vapor_functions_v1.hpp>
Expand Down Expand Up @@ -87,7 +87,9 @@ void Forcing(MeshBlock *pmb, Real const time, Real const dt,
for (int i = is; i <= ie; ++i) {
auto &&air = AirParcelHelper::gather_from_primitive(pmb, k, j, i);

Real cv = pthermo->GetCvMass(air, 0);
air.ToMoleFraction();
Real cv =
get_cv_mass(air, 0, pthermo->GetRd(), pthermo->GetCvRatioMass());

if (w(IPR, k, j, i) < prad) {
du(IEN, k, j, i) += dt * hrate * w(IDN, k, j, i) * cv *
Expand Down Expand Up @@ -159,8 +161,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {

// stop at just above P0
for (int i = is; i <= ie; ++i) {
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::PseudoAdiabat, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav);
if (air.w[IPR] < P0) break;
}

Expand All @@ -186,24 +187,20 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
air.w[IDN] = Ts;

// half a grid to cell center
pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2.,
Thermodynamics::Method::ReversibleAdiabat, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(is) / 2., "reversible", grav);

int i = is;
for (; i <= ie; ++i) {
if (air.w[IDN] < Tmin) break;
AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::PseudoAdiabat, grav,
1.e-5);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav, 1.e-5);
}

// Replace adiabatic atmosphere with isothermal atmosphere if temperature
// is too low
for (; i <= ie; ++i) {
AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::Isothermal, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav);
}
}
}
4 changes: 0 additions & 4 deletions examples/2024-CLi-juno-mwr/juno_mwr.inp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,10 +94,6 @@ Ptriple2 = 6060.
radiation_config = juno_mwr.yaml
outdir = (0,) (15,) (30,) (45,)

#<inversion>
#tasks = JunoProfileInversion
#control_file = mwr_inversion.yaml

<problem>
use_temperature_dependent_cp = true
use_fletcher16_cirs = false
Expand Down
12 changes: 9 additions & 3 deletions src/snap/thermodynamics/thermodynamics.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,21 +138,24 @@ class Thermodynamics {
//! \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_[0]; }
Real const *GetCvRatioMole() const { return cv_ratio_mole_.data(); }

//! \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]; }

//! const pointer to cloud_index_set_
IndexSet const *GetCloudIndexSet() const { return &cloud_index_set_[0]; }
IndexSet const *GetCloudIndexSet() const { return cloud_index_set_.data(); }

//! Reference specific heat capacity [J/(kg K)] at constant volume
Real GetCvMassRef(int n) const {
Expand All @@ -164,12 +167,15 @@ class Thermodynamics {
//! \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_[0]; }
Real const *GetCpRatioMole() const { return cp_ratio_mole_.data(); }

Real GetCpMassRef(int n) const {
Real cpd = Rd_ * gammad_ref_ / (gammad_ref_ - 1.);
Expand Down
6 changes: 2 additions & 4 deletions tests/test_convective_adjustment.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,8 +254,7 @@ TEST_F(TestConvectiveAdjustment, RandomProfile) {
auto pthermo = Thermodynamics::GetInstance();

// half a grid to cell center
pthermo->Extrapolate(&air, pcoord->dx1f(pmb->is) / 2.,
Thermodynamics::Method::ReversibleAdiabat, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(pmb->is) / 2., "reversible", grav);

AirColumn ac(pmb->ncells1);

Expand All @@ -270,8 +269,7 @@ TEST_F(TestConvectiveAdjustment, RandomProfile) {
<< std::endl;

ac[i] = air;
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::PseudoAdiabat, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav);
}
outFile1.close();

Expand Down

0 comments on commit 94f852e

Please sign in to comment.