-
Notifications
You must be signed in to change notification settings - Fork 12
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
improve convergence rate of the radiation Newton-Raphson solver (#720)
### Description
Changes:
1. Added a trick to make the Newton-Raphson iteration converge faster in
some situations: when $\Delta T > \max(T_{\rm gas}, T_{\rm rad})$, set
$T_{\rm gas}$ to $T_{\rm rad}$ and $R$ to 0. The motivation for doing
this trick is explained below.
2. Added `radiation.print_iteration_counts` parameter to enable counting
the mean and max number of Newton-Raphson iterations in each radiation
step. When set to `true`, the log file would look like:
```log
Coarse STEP 1 at t = 0 (0%) starts ...
time_subcycle = 0, total number of Newton-Raphson solvings is 128, (mean, max) number of Newton-Raphson iterations are 1, 1
Coarse STEP 2 at t = 0.0046875 (0.9375%) starts ...
time_subcycle = 0.0046875, total number of Newton-Raphson solvings is 128, (mean, max) number of Newton-Raphson iterations are 1, 1
```
3. Defined a new radiation problem, `RadMarshakDust`. Initially, the gas
is set with a density and temperature of 1, and the initial radiation
energy density is 0. Radiation with a temperature $T_r = 0.01$ streams
in from the left boundary. Before the dust model was introduced, we
assumed gas and dust were perfectly coupled, and the gas would be heated
by radiation where it could penetrate, and in regions without radiation,
the gas temperature would cool to 0.003159 (solution of $1 = a_R T^4 +
T$). In this test, however, we turned on the dust model and set the
dust-gas interaction coefficient to a very small number. As a result,
the dust and gas are decoupled: the dust temperature will reach the
radiation temperature and stay in equilibrium, while the gas remains at
its initial temperature.
In a future PR, I will extend this test by incorporating a multigroup
model. This will demonstrate how FUV radiation from the left boundary is
absorbed by the dust, which then re-radiates in IR.
### Related issues
None.
### 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)](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.
- [x] 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
1 parent
198cd1a
commit 76eff2f
Showing
7 changed files
with
386 additions
and
32 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,7 @@ | ||
| add_executable(test_radiation_marshak_dust test_radiation_marshak_dust.cpp ../../util/fextract.cpp ${QuokkaObjSources}) | ||
|
|
||
| if(AMReX_GPU_BACKEND MATCHES "CUDA") | ||
| setup_target_for_cuda_compilation(test_radiation_marshak_dust) | ||
| endif(AMReX_GPU_BACKEND MATCHES "CUDA") | ||
|
|
||
| add_test(NAME RadiationMarshakDust COMMAND test_radiation_marshak_dust RadMarshakDust.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) |
243 changes: 243 additions & 0 deletions
243
src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,243 @@ | ||
| //============================================================================== | ||
| // TwoMomentRad - a radiation transport library for patch-based AMR codes | ||
| // Copyright 2020 Benjamin Wibking. | ||
| // Released under the MIT license. See LICENSE file included in the GitHub repo. | ||
| //============================================================================== | ||
| /// \file test_radiation_marshak_dust.cpp | ||
| /// \brief Defines a test Marshak wave problem with weak coupling between dust and gas. | ||
| /// | ||
|
|
||
| #include "test_radiation_marshak_dust.hpp" | ||
| #include "AMReX.H" | ||
| #include "QuokkaSimulation.hpp" | ||
| #include "util/fextract.hpp" | ||
| #include "util/valarray.hpp" | ||
|
|
||
| struct StreamingProblem { | ||
| }; | ||
|
|
||
| constexpr double c = 1.0; // speed of light | ||
| constexpr double chat = 1.0; // reduced speed of light | ||
| constexpr double kappa0 = 10.0; // opacity | ||
| constexpr double rho0 = 1.0; | ||
| constexpr double CV = 1.0; | ||
| constexpr double mu = 1.5 / CV; // mean molecular weight | ||
| constexpr double initial_T = 1.0; | ||
| constexpr double a_rad = 1.0e10; | ||
| constexpr double erad_floor = 1.0e-10; | ||
| constexpr double initial_Erad = erad_floor; | ||
| constexpr double T_rad_L = 1.0e-2; // so EradL = 1e2 | ||
| constexpr double EradL = a_rad * T_rad_L * T_rad_L * T_rad_L * T_rad_L; | ||
| // constexpr double T_end_exact = 0.0031597766719577; // dust off; solution of 1 == a_rad * T^4 + T | ||
| constexpr double T_end_exact = initial_T; // dust on | ||
|
|
||
| template <> struct quokka::EOS_Traits<StreamingProblem> { | ||
| static constexpr double mean_molecular_weight = mu; | ||
| static constexpr double boltzmann_constant = 1.0; | ||
| static constexpr double gamma = 5. / 3.; | ||
| }; | ||
|
|
||
| template <> struct Physics_Traits<StreamingProblem> { | ||
| // cell-centred | ||
| static constexpr bool is_hydro_enabled = false; | ||
| static constexpr int numMassScalars = 0; // number of mass scalars | ||
| static constexpr int numPassiveScalars = numMassScalars + 0; // number of passive scalars | ||
| static constexpr bool is_radiation_enabled = true; | ||
| // face-centred | ||
| static constexpr bool is_mhd_enabled = false; | ||
| static constexpr int nGroups = 1; // number of radiation groups | ||
| }; | ||
|
|
||
| template <> struct RadSystem_Traits<StreamingProblem> { | ||
| static constexpr double c_light = c; | ||
| static constexpr double c_hat = chat; | ||
| static constexpr double radiation_constant = a_rad; | ||
| static constexpr double Erad_floor = erad_floor; | ||
| static constexpr int beta_order = 0; | ||
| static constexpr bool enable_dust_gas_thermal_coupling_model = true; | ||
| }; | ||
|
|
||
| template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real | ||
| { | ||
| return kappa0; | ||
| } | ||
|
|
||
| template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real | ||
| { | ||
| return kappa0; | ||
| } | ||
|
|
||
| template <> void QuokkaSimulation<StreamingProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem) | ||
| { | ||
| const amrex::Box &indexRange = grid_elem.indexRange_; | ||
| const amrex::Array4<double> &state_cc = grid_elem.array_; | ||
|
|
||
| const auto Erad0 = initial_Erad; | ||
| const auto Egas0 = initial_T * CV; | ||
|
|
||
| // loop over the grid and set the initial condition | ||
| amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { | ||
| for (int g = 0; g < Physics_Traits<StreamingProblem>::nGroups; ++g) { | ||
| state_cc(i, j, k, RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad0; | ||
| state_cc(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0; | ||
| state_cc(i, j, k, RadSystem<StreamingProblem>::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; | ||
| state_cc(i, j, k, RadSystem<StreamingProblem>::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; | ||
| } | ||
| state_cc(i, j, k, RadSystem<StreamingProblem>::gasEnergy_index) = Egas0; | ||
| state_cc(i, j, k, RadSystem<StreamingProblem>::gasDensity_index) = rho0; | ||
| state_cc(i, j, k, RadSystem<StreamingProblem>::gasInternalEnergy_index) = Egas0; | ||
| state_cc(i, j, k, RadSystem<StreamingProblem>::x1GasMomentum_index) = 0.; | ||
| state_cc(i, j, k, RadSystem<StreamingProblem>::x2GasMomentum_index) = 0.; | ||
| state_cc(i, j, k, RadSystem<StreamingProblem>::x3GasMomentum_index) = 0.; | ||
| }); | ||
| } | ||
|
|
||
| template <> | ||
| AMREX_GPU_DEVICE AMREX_FORCE_INLINE void | ||
| AMRSimulation<StreamingProblem>::setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4<amrex::Real> const &consVar, int /*dcomp*/, | ||
| int /*numcomp*/, amrex::GeometryData const &geom, const amrex::Real /*time*/, | ||
| const amrex::BCRec * /*bcr*/, int /*bcomp*/, int /*orig_comp*/) | ||
| { | ||
| #if (AMREX_SPACEDIM == 1) | ||
| auto i = iv.toArray()[0]; | ||
| int j = 0; | ||
| int k = 0; | ||
| #endif | ||
| #if (AMREX_SPACEDIM == 2) | ||
| auto [i, j] = iv.toArray(); | ||
| int k = 0; | ||
| #endif | ||
| #if (AMREX_SPACEDIM == 3) | ||
| auto [i, j, k] = iv.toArray(); | ||
| #endif | ||
|
|
||
| amrex::Box const &box = geom.Domain(); | ||
| amrex::GpuArray<int, 3> lo = box.loVect3d(); | ||
|
|
||
| if (i < lo[0]) { | ||
| // streaming inflow boundary | ||
| const double Erad = EradL; | ||
| const double Frad = c * Erad; | ||
|
|
||
| // multigroup radiation | ||
| // x1 left side boundary (Marshak) | ||
| for (int g = 0; g < Physics_Traits<StreamingProblem>::nGroups; ++g) { | ||
| consVar(i, j, k, RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erad; | ||
| consVar(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frad; | ||
| consVar(i, j, k, RadSystem<StreamingProblem>::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; | ||
| consVar(i, j, k, RadSystem<StreamingProblem>::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; | ||
| } | ||
| } | ||
|
|
||
| // gas boundary conditions are the same everywhere | ||
| const double Egas = initial_T * CV; | ||
| consVar(i, j, k, RadSystem<StreamingProblem>::gasEnergy_index) = Egas; | ||
| consVar(i, j, k, RadSystem<StreamingProblem>::gasDensity_index) = rho0; | ||
| consVar(i, j, k, RadSystem<StreamingProblem>::gasInternalEnergy_index) = Egas; | ||
| consVar(i, j, k, RadSystem<StreamingProblem>::x1GasMomentum_index) = 0.; | ||
| consVar(i, j, k, RadSystem<StreamingProblem>::x2GasMomentum_index) = 0.; | ||
| consVar(i, j, k, RadSystem<StreamingProblem>::x3GasMomentum_index) = 0.; | ||
| } | ||
|
|
||
| auto problem_main() -> int | ||
| { | ||
| // Problem parameters | ||
| // const int nx = 1000; | ||
| // const double Lx = 1.0; | ||
| const double CFL_number = 0.8; | ||
| const double dt_max = 1; | ||
| const int max_timesteps = 5000; | ||
|
|
||
| // Boundary conditions | ||
| constexpr int nvars = RadSystem<StreamingProblem>::nvar_; | ||
| amrex::Vector<amrex::BCRec> BCs_cc(nvars); | ||
| for (int n = 0; n < nvars; ++n) { | ||
| BCs_cc[n].setLo(0, amrex::BCType::ext_dir); // Dirichlet x1 | ||
| BCs_cc[n].setHi(0, amrex::BCType::foextrap); // extrapolate x1 | ||
| for (int i = 1; i < AMREX_SPACEDIM; ++i) { | ||
| BCs_cc[n].setLo(i, amrex::BCType::int_dir); // periodic | ||
| BCs_cc[n].setHi(i, amrex::BCType::int_dir); | ||
| } | ||
| } | ||
|
|
||
| // Problem initialization | ||
| QuokkaSimulation<StreamingProblem> sim(BCs_cc); | ||
|
|
||
| sim.radiationReconstructionOrder_ = 3; // PPM | ||
| // sim.stopTime_ = tmax; // set with runtime parameters | ||
| sim.radiationCflNumber_ = CFL_number; | ||
| sim.maxDt_ = dt_max; | ||
| sim.maxTimesteps_ = max_timesteps; | ||
| sim.plotfileInterval_ = -1; | ||
|
|
||
| // initialize | ||
| sim.setInitialConditions(); | ||
|
|
||
| // 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>(position.size()); | ||
|
|
||
| // compute error norm | ||
| std::vector<double> erad(nx); | ||
| std::vector<double> T(nx); | ||
| std::vector<double> T_exact(nx); | ||
| std::vector<double> xs(nx); | ||
| for (int i = 0; i < nx; ++i) { | ||
| amrex::Real const x = position[i]; | ||
| xs.at(i) = x; | ||
| double erad_sim = 0.0; | ||
| for (int g = 0; g < Physics_Traits<StreamingProblem>::nGroups; ++g) { | ||
| erad_sim += values.at(RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * g)[i]; | ||
| } | ||
| erad.at(i) = erad_sim; | ||
| const double e_gas = values.at(RadSystem<StreamingProblem>::gasInternalEnergy_index)[i]; | ||
| T.at(i) = quokka::EOS<StreamingProblem>::ComputeTgasFromEint(rho0, e_gas); | ||
| T_exact.at(i) = T_end_exact; | ||
| } | ||
|
|
||
| double err_norm = 0.; | ||
| double sol_norm = 0.; | ||
| for (int i = 0; i < nx; ++i) { | ||
| err_norm += std::abs(T[i] - T_exact[i]); | ||
| sol_norm += std::abs(T_exact[i]); | ||
| } | ||
|
|
||
| const double rel_err_norm = err_norm / sol_norm; | ||
| const double rel_err_tol = 0.02; | ||
| int status = 1; | ||
| if (rel_err_norm < rel_err_tol) { | ||
| status = 0; | ||
| } | ||
| amrex::Print() << "Relative L1 norm = " << rel_err_norm << std::endl; | ||
|
|
||
| #ifdef HAVE_PYTHON | ||
| // Plot results | ||
| matplotlibcpp::clf(); | ||
| std::map<std::string, std::string> plot_args; | ||
| plot_args["label"] = "numerical solution"; | ||
| matplotlibcpp::plot(xs, erad, plot_args); | ||
| matplotlibcpp::xlabel("x"); | ||
| matplotlibcpp::ylabel("E_rad"); | ||
| matplotlibcpp::legend(); | ||
| matplotlibcpp::title(fmt::format("t = {:f}", sim.tNew_[0])); | ||
| matplotlibcpp::tight_layout(); | ||
| matplotlibcpp::save("./radiation_marshak_dust_Erad.pdf"); | ||
|
|
||
| // plot temperature | ||
| matplotlibcpp::clf(); | ||
| matplotlibcpp::ylim(0.0, 1.1); | ||
| matplotlibcpp::plot(xs, T, plot_args); | ||
| matplotlibcpp::xlabel("x"); | ||
| matplotlibcpp::ylabel("Temperature"); | ||
| matplotlibcpp::title(fmt::format("t = {:f}", sim.tNew_[0])); | ||
| matplotlibcpp::tight_layout(); | ||
| matplotlibcpp::save("./radiation_marshak_dust_temperature.pdf"); | ||
| #endif // HAVE_PYTHON | ||
|
|
||
| // Cleanup and exit | ||
| amrex::Print() << "Finished." << std::endl; | ||
| return status; | ||
| } |
24 changes: 24 additions & 0 deletions
24
src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,24 @@ | ||
| #ifndef TEST_RADIATION_MARSHAK_DUST_HPP_ // NOLINT | ||
| #define TEST_RADIATION_MARSHAK_DUST_HPP_ | ||
| //============================================================================== | ||
| // TwoMomentRad - a radiation transport library for patch-based AMR codes | ||
| // Copyright 2020 Benjamin Wibking. | ||
| // Released under the MIT license. See LICENSE file included in the GitHub repo. | ||
| //============================================================================== | ||
| /// \file test_radiation_marshak_dust.hpp | ||
| /// \brief Defines a test Marshak wave problem with weak coupling between dust and gas. | ||
| /// | ||
|
|
||
| // external headers | ||
| #ifdef HAVE_PYTHON | ||
| #include "util/matplotlibcpp.h" | ||
| #endif | ||
| #include <fmt/format.h> | ||
|
|
||
| // internal headers | ||
|
|
||
| #include "radiation/radiation_system.hpp" | ||
|
|
||
| // function definitions | ||
|
|
||
| #endif // TEST_RADIATION_MARSHAK_DUST_HPP_ |
Oops, something went wrong.