Skip to content

Commit

Permalink
photoelectric heating (#766)
Browse files Browse the repository at this point in the history 
### Description
In this PR, I implement photoelectric heating within the current
gas-radiation-dust interaction framework. PE can be enabled by setting
`enable_photoelectric_heating = true` in `ISM_Traits`. The prerequisite
is `nGroups > 1` and `enable_dust_gas_thermal_coupling_model = true`.

Added two test problems, `RadiationMarshakDustPE-coupled` and
`RadiationMarshakDustPE-decoupled`, to test photoelectric heating in
cases where gas and dust are well coupled and decoupled. These two tests
share the same problem code `test_radiation_marshak_dust_and_PE.cpp` but
use different runtime parameters.

Here is the code structure:

```
- AddSourceTermsMultiGroup
  - outer loop:
    - compute `v_times_F_term` 
    - solve gas-dust-radiation energy exchange:
      - if constexpr (!enable_dust_gas_thermal_coupling_model_)
        - `SolveGasRadiationEnergyExchange`:
          - The matter-radiation coupling scheme from the multigroup paper
      - else if constexpr (!enable_photoelectric_heating_)
        - `SolveGasDustRadiationEnergyExchange()`:
          - The gas-dust-radiation coupling scheme proposed in #733, which has two cases: the well-coupled case and decoupled case
      - else
        - `SolveGasDustRadiationEnergyExchangeWithPE()`:
          - Based on the gas-dust-radiation coupling scheme, but with modified Jacobian in the well-coupled case to include PE heating. In the decoupled case, PE heating is added via a forward-Euler update using the updated radiation energy density after the gas-dust-radiation step.
    - compute updated `v_times_F_term` 
    - if `v_times_F_term` has converged, break
  - `UpdateFlux()`
```

### Related issues
A modification to #733 to include PE heating. 

### 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).
- [x] 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.
- [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
@chongchonghe @pre-commit-ci
chongchonghe and pre-commit-ci[bot] authored Oct 4, 2024
1 parent 8671c5e commit e69f403
Show file tree
Hide file tree
Showing 15 changed files with 1,509 additions and 327 deletions.
26 changes: 17 additions & 9 deletions src/QuokkaSimulation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1617,10 +1617,10 @@ void QuokkaSimulation<problem_t>::subcycleRadiationAtLevel(int lev, amrex::Real
// new radiation state is stored in state_new_cc_
// new hydro state is stored in state_new_cc_ (always the case during radiation update)

amrex::Gpu::Buffer<int> iteration_failure_counter(
{0, 0, 0}); // failure counter for: matter-radiation coupling, dust temperature, outer iteration
amrex::Gpu::Buffer<int> iteration_counter(
{0, 0, 0}); // iteration counter for: radiation update, Newton-Raphson iterations, max Newton-Raphson iterations
// failure counter for: matter-radiation coupling, dust temperature, outer iteration
amrex::Gpu::Buffer<int> iteration_failure_counter({0, 0, 0});
// iteration counter for: radiation update, Newton-Raphson iterations, max Newton-Raphson iterations, decoupled gas-dust update
amrex::Gpu::Buffer<int> iteration_counter({0, 0, 0, 0});
int *p_iteration_failure_counter = iteration_failure_counter.data();
int *p_iteration_counter = iteration_counter.data();

Expand Down Expand Up @@ -1659,13 +1659,15 @@ void QuokkaSimulation<problem_t>::subcycleRadiationAtLevel(int lev, amrex::Real

if (print_rad_counter_) {
auto h_iteration_counter = iteration_counter.copyToHost();
long global_solver_count = h_iteration_counter[0]; // number of Newton-Raphson solvings
long global_iteration_sum = h_iteration_counter[1]; // sum of Newton-Raphson iterations
int global_iteration_max = h_iteration_counter[2]; // max number of Newton-Raphson iterations
long global_solver_count = h_iteration_counter[0]; // number of Newton-Raphson solvings
long global_iteration_sum = h_iteration_counter[1]; // sum of Newton-Raphson iterations
int global_iteration_max = h_iteration_counter[2]; // max number of Newton-Raphson iterations
long global_decoupled_iteration_sum = h_iteration_counter[3]; // sum of decoupled gas-dust Newton-Raphson iterations

amrex::ParallelDescriptor::ReduceLongSum(global_solver_count);
amrex::ParallelDescriptor::ReduceLongSum(global_iteration_sum);
amrex::ParallelDescriptor::ReduceIntMax(global_iteration_max);
amrex::ParallelDescriptor::ReduceLongSum(global_decoupled_iteration_sum);

if (amrex::ParallelDescriptor::IOProcessor()) {
const auto n_cells = CountCells(lev);
Expand All @@ -1675,9 +1677,15 @@ void QuokkaSimulation<problem_t>::subcycleRadiationAtLevel(int lev, amrex::Real
static_cast<double>(global_iteration_sum) / static_cast<double>(global_solver_count);
const double global_solving_mean =
static_cast<double>(global_solver_count) / static_cast<double>(n_cells) / 2.0; // 2 stages
amrex::Print() << "average number of Newton-Raphson solvings per IMEX stage is " << global_solving_mean
const double global_decoupled_iteration_mean =
static_cast<double>(global_decoupled_iteration_sum) / static_cast<double>(global_solver_count);
amrex::Print() << "The average number of Newton-Raphson solvings per IMEX stage is " << global_solving_mean
<< ", (mean, max) number of Newton-Raphson iterations are " << global_iteration_mean << ", "
<< global_iteration_max << "\n";
<< global_iteration_max << ".\n";
if constexpr (RadSystem_Traits<problem_t>::enable_dust_gas_thermal_coupling_model) {
amrex::Print() << "The fraction of gas-dust interactions that are decoupled is "
<< global_decoupled_iteration_mean << "\n";
}
}
}
}
Expand Down
1 change: 1 addition & 0 deletions src/problems/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ add_subdirectory(RadTube)
add_subdirectory(RadDust)
add_subdirectory(RadDustMG)
add_subdirectory(RadMarshakDust)
add_subdirectory(RadMarshakDustPE)

add_subdirectory(RadhydroShell)
add_subdirectory(RadhydroShock)
Expand Down
2 changes: 1 addition & 1 deletion src/problems/RadDustMG/test_rad_dust_MG.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
// internal headers

#include "math/interpolate.hpp"
#include "radiation/radiation_system.hpp"
#include "radiation/radiation_dust_system.hpp"

// function definitions

Expand Down
89 changes: 47 additions & 42 deletions src/problems/RadMarshakDust/test_radiation_marshak_dust.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
#include "util/fextract.hpp"
#include "util/valarray.hpp"

struct StreamingProblem {
struct MarshakProblem {
};

AMREX_GPU_MANAGED double kappa1 = NAN; // dust opacity at IR
Expand Down Expand Up @@ -42,13 +42,13 @@ constexpr int n_group_ = 2;
static constexpr amrex::GpuArray<double, n_group_ + 1> radBoundaries_{1e-10, 100, 1e4};
static constexpr OpacityModel opacity_model_ = OpacityModel::piecewise_constant_opacity;

template <> struct quokka::EOS_Traits<StreamingProblem> {
template <> struct quokka::EOS_Traits<MarshakProblem> {
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> {
template <> struct Physics_Traits<MarshakProblem> {
// cell-centred
static constexpr bool is_hydro_enabled = false;
static constexpr int numMassScalars = 0; // number of mass scalars
Expand All @@ -59,7 +59,7 @@ template <> struct Physics_Traits<StreamingProblem> {
static constexpr int nGroups = n_group_; // number of radiation groups
};

template <> struct RadSystem_Traits<StreamingProblem> {
template <> struct RadSystem_Traits<MarshakProblem> {
static constexpr double c_light = c;
static constexpr double c_hat = chat;
static constexpr double radiation_constant = a_rad;
Expand All @@ -71,20 +71,25 @@ template <> struct RadSystem_Traits<StreamingProblem> {
static constexpr OpacityModel opacity_model = opacity_model_;
};

template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
template <> struct ISM_Traits<MarshakProblem> {
static constexpr double gas_dust_coupling_threshold = 1.0e-5;
static constexpr bool enable_photoelectric_heating = false;
};

template <> AMREX_GPU_HOST_DEVICE auto RadSystem<MarshakProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
{
return kappa1;
}

template <> AMREX_GPU_HOST_DEVICE auto RadSystem<StreamingProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
template <> AMREX_GPU_HOST_DEVICE auto RadSystem<MarshakProblem>::ComputeFluxMeanOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
{
return kappa1;
}

template <>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto
RadSystem<StreamingProblem>::DefineOpacityExponentsAndLowerValues(amrex::GpuArray<double, nGroups_ + 1> /*rad_boundaries*/, const double /*rho*/,
const double /*Tgas*/) -> amrex::GpuArray<amrex::GpuArray<double, nGroups_ + 1>, 2>
RadSystem<MarshakProblem>::DefineOpacityExponentsAndLowerValues(amrex::GpuArray<double, nGroups_ + 1> /*rad_boundaries*/, const double /*rho*/,
const double /*Tgas*/) -> amrex::GpuArray<amrex::GpuArray<double, nGroups_ + 1>, 2>
{
amrex::GpuArray<amrex::GpuArray<double, nGroups_ + 1>, 2> exponents_and_values{};
for (int i = 0; i < nGroups_ + 1; ++i) {
Expand All @@ -98,36 +103,36 @@ RadSystem<StreamingProblem>::DefineOpacityExponentsAndLowerValues(amrex::GpuArra
return exponents_and_values;
}

template <> void QuokkaSimulation<StreamingProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)
template <> void QuokkaSimulation<MarshakProblem>::setInitialConditionsOnGrid(quokka::grid const &grid_elem)
{
const amrex::Box &indexRange = grid_elem.indexRange_;
const amrex::Array4<double> &state_cc = grid_elem.array_;

const auto Egas0 = initial_T * CV;
const auto Erads = RadSystem<StreamingProblem>::ComputeThermalRadiationMultiGroup(initial_Trad, radBoundaries_);
const auto Erads = RadSystem<MarshakProblem>::ComputeThermalRadiationMultiGroup(initial_Trad, radBoundaries_);

// 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) = Erads[g];
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;
for (int g = 0; g < Physics_Traits<MarshakProblem>::nGroups; ++g) {
state_cc(i, j, k, RadSystem<MarshakProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g];
state_cc(i, j, k, RadSystem<MarshakProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
state_cc(i, j, k, RadSystem<MarshakProblem>::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
state_cc(i, j, k, RadSystem<MarshakProblem>::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.;
state_cc(i, j, k, RadSystem<MarshakProblem>::gasEnergy_index) = Egas0;
state_cc(i, j, k, RadSystem<MarshakProblem>::gasDensity_index) = rho0;
state_cc(i, j, k, RadSystem<MarshakProblem>::gasInternalEnergy_index) = Egas0;
state_cc(i, j, k, RadSystem<MarshakProblem>::x1GasMomentum_index) = 0.;
state_cc(i, j, k, RadSystem<MarshakProblem>::x2GasMomentum_index) = 0.;
state_cc(i, j, k, RadSystem<MarshakProblem>::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*/)
AMRSimulation<MarshakProblem>::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];
Expand All @@ -145,7 +150,7 @@ AMRSimulation<StreamingProblem>::setCustomBoundaryConditions(const amrex::IntVec
amrex::Box const &box = geom.Domain();
amrex::GpuArray<int, 3> lo = box.loVect3d();

// const auto Erads = RadSystem<StreamingProblem>::ComputeThermalRadiation(T_rad_L, radBoundaries_);
// const auto Erads = RadSystem<MarshakProblem>::ComputeThermalRadiation(T_rad_L, radBoundaries_);
quokka::valarray<double, 2> const Erads = {erad_floor, EradL};
const double c_light = c;
const auto Frads = Erads * c_light;
Expand All @@ -154,22 +159,22 @@ AMRSimulation<StreamingProblem>::setCustomBoundaryConditions(const amrex::IntVec
// streaming inflow boundary
// 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) = Erads[g];
consVar(i, j, k, RadSystem<StreamingProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frads[g];
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;
for (int g = 0; g < Physics_Traits<MarshakProblem>::nGroups; ++g) {
consVar(i, j, k, RadSystem<MarshakProblem>::radEnergy_index + Physics_NumVars::numRadVars * g) = Erads[g];
consVar(i, j, k, RadSystem<MarshakProblem>::x1RadFlux_index + Physics_NumVars::numRadVars * g) = Frads[g];
consVar(i, j, k, RadSystem<MarshakProblem>::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0;
consVar(i, j, k, RadSystem<MarshakProblem>::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.;
consVar(i, j, k, RadSystem<MarshakProblem>::gasEnergy_index) = Egas;
consVar(i, j, k, RadSystem<MarshakProblem>::gasDensity_index) = rho0;
consVar(i, j, k, RadSystem<MarshakProblem>::gasInternalEnergy_index) = Egas;
consVar(i, j, k, RadSystem<MarshakProblem>::x1GasMomentum_index) = 0.;
consVar(i, j, k, RadSystem<MarshakProblem>::x2GasMomentum_index) = 0.;
consVar(i, j, k, RadSystem<MarshakProblem>::x3GasMomentum_index) = 0.;
}

auto problem_main() -> int
Expand All @@ -187,7 +192,7 @@ auto problem_main() -> int
pp.query("kappa2", kappa2);

// Boundary conditions
constexpr int nvars = RadSystem<StreamingProblem>::nvar_;
constexpr int nvars = RadSystem<MarshakProblem>::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
Expand All @@ -199,7 +204,7 @@ auto problem_main() -> int
}

// Problem initialization
QuokkaSimulation<StreamingProblem> sim(BCs_cc);
QuokkaSimulation<MarshakProblem> sim(BCs_cc);

sim.radiationReconstructionOrder_ = 3; // PPM
// sim.stopTime_ = tmax; // set with runtime parameters
Expand Down Expand Up @@ -230,14 +235,14 @@ auto problem_main() -> int
for (int i = 0; i < nx; ++i) {
amrex::Real const x = position[i];
xs.at(i) = x;
erad1.at(i) = values.at(RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * 0)[i];
erad1.at(i) = values.at(RadSystem<MarshakProblem>::radEnergy_index + Physics_NumVars::numRadVars * 0)[i];
erad.at(i) = erad1.at(i);
if (n_group_ > 1) {
erad2.at(i) = values.at(RadSystem<StreamingProblem>::radEnergy_index + Physics_NumVars::numRadVars * 1)[i];
erad2.at(i) = values.at(RadSystem<MarshakProblem>::radEnergy_index + Physics_NumVars::numRadVars * 1)[i];
erad.at(i) += erad2.at(i);
}
const double e_gas = values.at(RadSystem<StreamingProblem>::gasInternalEnergy_index)[i];
T.at(i) = quokka::EOS<StreamingProblem>::ComputeTgasFromEint(rho0, e_gas);
const double e_gas = values.at(RadSystem<MarshakProblem>::gasInternalEnergy_index)[i];
T.at(i) = quokka::EOS<MarshakProblem>::ComputeTgasFromEint(rho0, e_gas);
T_exact.at(i) = T_end_exact;

erad2_exact.at(i) = x < sim.tNew_[0] ? EradL * std::exp(-x * rho0 * kappa2) : erad_floor;
Expand Down
2 changes: 1 addition & 1 deletion src/problems/RadMarshakDust/test_radiation_marshak_dust.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@

// internal headers

#include "radiation/radiation_system.hpp"
#include "radiation/radiation_dust_system.hpp"

// function definitions

Expand Down
8 changes: 8 additions & 0 deletions src/problems/RadMarshakDustPE/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
add_executable(test_radiation_marshak_dust_PE test_radiation_marshak_dust_and_PE.cpp ../../util/fextract.cpp ${QuokkaObjSources})

if(AMReX_GPU_BACKEND MATCHES "CUDA")
setup_target_for_cuda_compilation(test_radiation_marshak_dust_PE)
endif(AMReX_GPU_BACKEND MATCHES "CUDA")

add_test(NAME RadiationMarshakDustPE-coupled COMMAND test_radiation_marshak_dust_PE RadMarshakDustPEcoupled.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
add_test(NAME RadiationMarshakDustPE-decoupled COMMAND test_radiation_marshak_dust_PE RadMarshakDustPEdecoupled.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
Loading

0 comments on commit e69f403

Please sign in to comment.