Skip to content

Commit 56a7ec5

Browse files
BenWibkingBenWibking
committed
update HydrostaticAtmosphere test
1 parent 9504b66 commit 56a7ec5

2 files changed

Lines changed: 154 additions & 88 deletions

File tree

‎inputs/HydrostaticAtmosphere.in‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,27 @@
11
# Hydrostatic exponential atmosphere density floor unit test inputs
22

3+
# *****************************************************************
4+
# Problem size and geometry
5+
# *****************************************************************
6+
geometry.prob_lo = 0.0 0.0 0.0
7+
geometry.prob_hi = 1.0 1.0 1.0
8+
quokka.bc = ext_dir periodic periodic
9+
10+
# *****************************************************************
11+
# VERBOSITY
12+
# *****************************************************************
13+
amr.v = 0 # verbosity in Amr
14+
15+
# *****************************************************************
16+
# Resolution and refinement
17+
# *****************************************************************
18+
amr.n_cell = 8 8 8
19+
amr.max_level = 0 # number of levels = max_level + 1
20+
amr.blocking_factor = 1 # grid size must be divisible by this
21+
22+
do_reflux = 0
23+
do_subcycle = 0
24+
325
density_floor = 1.0
426
atmosphere_scale_height = 0.25
527
density_floor_expr = "1.0e-2*base_density_floor*exp(-x/0.25)"

‎src/problems/HydrostaticAtmosphere/testHydrostaticAtmosphere.cpp‎

Lines changed: 132 additions & 88 deletions
Original file line numberDiff line numberDiff line change
@@ -9,22 +9,21 @@
99

1010
#include <cmath>
1111

12-
#include "AMReX_Geometry.H"
13-
#include "AMReX_Math.H"
14-
#include "AMReX_ParallelDescriptor.H"
1512
#include "AMReX_ParmParse.H"
16-
#include "AMReX_Parser.H"
13+
#include "AMReX_Print.H"
1714
#include "AMReX_REAL.H"
18-
#include "AMReX_Reduce.H"
1915

20-
#include "eos.H"
21-
#include "extern_parameters.H"
16+
#include "QuokkaSimulation.hpp"
2217
#include "fundamental_constants.H"
2318
#include "hydro/hydro_system.hpp"
2419

2520
struct HydrostaticAtmosphereProblem {
2621
};
2722

23+
template <> struct SimulationData<HydrostaticAtmosphereProblem> {
24+
amrex::Real atmosphere_scale_height = NAN;
25+
};
26+
2827
template <> struct quokka::EOS_Traits<HydrostaticAtmosphereProblem> {
2928
static constexpr double gamma = 5. / 3.;
3029
static constexpr double mean_molecular_weight = C::m_u;
@@ -43,25 +42,138 @@ template <> struct Physics_Traits<HydrostaticAtmosphereProblem> {
4342
static constexpr UnitSystem unit_system = UnitSystem::CGS;
4443
};
4544

46-
auto problem_main() -> int
45+
constexpr amrex::Real kTgasInit = 1.0;
46+
constexpr amrex::Real kRhoInitFactor = 5.0e-3;
47+
48+
template <> void QuokkaSimulation<HydrostaticAtmosphereProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)
49+
{
50+
// extract variables required from the geom object
51+
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const dx = grid_elem.dx_;
52+
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const prob_lo = grid_elem.prob_lo_;
53+
const amrex::Box &indexRange = grid_elem.indexRange_;
54+
const amrex::Array4<double> &state_cc = grid_elem.array_;
55+
const int ncomp_cc = Physics_Indices<HydrostaticAtmosphereProblem>::nvarTotal_cc;
56+
57+
amrex::Real const base_density_floor = densityFloor_;
58+
amrex::Real const scale_height = userData_.atmosphere_scale_height;
59+
60+
// loop over the grid and set the initial condition
61+
amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
62+
amrex::Real const x = prob_lo[0] + (i + static_cast<amrex::Real>(0.5)) * dx[0];
63+
amrex::Real const rho_atm = base_density_floor * std::exp(-x / scale_height);
64+
amrex::Real const rho_init = kRhoInitFactor * rho_atm;
65+
amrex::Real const Eint_init = quokka::EOS<HydrostaticAtmosphereProblem>::ComputeEintFromTgas(rho_init, kTgasInit);
66+
67+
for (int n = 0; n < ncomp_cc; ++n) {
68+
state_cc(i, j, k, n) = 0.;
69+
}
70+
71+
state_cc(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::density_index) = rho_init;
72+
state_cc(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::energy_index) = Eint_init;
73+
state_cc(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::internalEnergy_index) = Eint_init;
74+
});
75+
}
76+
77+
template <>
78+
void QuokkaSimulation<HydrostaticAtmosphereProblem>::computeReferenceSolution(amrex::MultiFab &ref, amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &dx,
79+
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const &prob_lo)
4780
{
48-
init_extern_parameters();
49-
amrex::Real small_temp = 1.0e-10;
50-
amrex::Real small_dens = 1.0e-100;
51-
eos_init(small_temp, small_dens);
81+
amrex::Real const base_density_floor = densityFloor_;
82+
amrex::Real const scale_height = userData_.atmosphere_scale_height;
83+
const int ncomp_cc = ref.nComp();
84+
85+
if (useDensityFloorParser_) {
86+
auto const density_floor_parser = densityFloorParserExe_.value();
87+
for (amrex::MFIter iter(ref); iter.isValid(); ++iter) {
88+
const amrex::Box &indexRange = iter.validbox();
89+
auto const &state_ref = ref.array(iter);
90+
91+
amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
92+
amrex::Real const x = prob_lo[0] + (i + static_cast<amrex::Real>(0.5)) * dx[0];
93+
#if (AMREX_SPACEDIM >= 2)
94+
amrex::Real const y = prob_lo[1] + (j + static_cast<amrex::Real>(0.5)) * dx[1];
95+
#else
96+
amrex::Real const y = 0.0;
97+
#endif
98+
#if (AMREX_SPACEDIM == 3)
99+
amrex::Real const z = prob_lo[2] + (k + static_cast<amrex::Real>(0.5)) * dx[2];
100+
#else
101+
amrex::Real const z = 0.0;
102+
#endif
103+
amrex::Real const rho_atm = base_density_floor * std::exp(-x / scale_height);
104+
amrex::Real const rho_floor = density_floor_parser(x, y, z, base_density_floor);
105+
amrex::Real const rho_init = kRhoInitFactor * rho_atm;
106+
amrex::Real const Eint_init = quokka::EOS<HydrostaticAtmosphereProblem>::ComputeEintFromTgas(rho_init, kTgasInit);
107+
108+
for (int n = 0; n < ncomp_cc; ++n) {
109+
state_ref(i, j, k, n) = 0.;
110+
}
111+
112+
state_ref(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::density_index) = rho_floor;
113+
state_ref(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::energy_index) = Eint_init;
114+
state_ref(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::internalEnergy_index) = Eint_init;
115+
});
116+
}
117+
} else {
118+
auto const density_floor_func = [this] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z,
119+
amrex::Real base_floor) -> amrex::Real {
120+
return densityFloor(x, y, z, base_floor);
121+
};
122+
for (amrex::MFIter iter(ref); iter.isValid(); ++iter) {
123+
const amrex::Box &indexRange = iter.validbox();
124+
auto const &state_ref = ref.array(iter);
125+
126+
amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
127+
amrex::Real const x = prob_lo[0] + (i + static_cast<amrex::Real>(0.5)) * dx[0];
128+
#if (AMREX_SPACEDIM >= 2)
129+
amrex::Real const y = prob_lo[1] + (j + static_cast<amrex::Real>(0.5)) * dx[1];
130+
#else
131+
amrex::Real const y = 0.0;
132+
#endif
133+
#if (AMREX_SPACEDIM == 3)
134+
amrex::Real const z = prob_lo[2] + (k + static_cast<amrex::Real>(0.5)) * dx[2];
135+
#else
136+
amrex::Real const z = 0.0;
137+
#endif
138+
amrex::Real const rho_atm = base_density_floor * std::exp(-x / scale_height);
139+
amrex::Real const rho_floor = density_floor_func(x, y, z, base_density_floor);
140+
amrex::Real const rho_init = kRhoInitFactor * rho_atm;
141+
amrex::Real const Eint_init = quokka::EOS<HydrostaticAtmosphereProblem>::ComputeEintFromTgas(rho_init, kTgasInit);
142+
143+
for (int n = 0; n < ncomp_cc; ++n) {
144+
state_ref(i, j, k, n) = 0.;
145+
}
146+
147+
state_ref(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::density_index) = rho_floor;
148+
state_ref(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::energy_index) = Eint_init;
149+
state_ref(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::internalEnergy_index) = Eint_init;
150+
});
151+
}
152+
}
153+
}
52154

155+
auto problem_main() -> int
156+
{
53157
amrex::ParmParse const pp;
54158
amrex::Real base_density_floor = 0.0;
55159
if (pp.query("density_floor", base_density_floor) == 0) {
56160
amrex::Print() << "density_floor must be set for HydrostaticAtmosphere test.\n";
57161
return 1;
58162
}
163+
if (base_density_floor <= 0.0) {
164+
amrex::Print() << "density_floor must be positive for HydrostaticAtmosphere test.\n";
165+
return 1;
166+
}
59167

60168
amrex::Real scale_height = 0.0;
61169
if (pp.query("atmosphere_scale_height", scale_height) == 0) {
62170
amrex::Print() << "atmosphere_scale_height must be set for HydrostaticAtmosphere test.\n";
63171
return 1;
64172
}
173+
if (scale_height <= 0.0) {
174+
amrex::Print() << "atmosphere_scale_height must be positive for HydrostaticAtmosphere test.\n";
175+
return 1;
176+
}
65177

66178
std::string density_floor_expr;
67179
pp.query("density_floor_expr", density_floor_expr);
@@ -70,86 +182,18 @@ auto problem_main() -> int
70182
return 1;
71183
}
72184

73-
amrex::Parser parser(density_floor_expr);
74-
parser.registerVariables({"x", "y", "z", "base_density_floor"});
75-
auto const parser_exe = parser.compile<4>();
76-
77-
constexpr int nx = 4;
78-
constexpr int ny = 1;
79-
constexpr int nz = 1;
80-
amrex::IntVect const dom_lo(AMREX_D_DECL(0, 0, 0));
81-
amrex::IntVect const dom_hi(AMREX_D_DECL(nx - 1, ny - 1, nz - 1));
82-
amrex::Box const domain(dom_lo, dom_hi);
83-
amrex::RealBox const real_box({AMREX_D_DECL(0.0, 0.0, 0.0)}, {AMREX_D_DECL(1.0, 1.0, 1.0)});
84-
amrex::Array<int, AMREX_SPACEDIM> const is_periodic{AMREX_D_DECL(0, 0, 0)};
85-
amrex::Geometry const geom(domain, &real_box, amrex::CoordSys::cartesian, is_periodic.data());
86-
87-
amrex::BoxArray ba(domain);
88-
ba.maxSize(domain.size());
89-
amrex::DistributionMapping const dm(ba);
90-
int const ncomp = Physics_Indices<HydrostaticAtmosphereProblem>::nvarTotal_cc;
91-
amrex::MultiFab state(ba, dm, ncomp, 0);
92-
93-
state.setVal(0.0);
94-
95-
amrex::Real const Tgas_init = 1.0;
96-
amrex::Real const rho_init_factor = 5.0e-3;
97-
auto const *const prob_lo = geom.ProbLo();
98-
auto const *const dx = geom.CellSize();
99-
for (amrex::MFIter mfi(state); mfi.isValid(); ++mfi) {
100-
auto const &arr = state.array(mfi);
101-
amrex::Box const &bx = mfi.validbox();
102-
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
103-
amrex::Real const x = prob_lo[0] + (static_cast<amrex::Real>(i) + 0.5) * dx[0];
104-
amrex::Real const rho_atm = base_density_floor * std::exp(-x / scale_height);
105-
amrex::Real const rho_init = rho_init_factor * rho_atm;
106-
amrex::Real const Eint_init = quokka::EOS<HydrostaticAtmosphereProblem>::ComputeEintFromTgas(rho_init, Tgas_init);
107-
arr(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::density_index) = rho_init;
108-
arr(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::energy_index) = Eint_init;
109-
arr(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::internalEnergy_index) = Eint_init;
110-
});
111-
}
112-
113-
auto const density_floor_func = [=] AMREX_GPU_HOST_DEVICE(amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real base_floor) -> amrex::Real {
114-
return parser_exe(x, y, z, base_floor);
115-
};
116-
117-
HydroSystem<HydrostaticAtmosphereProblem>::EnforceLimits(base_density_floor, 0.0, state, geom.data(), density_floor_func);
118-
119-
amrex::ReduceOps<amrex::ReduceOpMax> reduce_op;
120-
amrex::ReduceData<amrex::Real> reduce_data(reduce_op);
121-
122-
for (amrex::MFIter mfi(state); mfi.isValid(); ++mfi) {
123-
amrex::Box const &bx = mfi.validbox();
124-
auto const &data = state.array(mfi);
125-
126-
reduce_op.eval(bx, reduce_data, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept -> amrex::GpuTuple<amrex::Real> {
127-
amrex::Real const x = prob_lo[0] + (static_cast<amrex::Real>(i) + 0.5) * dx[0];
128-
#if (AMREX_SPACEDIM >= 2)
129-
amrex::Real const y = prob_lo[1] + (static_cast<amrex::Real>(j) + 0.5) * dx[1];
130-
#else
131-
amrex::Real const y = 0.0;
132-
#endif
133-
#if (AMREX_SPACEDIM == 3)
134-
amrex::Real const z = prob_lo[2] + (static_cast<amrex::Real>(k) + 0.5) * dx[2];
135-
#else
136-
amrex::Real const z = 0.0;
137-
#endif
138-
amrex::ignore_unused(y, z);
139-
amrex::Real const rho_atm = base_density_floor * std::exp(-x / scale_height);
140-
amrex::Real const expected = 1.0e-2 * rho_atm;
141-
amrex::Real const actual = data(i, j, k, HydroSystem<HydrostaticAtmosphereProblem>::density_index);
142-
return {amrex::Math::abs(actual - expected)};
143-
});
144-
}
185+
QuokkaSimulation<HydrostaticAtmosphereProblem> sim;
186+
sim.userData_.atmosphere_scale_height = scale_height;
187+
sim.plotfileInterval_ = -1;
145188

146-
auto [max_err] = reduce_data.value();
147-
amrex::ParallelDescriptor::ReduceRealMax(max_err);
189+
sim.setInitialConditions();
190+
sim.FixupState(0);
148191

192+
amrex::Real const error_norm = sim.computeErrorNorm(false);
149193
amrex::Real const tol = 1.0e-12;
150194
int status = 0;
151-
if (!(max_err <= tol)) {
152-
amrex::Print() << "Max density floor error = " << max_err << "\n";
195+
if (!(error_norm <= tol)) {
196+
amrex::Print() << "Density floor error norm = " << error_norm << "\n";
153197
status = 1;
154198
}
155199

0 commit comments

Comments
 (0)