Skip to content

Commit 5b493c6

Browse files
mxchen39mxchen39
committed
add a DustDampingIteration test which runs both reference solution and numerical solution (the other 4 tests can be deleted later)
1 parent 83c783b commit 5b493c6

3 files changed

Lines changed: 362 additions & 0 deletions

File tree

‎inputs/DustDampingIteration.in‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
# *****************************************************************
2+
# Problem size and geometry
3+
# *****************************************************************
4+
geometry.prob_lo = 0.0 0.0 0.0
5+
geometry.prob_hi = 1.0 1.0 1.0
6+
quokka.bc = periodic periodic periodic
7+
8+
# *****************************************************************
9+
# VERBOSITY
10+
# *****************************************************************
11+
amr.v = 0 # verbosity in Amr
12+
13+
# *****************************************************************
14+
# Resolution and refinement
15+
# *****************************************************************
16+
amr.n_cell = 128 8 8
17+
amr.max_level = 0 # number of levels = max_level + 1
18+
19+
do_reflux = 0
20+
do_subcycle = 0
21+
suppress_output = 0
22+
23+
max_timesteps = 10000000
24+
25+
dust.omega = 1.0
26+
dust.enable_iter_stoptime = 1
27+
28+
# test A
29+
#stop_time = 10.0
30+
31+
# test B
32+
stop_time = 0.05
33+
34+
# test C
35+
#stop_time = 0.3

‎src/problems/DustDampingIteration/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
quokka_add_problem(JOB_NAME DustDampingIteration)

‎src/problems/DustDampingIteration/testDustDampingIteration.cpp‎

Lines changed: 326 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,326 @@
1+
/// \file testDustDampingIter.cpp
2+
/// \brief Defines a test problem for dust drag
3+
///
4+
5+
#include "QuokkaSimulation.hpp"
6+
#include "util/fextract.hpp"
7+
#include <cmath>
8+
#include <fmt/format.h>
9+
#include <fstream>
10+
#ifdef HAVE_PYTHON
11+
#include "util/matplotlibcpp.h"
12+
#endif
13+
14+
constexpr double rho_dust1 = 1.0;
15+
constexpr double rho_dust2 = 1.0;
16+
constexpr double P_INITIAL = 1.0;
17+
18+
struct DustDamping {
19+
};
20+
21+
template <> struct SimulationData<DustDamping> {
22+
std::vector<double> t_vec_;
23+
std::vector<double> v_gas_vec_;
24+
std::vector<double> v_dust1_vec_;
25+
std::vector<double> v_dust2_vec_;
26+
std::vector<double> E_gas_vec_;
27+
};
28+
29+
template <> struct quokka::EOS_Traits<DustDamping> {
30+
static constexpr double mean_molecular_weight = 1.0;
31+
static constexpr double gamma = 1.4;
32+
// static constexpr double cs_isothermal = 1.0; // only used when gamma = 1
33+
};
34+
35+
constexpr double rho = 1.0;
36+
constexpr double v0 = 1.0;
37+
constexpr double Egas0 = P_INITIAL / (quokka::EOS_Traits<DustDamping>::gamma - 1.0) + 0.5 * rho * v0 * v0;
38+
constexpr double Egas0_internal = P_INITIAL / (quokka::EOS_Traits<DustDamping>::gamma - 1.0);
39+
constexpr int numDustVars = Physics_NumVars::numDustVarsPerGroup;
40+
static constexpr amrex::GpuArray<amrex::Real, 2> dust_grain_radius = {0.02, 0.01};
41+
static constexpr amrex::GpuArray<amrex::Real, 2> dust_grain_density = {1.0, 1.0};
42+
static constexpr bool enable_supersonic_correction = true;
43+
44+
template <> struct Physics_Traits<DustDamping> {
45+
static constexpr bool is_self_gravity_enabled = false;
46+
static constexpr bool is_hydro_enabled = true;
47+
static constexpr int numMassScalars = 0; // number of mass scalars
48+
static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars
49+
static constexpr bool is_radiation_enabled = false;
50+
static constexpr bool is_dust_enabled = true;
51+
static constexpr int nDustGroups = 2; // number of dust groups
52+
static constexpr bool is_mhd_enabled = false;
53+
static constexpr int nGroups = 1; // number of radiation groups
54+
static constexpr UnitSystem unit_system = UnitSystem::CONSTANTS;
55+
static constexpr double boltzmann_constant = 1.0;
56+
static constexpr double gravitational_constant = 1.0;
57+
static constexpr double c_light = 1.0;
58+
static constexpr double radiation_constant = 1.0;
59+
};
60+
61+
template <>
62+
AMREX_GPU_HOST_DEVICE auto DustDrag<DustDamping>::ComputeReciprocalStoppingTime(
63+
amrex::Real rho_g, amrex::GpuArray<amrex::Real, Physics_Traits<DustDamping>::nDustGroups> rho_d,
64+
amrex::GpuArray<amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>, Physics_Traits<DustDamping>::nDustGroups + 1> vel, double cs)
65+
-> amrex::GpuArray<amrex::Real, Physics_Traits<DustDamping>::nDustGroups>
66+
{
67+
return ComputeReciprocalStoppingTimeHelper(rho_g, rho_d, vel, cs, dust_grain_radius, dust_grain_density, enable_supersonic_correction);
68+
}
69+
70+
template <> void QuokkaSimulation<DustDamping>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)
71+
{
72+
const amrex::Box &indexRange = grid_elem.indexRange_;
73+
const amrex::Array4<double> &state_cc = grid_elem.array_;
74+
75+
const auto vx0 = v0; // gas velocity
76+
const auto vx_dust1 = 2 * v0; // dust1 velocity
77+
const auto vx_dust2 = 10.0 * v0; // dust2 velocity
78+
79+
amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
80+
// for gas
81+
state_cc(i, j, k, HydroSystem<DustDamping>::density_index) = rho;
82+
state_cc(i, j, k, HydroSystem<DustDamping>::energy_index) = Egas0;
83+
state_cc(i, j, k, HydroSystem<DustDamping>::internalEnergy_index) = Egas0_internal;
84+
state_cc(i, j, k, HydroSystem<DustDamping>::x1Momentum_index) = rho * vx0;
85+
state_cc(i, j, k, HydroSystem<DustDamping>::x2Momentum_index) = 0.;
86+
state_cc(i, j, k, HydroSystem<DustDamping>::x3Momentum_index) = 0.;
87+
88+
// first-capture for CUDA
89+
const auto vx_dust1_local = vx_dust1;
90+
const auto vx_dust2_local = vx_dust2;
91+
92+
if constexpr (Physics_Traits<DustDamping>::is_dust_enabled) {
93+
// for dust1
94+
state_cc(i, j, k, HydroSystem<DustDamping>::dustDensity_index) = rho_dust1;
95+
state_cc(i, j, k, HydroSystem<DustDamping>::x1DustMomentum_index) = rho_dust1 * vx_dust1_local;
96+
state_cc(i, j, k, HydroSystem<DustDamping>::x2DustMomentum_index) = 0.;
97+
state_cc(i, j, k, HydroSystem<DustDamping>::x3DustMomentum_index) = 0.;
98+
// for dust2
99+
state_cc(i, j, k, HydroSystem<DustDamping>::dustDensity_index + numDustVars) = rho_dust2;
100+
state_cc(i, j, k, HydroSystem<DustDamping>::x1DustMomentum_index + numDustVars) = rho_dust2 * vx_dust2_local;
101+
state_cc(i, j, k, HydroSystem<DustDamping>::x2DustMomentum_index + numDustVars) = 0.;
102+
state_cc(i, j, k, HydroSystem<DustDamping>::x3DustMomentum_index + numDustVars) = 0.;
103+
}
104+
});
105+
}
106+
107+
template <> void QuokkaSimulation<DustDamping>::computeBeforeTimestep()
108+
{
109+
// extract initial physical quantities at t=0
110+
if (amrex::ParallelDescriptor::IOProcessor() && userData_.t_vec_.empty()) {
111+
auto [position, values] = fextract(state_new_cc_[0], Geom(0), 0, 0.5);
112+
113+
userData_.t_vec_.push_back(0.0); // initial time t=0
114+
115+
// extract physical quantities
116+
const double density = values.at(HydroSystem<DustDamping>::density_index)[0];
117+
const double momentum_x = values.at(HydroSystem<DustDamping>::x1Momentum_index)[0];
118+
const double Egas_total = values.at(HydroSystem<DustDamping>::energy_index)[0];
119+
120+
// store gas velocity
121+
const double v_gas = momentum_x / density;
122+
userData_.v_gas_vec_.push_back(v_gas);
123+
124+
// store gas total energy
125+
userData_.E_gas_vec_.push_back(Egas_total);
126+
127+
if constexpr (Physics_Traits<DustDamping>::is_dust_enabled) {
128+
// store dust1 velocity
129+
const double dust1_density = values.at(HydroSystem<DustDamping>::dustDensity_index)[0];
130+
const double dust1_momentum_x = values.at(HydroSystem<DustDamping>::x1DustMomentum_index)[0];
131+
const double v_dust1 = dust1_momentum_x / dust1_density;
132+
userData_.v_dust1_vec_.push_back(v_dust1);
133+
134+
// store dust2 velocity
135+
const double dust2_density = values.at(HydroSystem<DustDamping>::dustDensity_index + numDustVars)[0];
136+
const double dust2_momentum_x = values.at(HydroSystem<DustDamping>::x1DustMomentum_index + numDustVars)[0];
137+
const double v_dust2 = dust2_momentum_x / dust2_density;
138+
userData_.v_dust2_vec_.push_back(v_dust2);
139+
}
140+
}
141+
}
142+
143+
template <> void QuokkaSimulation<DustDamping>::computeAfterTimestep()
144+
{
145+
auto [position, values] = fextract(state_new_cc_[0], Geom(0), 0, 0.5);
146+
147+
if (amrex::ParallelDescriptor::IOProcessor()) {
148+
userData_.t_vec_.push_back(tNew_[0]); // store current time
149+
150+
// extract physical quantities
151+
const double density = values.at(HydroSystem<DustDamping>::density_index)[0];
152+
const double momentum_x = values.at(HydroSystem<DustDamping>::x1Momentum_index)[0];
153+
const double Egas_total = values.at(HydroSystem<DustDamping>::energy_index)[0];
154+
155+
// store gas velocity
156+
const double v_gas = momentum_x / density;
157+
userData_.v_gas_vec_.push_back(v_gas);
158+
159+
// store gas total energy
160+
userData_.E_gas_vec_.push_back(Egas_total);
161+
162+
if constexpr (Physics_Traits<DustDamping>::is_dust_enabled) {
163+
// store dust1 velocity
164+
const double dust1_density = values.at(HydroSystem<DustDamping>::dustDensity_index)[0];
165+
const double dust1_momentum_x = values.at(HydroSystem<DustDamping>::x1DustMomentum_index)[0];
166+
const double v_dust1 = dust1_momentum_x / dust1_density;
167+
userData_.v_dust1_vec_.push_back(v_dust1);
168+
169+
// store dust2 velocity
170+
const double dust2_density = values.at(HydroSystem<DustDamping>::dustDensity_index + numDustVars)[0];
171+
const double dust2_momentum_x = values.at(HydroSystem<DustDamping>::x1DustMomentum_index + numDustVars)[0];
172+
const double v_dust2 = dust2_momentum_x / dust2_density;
173+
userData_.v_dust2_vec_.push_back(v_dust2);
174+
}
175+
}
176+
}
177+
178+
auto run_simulation(double dt) -> SimulationData<DustDamping>
179+
{
180+
QuokkaSimulation<DustDamping> sim;
181+
182+
sim.reconstructionOrder_ = 3;
183+
sim.radiationReconstructionOrder_ = 3; // PPM
184+
sim.plotfileInterval_ = -1;
185+
sim.cflNumber_ = 1000000.0; // set large CFL to avoid CFL violation
186+
sim.constantDt_ = dt;
187+
188+
sim.setInitialConditions();
189+
190+
sim.evolve();
191+
192+
return sim.userData_;
193+
}
194+
195+
auto problem_main() -> int
196+
{
197+
amrex::Print() << "Running dust damping test with supersonic correction...\n";
198+
199+
const double dt_ref = 0.00001;
200+
const double dt_test = 0.005;
201+
202+
// step 1: run reference solution (dt = 0.00001)
203+
auto ref_data = run_simulation(dt_ref);
204+
205+
// step 2: run numerical solution (dt = 0.005)
206+
auto test_data = run_simulation(dt_test);
207+
208+
// step 3: calculate errors
209+
210+
// calculate sampling step size
211+
const auto step = static_cast<size_t>(std::round(dt_test / dt_ref));
212+
auto compute_relative_error = [](const std::vector<double> &sim,
213+
const std::vector<double> &ref, size_t step) {
214+
double err_sum = 0.0;
215+
double ref_sum = 0.0;
216+
int count = 0;
217+
218+
for (size_t i = 0; i < sim.size(); ++i) {
219+
size_t const ref_idx = i * step;
220+
if (ref_idx >= ref.size()) {
221+
break;
222+
}
223+
224+
err_sum += std::abs(sim[i] - ref[ref_idx]);
225+
ref_sum += std::abs(ref[ref_idx]);
226+
count++;
227+
}
228+
229+
if (count == 0 || ref_sum == 0.0) {
230+
return 1.0; // error value
231+
}
232+
return err_sum / ref_sum;
233+
};
234+
235+
double const rel_err_gas_vx = compute_relative_error(test_data.v_gas_vec_, ref_data.v_gas_vec_, step);
236+
double const rel_err_dust1_vx = compute_relative_error(test_data.v_dust1_vec_, ref_data.v_dust1_vec_, step);
237+
double const rel_err_dust2_vx = compute_relative_error(test_data.v_dust2_vec_, ref_data.v_dust2_vec_, step);
238+
double const rel_err_gas_E = compute_relative_error(test_data.E_gas_vec_, ref_data.E_gas_vec_, step);
239+
240+
amrex::Print() << "Comparison with reference solution (supersonic correction enabled):\n";
241+
amrex::Print() << "Relative L1 norm for gas vx = " << rel_err_gas_vx << "\n";
242+
amrex::Print() << "Relative L1 norm for dust1 vx = " << rel_err_dust1_vx << "\n";
243+
amrex::Print() << "Relative L1 norm for dust2 vx = " << rel_err_dust2_vx << "\n";
244+
amrex::Print() << "Relative L1 norm for gas E = " << rel_err_gas_E << "\n";
245+
246+
// determine whether the test has passed
247+
int status = 0;
248+
const double rel_err_tol = 0.01;
249+
250+
if ((rel_err_gas_vx > rel_err_tol) || (rel_err_dust1_vx > rel_err_tol) ||
251+
(rel_err_dust2_vx > rel_err_tol) || (rel_err_gas_E > rel_err_tol)) {
252+
status = 1;
253+
amrex::Print() << "Test FAILED: one or more errors exceed tolerance of "
254+
<< rel_err_tol << "\n";
255+
} else {
256+
amrex::Print() << "Test PASSED: all errors within tolerance of "
257+
<< rel_err_tol << "\n";
258+
}
259+
260+
#ifdef HAVE_PYTHON
261+
if (!ref_data.t_vec_.empty() && !test_data.t_vec_.empty()) {
262+
// gas velocity
263+
matplotlibcpp::clf();
264+
matplotlibcpp::plot(test_data.t_vec_, test_data.v_gas_vec_,
265+
{{"label", "numerical (iter, dt=0.005)"}, {"color", "r"},
266+
{"linestyle", "-"}, {"marker", "o"}, {"markersize", "3"}});
267+
matplotlibcpp::plot(ref_data.t_vec_, ref_data.v_gas_vec_,
268+
{{"label", "reference (non-iter, dt=0.00001)"},
269+
{"color", "r"}, {"linestyle", "--"}});
270+
matplotlibcpp::legend();
271+
matplotlibcpp::xlabel("t");
272+
matplotlibcpp::ylabel(R"($v_g$)");
273+
matplotlibcpp::title("Gas Velocity (with supersonic correction)");
274+
matplotlibcpp::tight_layout();
275+
matplotlibcpp::save("./dust_damping_iteration_gas_velocity.pdf");
276+
277+
// dust1 velocity
278+
matplotlibcpp::clf();
279+
matplotlibcpp::plot(test_data.t_vec_, test_data.v_dust1_vec_,
280+
{{"label", "numerical (iter, dt=0.005)"}, {"color", "b"},
281+
{"linestyle", "-"}, {"marker", "o"}, {"markersize", "3"}});
282+
matplotlibcpp::plot(ref_data.t_vec_, ref_data.v_dust1_vec_,
283+
{{"label", "reference (non-iter, dt=0.00001)"},
284+
{"color", "b"}, {"linestyle", "--"}});
285+
matplotlibcpp::legend();
286+
matplotlibcpp::xlabel("t");
287+
matplotlibcpp::ylabel(R"($v_{d,1}$)");
288+
matplotlibcpp::title("Dust1 Velocity (with supersonic correction)");
289+
matplotlibcpp::tight_layout();
290+
matplotlibcpp::save("./dust_damping_iteration_dust1_velocity.pdf");
291+
292+
// dust2 velocity
293+
matplotlibcpp::clf();
294+
matplotlibcpp::plot(test_data.t_vec_, test_data.v_dust2_vec_,
295+
{{"label", "numerical (iter, dt=0.005)"}, {"color", "g"},
296+
{"linestyle", "-"}, {"marker", "o"}, {"markersize", "3"}});
297+
matplotlibcpp::plot(ref_data.t_vec_, ref_data.v_dust2_vec_,
298+
{{"label", "reference (non-iter, dt=0.00001)"},
299+
{"color", "g"}, {"linestyle", "--"}});
300+
matplotlibcpp::legend();
301+
matplotlibcpp::xlabel("t");
302+
matplotlibcpp::ylabel(R"($v_{d,2}$)");
303+
matplotlibcpp::title("Dust2 Velocity (with supersonic correction)");
304+
matplotlibcpp::tight_layout();
305+
matplotlibcpp::save("./dust_damping_iteration_dust2_velocity.pdf");
306+
307+
// gas energy
308+
matplotlibcpp::clf();
309+
matplotlibcpp::plot(test_data.t_vec_, test_data.E_gas_vec_,
310+
{{"label", "numerical (iter, dt=0.005)"}, {"color", "m"},
311+
{"linestyle", "-"}, {"marker", "o"}, {"markersize", "3"}});
312+
matplotlibcpp::plot(ref_data.t_vec_, ref_data.E_gas_vec_,
313+
{{"label", "reference (non-iter, dt=0.00001)"},
314+
{"color", "m"}, {"linestyle", "--"}});
315+
matplotlibcpp::legend();
316+
matplotlibcpp::xlabel("t");
317+
matplotlibcpp::ylabel(R"($E_g$)");
318+
matplotlibcpp::title("Gas Energy (with supersonic correction)");
319+
matplotlibcpp::tight_layout();
320+
matplotlibcpp::save("./dust_damping_iteration_gas_energy.pdf");
321+
}
322+
#endif
323+
324+
amrex::Print() << "\nTest complete.\n";
325+
return status;
326+
}

0 commit comments

Comments
 (0)