Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

New thermodynamics #180

Merged
merged 15 commits into from
Dec 7, 2024
Prev Previous commit
Next Next commit
Revise thermodynamcis - bryan (#173)
- bryan case passes
- use cantera for thermodynamics calculation
- fix implicit correction for vapors using consistent diffusion
chengcli committed Dec 7, 2024

Verified

This commit was signed with the committer’s verified signature. The key has expired.
commit 075df46096ff763011cf7fe59424e9f1c23d523b
12 changes: 2 additions & 10 deletions cmake/athena.cmake
Original file line number Diff line number Diff line change
@@ -19,14 +19,6 @@ set(patch_command

set(PACKAGE_NAME athenapp)
set(REPO_URL "https://github.com/chengcli/athenapp")
set(REPO_TAG "12b11a6004827ba4732c248c021be55b2906ab82")
# set(REPO_TAG "12b11a6004827ba4732c248c021be55b2906ab82")
set(REPO_TAG "1750ec93a0969a9d572d8e0d6a094583b7007bb0")
add_package(${PACKAGE_NAME} ${REPO_URL} ${REPO_TAG} "${patch_command}" ON)

# FetchContent_Declare( athenapp GIT_REPOSITORY
# https://github.com/chengcli/athenapp/ GIT_TAG snap-mods PATCH_COMMAND
# ${patch_command} UPDATE_DISCONNECTED TRUE) DOWNLOAD_EXTRACT_TIMESTAMP TRUE URL
# https://github.com/chengcli/athenapp/archive/refs/tags/v0.8.tar.gz)

# FetchContent_MakeAvailable(athenapp)

# include_directories(${athenapp_SOURCE_DIR})
1 change: 1 addition & 0 deletions cmake/examples/bryan.cmake
Original file line number Diff line number Diff line change
@@ -16,5 +16,6 @@ set(NPHASE_LEGACY 2)
set(NETCDF OFF)
set(PNETCDF ON)
set(MPI ON)
set(EQUATION_OF_STATE ideal_moist)
# set(RSOLVER lmars)
set(RSOLVER hllc_transform)
2 changes: 1 addition & 1 deletion configure.hpp.in
Original file line number Diff line number Diff line change
@@ -9,7 +9,7 @@ enum {
NPHASE_LEGACY = @NPHASE_LEGACY@,

// air-parcle data
//NVAPOR = @NVAPOR@,
NVAPOR = @NVAPOR@,
NCLOUD = @NCLOUD@,
NMASS = @NMASS@,
NCHEMISTRY = @NCHEMISTRY@,
86 changes: 48 additions & 38 deletions examples/2019-Li-snap/bryan.cpp
Original file line number Diff line number Diff line change
@@ -37,23 +37,27 @@

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

// special includes
// #include "bryan_vapor_functions.hpp"

int iH2O, iH2Oc;
Real p0, grav;

void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
AllocateUserOutputVariables(7);
AllocateUserOutputVariables(11);
SetUserOutputVariableName(0, "temp");
SetUserOutputVariableName(1, "theta");
SetUserOutputVariableName(2, "theta_v");
SetUserOutputVariableName(3, "mse");
SetUserOutputVariableName(4, "theta_e");
SetUserOutputVariableName(5, "rh");
SetUserOutputVariableName(6, "qtol");
SetUserOutputVariableName(1, "tempv");
SetUserOutputVariableName(2, "enthalpy");
SetUserOutputVariableName(3, "entropy");
SetUserOutputVariableName(4, "intEng");

SetUserOutputVariableName(5, "theta");
SetUserOutputVariableName(6, "theta_v");
SetUserOutputVariableName(7, "mse");

SetUserOutputVariableName(8, "rh");
SetUserOutputVariableName(9, "theta_e");
SetUserOutputVariableName(10, "qtol");
}

void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
@@ -64,32 +68,41 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
for (int j = js; j <= je; ++j)
for (int i = is; i <= ie; ++i) {
user_out_var(0, k, j, i) = pthermo->GetTemp(w.at(k, j, i));
user_out_var(1, k, j, i) = pthermo->PotentialTemp(w.at(k, j, i), p0);
user_out_var(1, k, j, i) =
user_out_var(0, k, j, i) * pthermo->RovRd(w.at(k, j, i));
user_out_var(2, k, j, i) = pthermo->GetEnthalpy(w.at(k, j, i));
user_out_var(3, k, j, i) = pthermo->GetEntropy(w.at(k, j, i));
user_out_var(4, k, j, i) = pthermo->GetInternalEnergy(w.at(k, j, i));

// theta
user_out_var(5, k, j, i) = potential_temp(pthermo, w.at(k, j, i), p0);
// theta_v
user_out_var(2, k, j, i) =
user_out_var(1, j, i) * pthermo->RovRd(w.at(k, j, i));
/* mse
user_out_var(3, k, j, i) =
pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i);
user_out_var(6, k, j, i) =
user_out_var(5, k, j, i) * pthermo->RovRd(w.at(k, j, i));

// mse
user_out_var(7, k, j, i) =
moist_static_energy(pthermo, w.at(k, j, i), grav * pcoord->x1v(i));

// rh
user_out_var(8, k, j, i) = relative_humidity(pthermo, w.at(k, j, i))[1];
// theta_e
user_out_var(4, k, j, i) =
pthermo->EquivalentPotentialTemp(this, p0, iH2O, k, j, i);
user_out_var(5, k, j, i) =
pthermo->RelativeHumidity(this, iH2O, k, j, i);*/
user_out_var(9, k, j, i) = equivalent_potential_temp(
pthermo, w.at(k, j, i), user_out_var(8, k, j, i), p0);

// total mixing ratio
user_out_var(6, k, j, i) = w(iH2O, k, j, i) + w(iH2Oc, k, j, i);
user_out_var(10, k, j, i) = w(iH2O, k, j, i) + w(iH2Oc, k, j, i);
}
}

void Mesh::InitUserMeshData(ParameterInput *pin) {
auto pindex = IndexMap::GetInstance();

auto pthermo = Thermodynamics::GetInstance();
grav = -pin->GetReal("hydro", "grav_acc1");
p0 = pin->GetReal("problem", "p0");

// index
iH2O = pindex->GetVaporId("H2O");
iH2Oc = pindex->GetCloudId("H2O(c)");
iH2O = pthermo->SpeciesIndex("H2O");
iH2Oc = pthermo->SpeciesIndex("H2O(l)");
}

void MeshBlock::ProblemGenerator(ParameterInput *pin) {
@@ -107,14 +120,12 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
Real qt = pin->GetReal("problem", "qt");

// construct a reversible adiabat
std::vector<Real> yfrac({1. - qt, qt, 0.});
for (int k = ks; k <= ke; ++k)
for (int j = js; j <= je; ++j) {
w(IDN, k, j, is) = Ts;
w(IPR, k, j, is) = Ps;
w(iH2O, k, j, is) = qt;
w(iH2Oc, k, j, is) = 0.;
pthermo->SetMassFractions<Real>(yfrac.data());
pthermo->EquilibrateTP(Ts, Ps);

pthermo->SetPrimitive(w.at(k, j, is));
// half a grid to cell center
pthermo->Extrapolate_inplace(pcoord->dx1f(is) / 2., "reversible", grav);

@@ -135,21 +146,20 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
if (L < 1.) {
pthermo->SetPrimitive(w.at(k, j, i));
Real temp = pthermo->GetTemp();
Real pres = pthermo->GetPres();

Real temp_v = temp * pthermo->RovRd();
temp_v *= 1. + dT * sqr(cos(M_PI * L / 2.)) / 300.;

int err = root(temp, temp + dT, 1.E-8, &temp,
[&pthermo, temp_v](Real temp) {
pthermo->SetTemperature(temp);
pthermo->EquilibrateTP();
Real rovrd = pthermo->RovRd();
return temp * rovrd - temp_v;
int err = root(temp, temp + dT * Ts / 300., 1.E-8, &temp,
[&pthermo, temp_v, pres](Real temp) {
pthermo->EquilibrateTP(temp, pres);
return temp * pthermo->RovRd() - temp_v;
});

// if (err) throw RuntimeError("pgen", "TVSolver doesn't converge");

pthermo->SetTemperature(temp);
pthermo->EquilibrateTP();
pthermo->EquilibrateTP(temp, pres);
pthermo->GetPrimitive(w.at(k, j, i));
}
}
21 changes: 4 additions & 17 deletions examples/2019-Li-snap/bryan.inp
Original file line number Diff line number Diff line change
@@ -29,7 +29,7 @@ xorder = 5 # horizontal reconstruction order
integrator = rk3 # integration method

<mesh>
nx1 = 200. # Number of zones in X1-direction
nx1 = 200 # Number of zones in X1-direction
x1min = 0. # minimum value of X1
x1max = 10.E3 # maximum value of X1
ix1_bc = reflecting # inner-X1 boundary flag
@@ -56,23 +56,9 @@ nx3 = 1
grav_acc1 = -9.81
gamma = 1.4 # gamma = C_p/C_v

<species>
vapor = H2O
cloud = H2O(c)

<thermodynamics>
Rd = 287.
eps1 = 0.621 0.621
rcp1 = 1. 0.
beta1 = 0. 24.845
Ttriple1 = 273.16
Ptriple1 = 611.7

sa.max_iter = 2
sa.relax = 1.0
sa.ftol = 1.e-3

<problem>
thermodynamics_config = earth-moist.yaml

dT = 2.0
xc = 10.E3
zc = 2.E3
@@ -81,3 +67,4 @@ zr = 2.E3
p0 = 1.E5
Ts = 289.85
qt = 0.0196
#qt = 0.0
Original file line number Diff line number Diff line change
@@ -10,25 +10,28 @@ species:
composition: {O: 0.42, N: 1.56, Ar: 0.01}
thermo:
model: constant-cp
T0: 273.16
cp0: 29.1 J/mol/K

- name: H2O
composition: {H: 2, O: 1}
thermo:
model: constant-cp
cp0: 37.4 J/mol/K
T0: 273.16
cp0: 21.1 J/mol/K

- name: H2O(l)
composition: {H: 2, O: 1}
thermo:
model: constant-cp
cp0: 75.3 J/mol/K
h0: -28.5 kJ/mol
T0: 273.16
cp0: 62.556 J/mol/K
h0: -45.103 kJ/mol # h0 = - (beta - delta) * R * T0
equation-of-state:
model: constant-volume
density: 1.0 g/cm^3
molar-volume: 0.0 cm^3/mol

reactions:
- equation: H2O <=> H2O(l)
type: nucleation
rate-constant: {formula: ideal, T3: 273.16, P3: 611.7, beta: 24.845, delta: 4.986, minT: 273.16}
rate-constant: {formula: ideal, T3: 273.16, P3: 611.7, beta: 24.845, delta: 4.986, minT: 200.}
14 changes: 9 additions & 5 deletions examples/2019-Li-snap/jup-thermo.yaml
Original file line number Diff line number Diff line change
@@ -9,6 +9,10 @@
# h = h0 + cp0 * (t - t0)
# s = s0 + cp0 * ln(t/t0)

description: |-
This input file describes the properties of species
and their thermodynamics relations in Jupiter's atmosphere.

phases:
- name: gas
thermo: ideal-moist
@@ -48,7 +52,7 @@ species:
h0: -28.5 kJ/mol
equation-of-state:
model: constant-volume
density: 1.0 g/cm^3
molar-volume: 0.0 cm^3/mol

- name: NH3(l)
composition: {N: 1, H: 3}
@@ -58,7 +62,7 @@ species:
h0: -45.9 kJ/mol
equation-of-state:
model: constant-volume
density: 1.0 g/cm^3
molar-volume: 0.0 cm^3/mol

- name: H2O(s)
composition: {H: 2, O: 1}
@@ -68,7 +72,7 @@ species:
h0: -28.5 kJ/mol
equation-of-state:
model: constant-volume
density: 0.917 g/cm^3
molar-volume: 0.0 cm^3/mol

- name: NH3(s)
composition: {N: 1, H: 3}
@@ -78,7 +82,7 @@ species:
h0: -45.9 kJ/mol
equation-of-state:
model: constant-volume
density: 0.917 g/cm^3
molar-volume: 0.0 cm^3/mol

- name: NH4SH(s)
composition: {N: 1, H: 5, S: 1}
@@ -88,7 +92,7 @@ species:
h0: -105.9 kJ/mol
equation-of-state:
model: constant-volume
density: 1.5 g/cm^3
molar-volume: 0.0 cm^3/mol

reactions:
- equation: H2O <=> H2O(l)
4 changes: 2 additions & 2 deletions examples/2019-Li-snap/robert.cpp
Original file line number Diff line number Diff line change
@@ -27,7 +27,7 @@
#include <impl.hpp>

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

// @sect3{Preamble}
@@ -53,7 +53,7 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
for (int j = js; j <= je; ++j)
for (int i = is; i <= ie; ++i) {
user_out_var(0, k, j, i) = pthermo->GetTemp(w.at(k, j, i));
user_out_var(1, k, j, i) = pthermo->PotentialTemp(w.at(k, j, i), p0);
user_out_var(1, k, j, i) = potential_temp(pthermo, w.at(k, j, i), p0);
}
}

14 changes: 7 additions & 7 deletions examples/2019-Li-snap/straka.cpp
Original file line number Diff line number Diff line change
@@ -48,7 +48,7 @@
// Finally, the Thermodynamics class works with thermodynamic aspects of the
// problem such as the temperature, potential temperature, condensation of
// vapor, etc.
#include <snap/stride_iterator.hpp>
#include <snap/thermodynamics/atm_thermodynamics.hpp>
#include <snap/thermodynamics/thermodynamics.hpp>

// Functions in the math library are protected by a specific namespace because
@@ -103,7 +103,7 @@ void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
// <code>w</code> that stores density, pressure, and velocities at each
// grid.
user_out_var(0, j, i) = pthermo->GetTemp(w.at(k, j, i));
user_out_var(1, j, i) = pthermo->PotentialTemp(w.at(k, j, i), p0);
user_out_var(1, j, i) = potential_temp(pthermo, w.at(k, j, i), p0);
}
}

@@ -141,15 +141,15 @@ void Diffusion(MeshBlock *pmb, Real const time, Real const dt,
// the Thermodynamics class to calculate temperature and potential
// temperature.
Real temp = pthermo->GetTemp(w.at(pmb->ks, j, i));
Real theta = pthermo->PotentialTemp(w.at(pmb->ks, j, i), p0);
Real theta = potential_temp(pthermo, w.at(pmb->ks, j, i), p0);

// The thermal diffusion is applied to the potential temperature field,
// which is not exactly correct. But this is the setting of the test
// program.
Real theta_ip1_j = pthermo->PotentialTemp(w.at(k, j + 1, i), p0);
Real theta_im1_j = pthermo->PotentialTemp(w.at(k, j - 1, i), p0);
Real theta_i_jp1 = pthermo->PotentialTemp(w.at(k, j, i + 1), p0);
Real theta_i_jm1 = pthermo->PotentialTemp(w.at(k, j, i - 1), p0);
Real theta_ip1_j = potential_temp(pthermo, w.at(k, j + 1, i), p0);
Real theta_im1_j = potential_temp(pthermo, w.at(k, j - 1, i), p0);
Real theta_i_jp1 = potential_temp(pthermo, w.at(k, j, i + 1), p0);
Real theta_i_jm1 = potential_temp(pthermo, w.at(k, j, i - 1), p0);

// Add viscous dissipation to the velocities. Now you encounter another
// variable called <code>u</code>. <code>u</code> stands for `conserved
2 changes: 1 addition & 1 deletion patches/24.time_integrator.patch
Original file line number Diff line number Diff line change
@@ -169,7 +169,7 @@ index 2bfc5178..173b2566 100644
ave_wghts[2] = stage_wghts[stage-1].gamma_3;
if (ave_wghts[0] == 0.0 && ave_wghts[1] == 1.0 && ave_wghts[2] == 0.0) {
ps->s.SwapAthenaArray(ps->s1);
+ pmb->pimpl->MapScalarsConserved(ps->s);
+ // pmb->pimpl->MapScalarsConserved(ps->s);
} else {
pmb->WeightedAve(ps->s, ps->s1, ps->s2, ps->s0, ps->s_fl_div, ave_wghts);
}
Loading