Skip to content

Commit

Permalink
Fix all thermodynamic methods
Browse files Browse the repository at this point in the history
  • Loading branch information
@chengcli
chengcli committed Mar 30, 2024
1 parent 67fd23e commit 804f230
Show file tree
Hide file tree
Showing 14 changed files with 46 additions and 82 deletions.
15 changes: 6 additions & 9 deletions examples/2018-Li-harp/giants_re.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,14 +129,12 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {

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

// extrapolate down to where var is
pthermo->Extrapolate(&air, log(P0 / air.w[IPR]),
Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, log(P0 / air.w[IPR]), "dry");

// make up for the difference
Ts += T0 - air.w[IDN];
Expand All @@ -163,15 +161,15 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
for (; i <= ie; ++i) {
AirParcelHelper::distribute_to_primitive(this, k, j, i, air);

pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, -dlnp, "dry");
if (air.w[IDN] < Tmin) break;
}

// Replace adiabatic atmosphere with isothermal atmosphere if temperature
// is too low
pthermo->Extrapolate(&air, dlnp, Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, dlnp, "dry");
for (; i <= ie; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::Isothermal);
pthermo->Extrapolate(&air, -dlnp, "isothermal");
AirParcelHelper::distribute_to_primitive(this, k, j, i, air);
}
}
Expand All @@ -193,8 +191,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
air.ToMoleFraction();

for (int i = ibegin; i < iend; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat,
0., adlnTdlnP);
pthermo->Extrapolate(&air, -dlnp, "dry", 0., adlnTdlnP);
AirParcelHelper::distribute_to_primitive(this, k, j, i + 1, air);
}
}
Expand Down
9 changes: 3 additions & 6 deletions examples/2019-Li-snap/sedimentation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,8 +134,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 Down Expand Up @@ -165,16 +164,14 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
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);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav);
}

// 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);
}
}

Expand Down
3 changes: 1 addition & 2 deletions examples/2023-Chen-exo3/hot_jupiter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -244,8 +244,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
int i = is;
for (; i <= ie; ++i) {
AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::Isothermal, grav, 0.001);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav, 0.001);
// add noise
air.w[IVY] = 10. * distribution(generator);
air.w[IVZ] = 10. * distribution(generator);
Expand Down
6 changes: 2 additions & 4 deletions examples/2023-Chen-exo3/hs94.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -277,8 +277,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
if (pcoord->x1v(i) - Rp > z_iso) break;

AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::DryAdiabat, grav, 0.001);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav, 0.001);
// add noise
air.w[IVY] = 10. * distribution(generator);
air.w[IVZ] = 10. * distribution(generator);
Expand All @@ -287,8 +286,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
// construct isothermal atmosphere
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);
}
}
}
Expand Down
6 changes: 2 additions & 4 deletions examples/2023-Chen-exo3/hs94_polar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,17 +174,15 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
if (pcoord->x1v(i) - Rp > z_iso) break;

pimpl->DistributeToConserved(air, k, j, i);
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::DryAdiabat, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav);
// add noise
air.w[IVX] = 0.01 * distribution(generator);
}

// construct isothermal atmosphere
for (; i <= ie; ++i) {
pimpl->DistributeToConserved(air, k, j, i);
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::Isothermal, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav);
}
}
}
Expand Down
6 changes: 2 additions & 4 deletions examples/2023-Chen-exo3/polar_dry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -224,8 +224,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
// if (pcoord->x1v(i) - Rp > z_iso) break;

// pimpl->DistributeToConserved(air, k, j, i);
// pthermo->Extrapolate(&air, pcoord->dx1f(i),
// Thermodynamics::Method::DryAdiabat, grav,
// pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav,
// 0.001);
// // add noise
// // air.w[IVX] = 0.01 * distribution(generator);
Expand All @@ -234,8 +233,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
// construct isothermal atmosphere
for (; i <= ie; ++i) {
AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::Isothermal, grav, 0.001);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "isothermal", grav, 0.001);
// add noise
air.w[IVY] = 0.00001 * distribution(generator);
air.w[IVZ] = 0.00001 * distribution(generator);
Expand Down
3 changes: 1 addition & 2 deletions examples/2023-Chen-exo3/test_adiabat.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,8 +66,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
air.w[IDN] = Ts;
for (int i = is; i <= ie; ++i) {
AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::DryAdiabat, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav);
}
}
}
17 changes: 7 additions & 10 deletions examples/2023-Li-saturn-vla/saturn_radio.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,14 +181,12 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {

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

// extrapolate down to where var is
pthermo->Extrapolate(&air, log(P0 / air.w[IPR]),
Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, log(P0 / air.w[IPR]), "dry");

// make up for the difference
Ts += T0 - air.w[IDN];
Expand Down Expand Up @@ -220,16 +218,16 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {

AirParcelHelper::distribute_to_primitive(this, k, j, i, air);

pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, -dlnp, "dry");

if (air.w[IDN] < Tmin) break;
}

// Replace adiabatic atmosphere with isothermal atmosphere if temperature
// is too low
pthermo->Extrapolate(&air, dlnp, Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, dlnp, "dry");
for (; i <= ie; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::Isothermal);
pthermo->Extrapolate(&air, -dlnp, "isothermal");
AirParcelHelper::distribute_to_primitive(this, k, j, i, air);
}
}
Expand All @@ -252,8 +250,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
air.ToMoleFraction();

for (int i = ibegin; i < iend; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat,
0., adlnTdlnP);
pthermo->Extrapolate(&air, -dlnp, "dry", 0., adlnTdlnP);
AirParcelHelper::distribute_to_primitive(this, k, j, i + 1, air);
}
}
Expand All @@ -272,7 +269,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
air.ToMoleFraction();

for (int i = ibegin; i < iend; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, -dlnp, "dry");
air.w[iNH3] += adlnNH3dlnP * air.w[iNH3] * dlnp;
auto rates = pthermo->TryEquilibriumTP_VaporCloud(air, iNH3);
air.w[iNH3] += rates[0];
Expand Down
17 changes: 7 additions & 10 deletions examples/2023-Li-uranus/uranus_mwr.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,14 +180,12 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {

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

// extrapolate down to where var is
pthermo->Extrapolate(&air, log(P0 / air.w[IPR]),
Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, log(P0 / air.w[IPR]), "dry");

// make up for the difference
Ts += T0 - air.w[IDN];
Expand Down Expand Up @@ -220,16 +218,16 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {

AirParcelHelper::distribute_to_primitive(this, k, j, i, air);

pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, -dlnp, "dry");

if (air.w[IDN] < Tmin) break;
}

// Replace adiabatic atmosphere with isothermal atmosphere if temperature
// is too low
pthermo->Extrapolate(&air, dlnp, Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, dlnp, "dry");
for (; i <= ie; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::Isothermal);
pthermo->Extrapolate(&air, -dlnp, "isothermal");
AirParcelHelper::distribute_to_primitive(this, k, j, i, air);
}
}
Expand All @@ -252,8 +250,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
air.ToMoleFraction();

for (int i = ibegin; i < iend; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat,
0., adlnTdlnP);
pthermo->Extrapolate(&air, -dlnp, "dry", 0., adlnTdlnP);
AirParcelHelper::distribute_to_primitive(this, k, j, i + 1, air);
}
}
Expand All @@ -272,7 +269,7 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
air.ToMoleFraction();

for (int i = ibegin; i < iend; ++i) {
pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat);
pthermo->Extrapolate(&air, -dlnp, "dry");
air.w[iNH3] += adlnNH3dlnP * air.w[iNH3] * dlnp;
auto rates = pthermo->TryEquilibriumTP_VaporCloud(air, iNH3);
air.w[iNH3] += rates[0];
Expand Down
12 changes: 4 additions & 8 deletions examples/2024-FDing-jupiter-rt/jupiter1d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,8 +126,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 @@ -153,23 +152,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::DryAdiabat, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "dry", grav);
}

// 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);
}
}
}
6 changes: 2 additions & 4 deletions examples/2024-Li-plume/plume.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,15 +152,13 @@ void MeshBlock::ProblemGenerator(ParameterInput *pin) {
air.c[iH2Oc] = 0.;

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

for (int i = is; i <= ie; ++i) {
// add noise
air.w[IVX] = 0.01 * (1. * rand() / RAND_MAX - 0.5);
AirParcelHelper::distribute_to_conserved(this, k, j, i, air);
pthermo->Extrapolate(&air, pcoord->dx1f(i),
Thermodynamics::Method::PseudoAdiabat, grav);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav);
}
}
}
9 changes: 3 additions & 6 deletions examples/2024-XZhang-cloud-rt/hjupiter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,23 +181,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);
pthermo->Extrapolate(&air, pcoord->dx1f(i), "pseudo", grav);
}

// 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);
}
}
}
13 changes: 4 additions & 9 deletions examples/2024-XZhang-cloud-rt/hywater.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,8 +146,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 @@ -172,25 +171,21 @@ 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;
air.w[IVX] = 0.1 * sin(2. * M_PI * rand() / RAND_MAX);
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);
}

peos->ConservedToPrimitive(phydro->u, phydro->w, pfield->b, phydro->w,
Expand Down
Loading

0 comments on commit 804f230

Please sign in to comment.