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.

Sign up for GitHub

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

change initial condition of RadForce #627

Merged
merged 17 commits into from
May 12, 2024
Merged

change initial condition of RadForce #627

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
100f4fd
change initial condition for RadForce
chongchonghe May 4, 2024
9e7adb1
set rho to uniform
chongchonghe May 4, 2024
0201840
fixed tidy warning
chongchonghe May 4, 2024
c642652
fixed tidy warning
chongchonghe May 4, 2024
a715461
on macOS, at 2024-05-04 21:04:41
chongchonghe May 4, 2024
a0a3260
remove unused variable
chongchonghe May 5, 2024
67d434a
remove unused variable
chongchonghe May 5, 2024
dd06587
attempted to use setOutflowBoundaryLowOrder
chongchonghe May 7, 2024
0e529e6
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] May 7, 2024
d61c692
revert back to foextrap
chongchonghe May 7, 2024
a2885dc
fix GPU error
chongchonghe May 10, 2024
77d9df9
tidy
chongchonghe May 10, 2024
a1e0133
fixe division by 0 error
chongchonghe May 10, 2024
8e05ec6
fixe division by 0 error
chongchonghe May 10, 2024
cc0121a
Merge branch 'chongchong/convergence-test' of https://github.com/quok…
chongchonghe May 10, 2024
42f13d5
clean up force.cpp
chongchonghe May 12, 2024
32a96dc
address tidy comments and clean up
chongchonghe May 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 47 additions & 40 deletions src/RadForce/test_radiation_force.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,8 +88,9 @@ 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;
for (int g = 0; g < nGroups_; ++g) {
kappaFVec[g] = kappa0;
}
return kappaFVec;
}

Expand Down Expand Up @@ -141,16 +141,9 @@ template <> void RadhydroSimulation<TubeProblem>::preCalculateInitialConditions(
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 +152,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 +163,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 +192,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 +207,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 +229,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 +239,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); // extrapolate
// for y-, z- axes:
for (int i = 1; i < AMREX_SPACEDIM; ++i) {
// periodic
Expand All @@ -267,6 +248,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,6 +262,7 @@ auto problem_main() -> int
sim.radiationCflNumber_ = CFL_number;
sim.maxTimesteps_ = max_timesteps;
sim.plotfileInterval_ = -1;
sim.maxDt_ = max_dt;

// initialize
sim.setInitialConditions();
Expand All @@ -287,7 +273,15 @@ auto problem_main() -> int

// 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_tot = static_cast<int>(position0.size());

// cut off at x = 1
int nx = 0;
for (int i = 0; i < nx_tot; ++i) {
if (position[i] < 0.98e16) {
nx++;
}
}

// compute error norm
std::vector<double> xs(nx);
Expand All @@ -298,15 +292,26 @@ auto problem_main() -> int
std::vector<double> vx_exact_arr(nx);
std::vector<double> Frad_err(nx);

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 const x_size = static_cast<int>(x_arr_g.size());

for (int i = 0; i < nx; ++i) {
amrex::Real const x_ = position[i] / 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;

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;
double const rho_exact = _rho;
// double x1GasMom_exact = values0.at(RadSystem<TubeProblem>::x1GasMomentum_index)[i];
double const rho = values.at(RadSystem<TubeProblem>::gasDensity_index)[i];
double const Frad = values.at(RadSystem<TubeProblem>::x1RadFlux_index)[i];
double const x1GasMom = values.at(RadSystem<TubeProblem>::x1GasMomentum_index)[i];
double const vx = x1GasMom / rho;
double const vx_exact = _vel;

vx_arr.at(i) = vx / a0;
vx_exact_arr.at(i) = vx_exact / a0;
Expand All @@ -319,8 +324,8 @@ auto problem_main() -> int
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(vx_arr[i] - vx_exact_arr[i]);
sol_norm += std::abs(vx_exact_arr[i]);
}

const double rel_err_norm = err_norm / sol_norm;
Expand All @@ -336,7 +341,9 @@ auto problem_main() -> int
std::map<std::string, std::string> rho_args;
std::unordered_map<std::string, std::string> rhoexact_args;
rho_args["label"] = "simulation";
rho_args["color"] = "C0";
rhoexact_args["label"] = "exact solution";
rhoexact_args["color"] = "C1";
matplotlibcpp::plot(xs, rho_arr, rho_args);
matplotlibcpp::scatter(xs, rho_exact_arr, 1.0, rhoexact_args);
matplotlibcpp::legend();
Expand All @@ -347,14 +354,14 @@ auto problem_main() -> int
matplotlibcpp::save("./radiation_force_tube.pdf");

// plot velocity
int s = 4; // stride
int const s = nx_tot / 64; // stride
std::map<std::string, std::string> vx_args;
std::unordered_map<std::string, std::string> vxexact_args;
vxexact_args["label"] = "exact solution";
vx_args["label"] = "simulation";
vx_args["color"] = "C3";
vxexact_args["color"] = "C3";
vx_args["color"] = "C0";
vxexact_args["marker"] = "o";
vxexact_args["color"] = "C1";
// vxexact_args["edgecolors"] = "k";
matplotlibcpp::clf();
matplotlibcpp::plot(xs, vx_arr, vx_args);
Expand Down
8 changes: 7 additions & 1 deletion tests/RadForce.in
View file
Open in desktop
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