Skip to content

Commit

Permalink
change initial condition of RadForce (#627)
Browse files Browse the repository at this point in the history 
### Description

The RadForce test can be improved to be more realistic and less `ad
hoc`. Instead of initializing with the exact solution and letting the
simulation maintain a static state, we start the simulation with a
uniform density of rho0 and zero velocity. The left boundary condition
is kept the same as before, but the right boundary is set to be
foextrap. The end result is exactly the same as before: the radiation
blows away the gas and the system comes to the static state given by the
differential equation.

### Checklist
_Before this pull request can be reviewed, all of these tasks should be
completed. Denote completed tasks with an `x` inside the square brackets
`[ ]` in the Markdown source below:_
- [x] I have added a description (see above).
- [ ] I have added a link to any related issues see (see above).
- [x] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [x] I have added tests for any new physics that this PR adds to the
code.
- [x] I have tested this PR on my local computer and all tests pass.
- [ ] I have manually triggered the GPU tests with the magic comment
`/azp run`.
- [x] I have requested a reviewer for this PR.

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
  • Loading branch information
@chongchonghe @pre-commit-ci
chongchonghe and pre-commit-ci[bot] authored May 12, 2024
1 parent 4ab068e commit d7f4b43
Show file tree
Hide file tree
Showing 2 changed files with 82 additions and 110 deletions.
184 changes: 75 additions & 109 deletions src/RadForce/test_radiation_force.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
/// \brief Defines a test problem for radiation force terms.
///

#include <cstdint>
#include <string>

#include "AMReX.H"
Expand Down Expand Up @@ -38,7 +39,6 @@ constexpr double tau = 1.0e-6; // optical depth (dimensionless)
constexpr double rho0 = 1.0e5 * mu; // g cm^-3
constexpr double Mach0 = 1.1; // Mach number at wind base
constexpr double Mach1 = 2.128410288469465339;
constexpr double rho1 = (Mach0 / Mach1) * rho0;

constexpr double Frad0 = rho0 * a0 * c_light_cgs_ / tau; // erg cm^-2 s^-1
constexpr double g0 = kappa0 * Frad0 / c_light_cgs_; // cm s^{-2}
Expand All @@ -60,7 +60,7 @@ template <> struct Physics_Traits<TubeProblem> {
// face-centred
static constexpr bool is_mhd_enabled = false;
// number of radiation groups
static constexpr int nGroups = 2;
static constexpr int nGroups = 1;
};

template <> struct RadSystem_Traits<TubeProblem> {
Expand All @@ -69,7 +69,6 @@ template <> struct RadSystem_Traits<TubeProblem> {
static constexpr double radiation_constant = radiation_constant_cgs_;
static constexpr double Erad_floor = 0.;
static constexpr double energy_unit = C::ev2erg;
static constexpr amrex::GpuArray<double, Physics_Traits<TubeProblem>::nGroups + 1> radBoundaries{0., 13.6, inf}; // eV
static constexpr int beta_order = 1;
};

Expand All @@ -89,68 +88,18 @@ AMREX_GPU_HOST_DEVICE auto RadSystem<TubeProblem>::ComputeFluxMeanOpacity(const
-> quokka::valarray<double, Physics_Traits<TubeProblem>::nGroups>
{
quokka::valarray<double, Physics_Traits<TubeProblem>::nGroups> kappaFVec{};
kappaFVec[0] = kappa0 * 1.5;
kappaFVec[1] = kappa0 * 0.5;
return kappaFVec;
}

// declare global variables
// initial conditions read from file
amrex::Gpu::HostVector<double> x_arr;
amrex::Gpu::HostVector<double> rho_arr;
amrex::Gpu::HostVector<double> Mach_arr;

amrex::Gpu::DeviceVector<double> x_arr_g;
amrex::Gpu::DeviceVector<double> rho_arr_g;
amrex::Gpu::DeviceVector<double> Mach_arr_g;

template <> void RadhydroSimulation<TubeProblem>::preCalculateInitialConditions()
{
std::string filename = "../extern/pressure_tube/optically_thin_wind.txt";
std::ifstream fstream(filename, std::ios::in);
AMREX_ALWAYS_ASSERT(fstream.is_open());
std::string header;
std::getline(fstream, header);

for (std::string line; std::getline(fstream, line);) {
std::istringstream iss(line);
std::vector<double> values;
for (double value = NAN; iss >> value;) {
values.push_back(value);
}
auto x = values.at(0); // position
auto rho = values.at(1); // density
auto Mach = values.at(2); // Mach number

x_arr.push_back(x);
rho_arr.push_back(rho);
Mach_arr.push_back(Mach);
for (int g = 0; g < nGroups_; ++g) {
kappaFVec[g] = kappa0;
}

// copy to device
x_arr_g.resize(x_arr.size());
rho_arr_g.resize(rho_arr.size());
Mach_arr_g.resize(Mach_arr.size());

amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, x_arr.begin(), x_arr.end(), x_arr_g.begin());
amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, rho_arr.begin(), rho_arr.end(), rho_arr_g.begin());
amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, Mach_arr.begin(), Mach_arr.end(), Mach_arr_g.begin());
amrex::Gpu::streamSynchronizeAll();
return kappaFVec;
}

template <> void RadhydroSimulation<TubeProblem>::setInitialConditionsOnGrid(quokka::grid grid_elem)
{
// extract variables required from the geom object
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx = grid_elem.dx_;
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> prob_lo = grid_elem.prob_lo_;
const amrex::Box &indexRange = grid_elem.indexRange_;
const amrex::Array4<double> &state_cc = grid_elem.array_;

auto const &x_ptr = x_arr_g.dataPtr();
auto const &rho_ptr = rho_arr_g.dataPtr();
auto const &Mach_ptr = Mach_arr_g.dataPtr();
int x_size = static_cast<int>(x_arr_g.size());

// calculate radEnergyFractions
quokka::valarray<amrex::Real, Physics_Traits<TubeProblem>::nGroups> radEnergyFractions{};
for (int g = 0; g < Physics_Traits<TubeProblem>::nGroups; ++g) {
Expand All @@ -159,12 +108,7 @@ template <> void RadhydroSimulation<TubeProblem>::setInitialConditionsOnGrid(quo

// loop over the grid and set the initial condition
amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
amrex::Real const x = (prob_lo[0] + (i + amrex::Real(0.5)) * dx[0]) / Lx;
amrex::Real const D = interpolate_value(x, x_ptr, rho_ptr, x_size);
amrex::Real const Mach = interpolate_value(x, x_ptr, Mach_ptr, x_size);

amrex::Real const rho = D * rho0;
amrex::Real const vel = Mach * a0;
amrex::Real const rho = rho0;

for (int g = 0; g < Physics_Traits<TubeProblem>::nGroups; ++g) {
state_cc(i, j, k, RadSystem<TubeProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) =
Expand All @@ -175,7 +119,7 @@ template <> void RadhydroSimulation<TubeProblem>::setInitialConditionsOnGrid(quo
}

state_cc(i, j, k, RadSystem<TubeProblem>::gasDensity_index) = rho;
state_cc(i, j, k, RadSystem<TubeProblem>::x1GasMomentum_index) = rho * vel;
state_cc(i, j, k, RadSystem<TubeProblem>::x1GasMomentum_index) = 0;
state_cc(i, j, k, RadSystem<TubeProblem>::x2GasMomentum_index) = 0;
state_cc(i, j, k, RadSystem<TubeProblem>::x3GasMomentum_index) = 0;
state_cc(i, j, k, RadSystem<TubeProblem>::gasEnergy_index) = 0;
Expand Down Expand Up @@ -204,7 +148,6 @@ AMRSimulation<TubeProblem>::setCustomBoundaryConditions(const amrex::IntVect &iv

amrex::Box const &box = geom.Domain();
amrex::GpuArray<int, 3> lo = box.loVect3d();
amrex::GpuArray<int, 3> hi = box.hiVect3d();

amrex::Real const Erad = Frad0 / c_light_cgs_;
amrex::Real const Frad = Frad0;
Expand All @@ -220,13 +163,6 @@ AMRSimulation<TubeProblem>::setCustomBoundaryConditions(const amrex::IntVect &iv
// left side
rho = rho0;
vel = Mach0 * a0;
} else if (i > hi[0]) {
// right side
rho = rho1;
vel = Mach1 * a0;
}

if ((i < lo[0]) || (i > hi[0])) {
// Dirichlet
for (int g = 0; g < Physics_Traits<TubeProblem>::nGroups; ++g) {
consVar(i, j, k, RadSystem<TubeProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad * radEnergyFractions[g];
Expand All @@ -249,6 +185,7 @@ auto problem_main() -> int
// Problem parameters
// const int nx = 128;
constexpr double CFL_number = 0.4;
double max_dt = 1.0e10;
constexpr double tmax = 10.0 * (Lx / a0);
constexpr int max_timesteps = 1e6;

Expand All @@ -258,7 +195,7 @@ auto problem_main() -> int
for (int n = 0; n < nvars; ++n) {
// for x-axis:
BCs_cc[n].setLo(0, amrex::BCType::ext_dir);
BCs_cc[n].setHi(0, amrex::BCType::ext_dir);
BCs_cc[n].setHi(0, amrex::BCType::foextrap);
// for y-, z- axes:
for (int i = 1; i < AMREX_SPACEDIM; ++i) {
// periodic
Expand All @@ -267,6 +204,10 @@ auto problem_main() -> int
}
}

// Read max_dt from parameter file
amrex::ParmParse const pp;
pp.query("max_dt", max_dt);

// Problem initialization
RadhydroSimulation<TubeProblem> sim(BCs_cc);

Expand All @@ -277,50 +218,74 @@ auto problem_main() -> int
sim.radiationCflNumber_ = CFL_number;
sim.maxTimesteps_ = max_timesteps;
sim.plotfileInterval_ = -1;
sim.maxDt_ = max_dt;

// initialize
sim.setInitialConditions();
auto [position0, values0] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0);

// evolve
sim.evolve();

// read output variables
auto [position, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0);
const int nx = static_cast<int>(position0.size());
const int nx = static_cast<int>(position.size());

// compute error norm
std::vector<double> xs(nx);
std::vector<double> xs_norm(nx);
std::vector<double> rho_arr(nx);
std::vector<double> rho_exact_arr(nx);
std::vector<double> rho_err(nx);
std::vector<double> vx_arr(nx);
std::vector<double> vx_exact_arr(nx);
std::vector<double> Frad_err(nx);
std::vector<double> Mach_arr(nx);

for (int i = 0; i < nx; ++i) {
xs.at(i) = position[i];
double rho_exact = values0.at(RadSystem<TubeProblem>::gasDensity_index)[i];
double x1GasMom_exact = values0.at(RadSystem<TubeProblem>::x1GasMomentum_index)[i];
double rho = values.at(RadSystem<TubeProblem>::gasDensity_index)[i];
double Frad = values.at(RadSystem<TubeProblem>::x1RadFlux_index)[i];
double x1GasMom = values.at(RadSystem<TubeProblem>::x1GasMomentum_index)[i];
double vx = x1GasMom / rho;
double vx_exact = x1GasMom_exact / rho_exact;

vx_arr.at(i) = vx / a0;
vx_exact_arr.at(i) = vx_exact / a0;
Frad_err.at(i) = (Frad - Frad0) / Frad0;
rho_err.at(i) = (rho - rho_exact) / rho_exact;
rho_exact_arr.at(i) = rho_exact;
xs_norm.at(i) = position[i] / Lx;

double const rho = values.at(RadSystem<TubeProblem>::gasDensity_index)[i];
double const x1GasMom = values.at(RadSystem<TubeProblem>::x1GasMomentum_index)[i];
double const vx = x1GasMom / rho;

rho_arr.at(i) = rho;
Mach_arr.at(i) = vx / a0;
}

// read in exact solution
std::vector<double> x_exact;
std::vector<double> x_exact_scaled;
std::vector<double> rho_exact;
std::vector<double> Mach_exact;

std::string const filename = "../extern/pressure_tube/optically_thin_wind.txt";
std::ifstream fstream(filename, std::ios::in);
AMREX_ALWAYS_ASSERT(fstream.is_open());
std::string header;
std::getline(fstream, header);

for (std::string line; std::getline(fstream, line);) {
std::istringstream iss(line);
std::vector<double> values;
for (double value = NAN; iss >> value;) {
values.push_back(value);
}
auto x = values.at(0); // position
auto rho = values.at(1); // density
auto Mach = values.at(2); // Mach number

x_exact.push_back(x);
x_exact_scaled.push_back(x * Lx);
rho_exact.push_back(rho * rho0);
Mach_exact.push_back(Mach);
}

// interpolate exact solution to simulation grid
std::vector<double> Mach_interp(nx);

interpolate_arrays(xs_norm.data(), Mach_interp.data(), nx, x_exact.data(), Mach_exact.data(), static_cast<int>(x_exact.size()));

double err_norm = 0.;
double sol_norm = 0.;
for (int i = 0; i < nx; ++i) {
err_norm += std::abs(rho_arr[i] - rho_exact_arr[i]);
sol_norm += std::abs(rho_exact_arr[i]);
err_norm += std::abs(Mach_arr[i] - Mach_interp[i]);
sol_norm += std::abs(Mach_interp[i]);
}

const double rel_err_norm = err_norm / sol_norm;
Expand All @@ -333,12 +298,14 @@ auto problem_main() -> int

#ifdef HAVE_PYTHON
// Plot density
std::map<std::string, std::string> rho_args;
std::unordered_map<std::string, std::string> rhoexact_args;
rho_args["label"] = "simulation";
std::map<std::string, std::string> rhoexact_args;
std::unordered_map<std::string, std::string> rho_args;
rhoexact_args["label"] = "exact solution";
matplotlibcpp::plot(xs, rho_arr, rho_args);
matplotlibcpp::scatter(xs, rho_exact_arr, 1.0, rhoexact_args);
rhoexact_args["color"] = "C0";
rho_args["label"] = "simulation";
rho_args["color"] = "C1";
matplotlibcpp::plot(x_exact_scaled, rho_exact, rhoexact_args);
matplotlibcpp::scatter(xs, rho_arr, 1.0, rho_args);
matplotlibcpp::legend();
matplotlibcpp::title(fmt::format("t = {:.4g} s", sim.tNew_[0]));
matplotlibcpp::xlabel("x (cm)");
Expand All @@ -347,18 +314,17 @@ auto problem_main() -> int
matplotlibcpp::save("./radiation_force_tube.pdf");

// plot velocity
int s = 4; // stride
std::map<std::string, std::string> vx_args;
std::unordered_map<std::string, std::string> vxexact_args;
vxexact_args["label"] = "exact solution";
int const s = nx / 64; // stride
std::map<std::string, std::string> vx_exact_args;
std::unordered_map<std::string, std::string> vx_args;
vx_exact_args["label"] = "exact solution";
vx_exact_args["color"] = "C0";
vx_args["label"] = "simulation";
vx_args["color"] = "C3";
vxexact_args["color"] = "C3";
vxexact_args["marker"] = "o";
// vxexact_args["edgecolors"] = "k";
vx_args["marker"] = "o";
vx_args["color"] = "C1";
matplotlibcpp::clf();
matplotlibcpp::plot(xs, vx_arr, vx_args);
matplotlibcpp::scatter(strided_vector_from(xs, s), strided_vector_from(vx_exact_arr, s), 10.0, vxexact_args);
matplotlibcpp::plot(x_exact_scaled, Mach_exact, vx_exact_args);
matplotlibcpp::scatter(strided_vector_from(xs, s), strided_vector_from(Mach_arr, s), 10.0, vx_args);
matplotlibcpp::legend();
// matplotlibcpp::title(fmt::format("t = {:.4g} s", sim.tNew_[0]));
matplotlibcpp::xlabel("length x (cm)");
Expand Down
8 changes: 7 additions & 1 deletion tests/RadForce.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,9 @@ geometry.prob_lo = 0.0 0.0 0.0
geometry.prob_hi = 1.0263747986171498e16 1.0263747986171498e16 1.0263747986171498e16
geometry.is_periodic = 0 1 1

# for convergence test, start with max_dt = 1.0e8
max_dt = 1.0e10

# *****************************************************************
# VERBOSITY
# *****************************************************************
Expand All @@ -20,4 +23,7 @@ amr.max_grid_size_x = 128

do_reflux = 0
do_subcycle = 0
suppress_output = 0
suppress_output = 1

cfl = 0.4
radiation.cfl = 0.4

0 comments on commit d7f4b43

Please sign in to comment.