Skip to content

Commit

Permalink
PPL opacity models (#678)
Browse files Browse the repository at this point in the history
### Description

Main updates:

- Defined three opacity models described in Paper III: piecewise
constant (`piecewise_constant_opacity`), piecewise powerlaw with fixed
slope (`PPL_opacity_fixed_slope_spectrum`), piecewise powerlaw with full
spectrum fitting (`PPL_opacity_full_spectrum`).

- Added tests for each of these methods. 
- RadhydroBB tests the PC method's ability to resolve the
Doppler-shifted spectrum in a moving blackbody.
- RadMarshakVaytet tests the ability of all three methods (by default
`PPL_opacity_full_spectrum`) to deal with frequency-dependent velocities
within an energy bin.
- RadTube tests that in the special case of opacity being independent of
frequency, the `PPL_opacity_fixed_slope_spectrum` method provides the
same result as the single-group method.
- RadhydroShockMultigroup tests `PPL_opacity_full_spectrum` in the
special case of constant opacity.

More updates:

- Removed `ComputePlanckOpacityTempDerivative`
- Defined `ComputeDiffusionFluxMeanOpacity` and removed the calculation
of component-specific flux-mean opacity. See Paper III.
- Defined `ComputeBinCenterOpacity` to use when `opacity_model_ ==
OpacityModel::PPL_opacity_full_spectrum` and
`use_diffuse_flux_mean_opacity = false`.
- Changed the structure of the Newton-Raphson iteration loop. We check
the residuals first, then calculate opacities and the Jacobian.
- Fixed a bug in the Newton iteration where the calculation of alpha_E
uses `Erad0Vec` but `EradVec_guess` should be used.
- More robust break condition for the Newton iteration.
- Stop updating `alpha_E` and `alpha_B` after
`max_ite_to_update_alpha_E` steps in the Newton iteration.
- Set `reconstructionOrder_=3` in the RadTube and
RadhydroShockMultigroup tests to avoid "erroneous arithmetic" error
after using `-DCMAKE_CXX_FLAGS=-ffp-exception-behavior=strict`.



### 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)](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 and pre-commit-ci[bot] authored Jul 24, 2024
1 parent d5210ff commit 4686c6a
Show file tree
Hide file tree
Showing 21 changed files with 2,888 additions and 502 deletions.
1,025 changes: 1,025 additions & 0 deletions extern/Doppler-spectrum/exact_flux_density.csv

Large diffs are not rendered by default.

4 changes: 3 additions & 1 deletion src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -166,6 +166,7 @@ add_subdirectory(RadForce)
add_subdirectory(RadMarshak)
add_subdirectory(RadMarshakCGS)
add_subdirectory(RadMarshakAsymptotic)
add_subdirectory(RadMarshakVaytet)
add_subdirectory(RadMatterCoupling)
add_subdirectory(RadMatterCouplingRSLA)
add_subdirectory(RadPulse)
Expand All @@ -184,8 +185,9 @@ add_subdirectory(RadhydroUniformAdvecting)
add_subdirectory(RadhydroPulse)
add_subdirectory(RadhydroPulseDyn)
add_subdirectory(RadhydroPulseGrey)
add_subdirectory(RadhydroPulseMG)
add_subdirectory(RadhydroPulseMGconst)
add_subdirectory(RadhydroPulseMGint)
add_subdirectory(RadhydroBB)

add_subdirectory(BinaryOrbitCIC)
add_subdirectory(Cooling)
Expand Down
12 changes: 0 additions & 12 deletions src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -71,18 +71,6 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem<SuOlsonProblemCgs>::Comp
return ComputePlanckOpacity(rho, Tgas);
}

template <>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto
RadSystem<SuOlsonProblemCgs>::ComputePlanckOpacityTempDerivative(const double rho, const double Tgas) -> quokka::valarray<double, nGroups_>
{
quokka::valarray<double, nGroups_> opacity_deriv{};
auto sigma_dT = (-3.0 * kappa / Tgas) * std::pow(Tgas / T_hohlraum, -3); // cm^-1
for (int i = 0; i < nGroups_; ++i) {
opacity_deriv[i] = sigma_dT / rho;
}
return opacity_deriv;
}

template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem<SuOlsonProblemCgs>::ComputeEddingtonFactor(double /*f*/) -> double
{
return (1. / 3.); // Eddington approximation
Expand Down
7 changes: 7 additions & 0 deletions src/RadMarshakVaytet/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
add_executable(test_radiation_marshak_Vaytet test_radiation_marshak_Vaytet.cpp ../fextract.cpp ../interpolate.cpp ${QuokkaObjSources})

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

add_test(NAME MarshakWaveVaytet COMMAND test_radiation_marshak_Vaytet MarshakVaytet.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
486 changes: 486 additions & 0 deletions src/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp

Large diffs are not rendered by default.

26 changes: 26 additions & 0 deletions src/RadMarshakVaytet/test_radiation_marshak_Vaytet.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
#ifndef TEST_RADIATION_MARSHAK_VAYTET_HPP_ // NOLINT
#define TEST_RADIATION_MARSHAK_VAYTET_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_vaytet.hpp
/// \brief Defines a test problem for radiation in the static diffusion regime.
///

// external headers
#ifdef HAVE_PYTHON
#include "matplotlibcpp.h"
#endif
#include <fmt/format.h>
#include <fstream>

// internal headers

#include "interpolate.hpp"
#include "radiation_system.hpp"

// function definitions

#endif // TEST_RADIATION_MARSHAK_VAYTET_HPP_
13 changes: 0 additions & 13 deletions src/RadPulse/test_radiation_pulse.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -81,19 +81,6 @@ AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::ComputeFluxMeanOpacity(const
return ComputePlanckOpacity(rho, Tgas);
}

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::ComputePlanckOpacityTempDerivative(const double rho,
const double Tgas) -> quokka::valarray<double, nGroups_>
{
quokka::valarray<double, nGroups_> opacity_deriv{};
double opacity_deriv_scalar = 0.;
if (Tgas > 1.0) {
opacity_deriv_scalar = (kappa0 / rho) * (3.0 / T0) * std::pow(Tgas / T0, 2);
}
opacity_deriv.fillin(opacity_deriv_scalar);
return opacity_deriv;
}

template <> void RadhydroSimulation<PulseProblem>::setInitialConditionsOnGrid(quokka::grid grid_elem)
{
// extract variables required from the geom object
Expand Down
76 changes: 47 additions & 29 deletions src/RadTube/test_radiation_tube.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
// Released under the MIT license. See LICENSE file included in the GitHub repo.
//==============================================================================
/// \file test_radiation_tube.cpp
/// \brief Defines a test problem for radiation pressure terms.
/// \brief Defines a test problem for radiation pressure terms. This is also a trivial test for the PPL_fixed_slope opacity model.
///

#include <string>
Expand Down Expand Up @@ -36,6 +36,7 @@ constexpr double rho0 = 1.0; // g cm^-3
constexpr double T0 = 2.75e7; // K
constexpr double rho1 = 2.1940476649492044; // g cm^-3
constexpr double T1 = 2.2609633884436745e7; // K
constexpr double a_rad = C::a_rad;

constexpr double a0 = 4.0295519855200705e7; // cm s^-1

Expand All @@ -60,28 +61,45 @@ template <> struct Physics_Traits<TubeProblem> {
template <> struct RadSystem_Traits<TubeProblem> {
static constexpr double c_light = c_light_cgs_;
static constexpr double c_hat = 10.0 * a0;
static constexpr double radiation_constant = radiation_constant_cgs_;
static constexpr double radiation_constant = a_rad;
static constexpr double Erad_floor = 0.;
static constexpr double energy_unit = C::k_B;
static constexpr amrex::GpuArray<double, Physics_Traits<TubeProblem>::nGroups + 1> radBoundaries{0., 3.3 * T0, inf}; // Kelvin
static constexpr amrex::GpuArray<double, Physics_Traits<TubeProblem>::nGroups + 1> radBoundaries{0.01 * T0, 3.3 * T0, 1000. * T0}; // Kelvin
// static constexpr amrex::GpuArray<double, Physics_Traits<TubeProblem>::nGroups + 1> radBoundaries{0.01 * T0, 1000. * T0}; // Kelvin
static constexpr int beta_order = 1;
static constexpr OpacityModel opacity_model = OpacityModel::user;
// static constexpr OpacityModel opacity_model = OpacityModel::single_group;
static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity;
// static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum;
};

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<TubeProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto
RadSystem<TubeProblem>::DefineOpacityExponentsAndLowerValues(amrex::GpuArray<double, nGroups_ + 1> /*rad_boundaries*/, const double /*rho*/,
const double /*Tgas*/) -> amrex::GpuArray<amrex::GpuArray<double, nGroups_ + 1>, 2>
{
quokka::valarray<double, nGroups_> kappaPVec{};
for (int g = 0; g < nGroups_; ++g) {
kappaPVec[g] = kappa0;
amrex::GpuArray<amrex::GpuArray<double, nGroups_ + 1>, 2> exponents_and_values{};
for (int i = 0; i < nGroups_ + 1; ++i) {
exponents_and_values[0][i] = 0.0;
exponents_and_values[1][i] = kappa0;
}
return kappaPVec;
return exponents_and_values;
}

template <> AMREX_GPU_HOST_DEVICE auto RadSystem<TubeProblem>::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray<double, nGroups_>
{
return ComputePlanckOpacity(rho, Tgas);
}
// template <>
// AMREX_GPU_HOST_DEVICE auto RadSystem<TubeProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> quokka::valarray<double, nGroups_>
// {
// quokka::valarray<double, nGroups_> kappaPVec{};
// for (int g = 0; g < nGroups_; ++g) {
// kappaPVec[g] = kappa0;
// }
// return kappaPVec;
// }

// template <> AMREX_GPU_HOST_DEVICE auto RadSystem<TubeProblem>::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> quokka::valarray<double,
// nGroups_>
// {
// return ComputePlanckOpacity(rho, Tgas);
// }

// declare global variables
// initial conditions read from file
Expand Down Expand Up @@ -276,8 +294,8 @@ auto problem_main() -> int
// Problem initialization
RadhydroSimulation<TubeProblem> sim(BCs_cc);

sim.radiationReconstructionOrder_ = 2; // PLM
sim.reconstructionOrder_ = 2; // PLM
sim.radiationReconstructionOrder_ = 3; // PPM
sim.reconstructionOrder_ = 3; // PPM
sim.stopTime_ = tmax;
sim.cflNumber_ = CFL_number;
sim.radiationCflNumber_ = CFL_number;
Expand Down Expand Up @@ -323,27 +341,27 @@ auto problem_main() -> int
}
Erad_exact_arr[i] = Erad_0;
Erad_arr[i] = Erad_t;
const double Trad_exact = std::pow(Erad_0 / radiation_constant_cgs_, 1. / 4.);
const double Trad = std::pow(Erad_t / radiation_constant_cgs_, 1. / 4.);
const double Trad_exact = std::pow(Erad_0 / a_rad, 1. / 4.);
const double Trad = std::pow(Erad_t / a_rad, 1. / 4.);
Trad_arr[i] = Trad;
Trad_exact_arr[i] = Trad_exact;
Trad_err[i] = (Trad - Trad_exact) / Trad_exact;

double Egas_exact = values0.at(RadSystem<TubeProblem>::gasEnergy_index)[i];
double x1GasMom_exact = values0.at(RadSystem<TubeProblem>::x1GasMomentum_index)[i];
double x2GasMom_exact = values0.at(RadSystem<TubeProblem>::x2GasMomentum_index)[i];
double x3GasMom_exact = values0.at(RadSystem<TubeProblem>::x3GasMomentum_index)[i];
double const Egas_exact = values0.at(RadSystem<TubeProblem>::gasEnergy_index)[i];
double const x1GasMom_exact = values0.at(RadSystem<TubeProblem>::x1GasMomentum_index)[i];
double const x2GasMom_exact = values0.at(RadSystem<TubeProblem>::x2GasMomentum_index)[i];
double const x3GasMom_exact = values0.at(RadSystem<TubeProblem>::x3GasMomentum_index)[i];

double Egas = values.at(RadSystem<TubeProblem>::gasEnergy_index)[i];
double x1GasMom = values.at(RadSystem<TubeProblem>::x1GasMomentum_index)[i];
double x2GasMom = values.at(RadSystem<TubeProblem>::x2GasMomentum_index)[i];
double x3GasMom = values.at(RadSystem<TubeProblem>::x3GasMomentum_index)[i];
double const Egas = values.at(RadSystem<TubeProblem>::gasEnergy_index)[i];
double const x1GasMom = values.at(RadSystem<TubeProblem>::x1GasMomentum_index)[i];
double const x2GasMom = values.at(RadSystem<TubeProblem>::x2GasMomentum_index)[i];
double const x3GasMom = values.at(RadSystem<TubeProblem>::x3GasMomentum_index)[i];

double Eint_exact = RadSystem<TubeProblem>::ComputeEintFromEgas(rho_exact, x1GasMom_exact, x2GasMom_exact, x3GasMom_exact, Egas_exact);
double Tgas_exact = quokka::EOS<TubeProblem>::ComputeTgasFromEint(rho_exact, Eint_exact);
double const Eint_exact = RadSystem<TubeProblem>::ComputeEintFromEgas(rho_exact, x1GasMom_exact, x2GasMom_exact, x3GasMom_exact, Egas_exact);
double const Tgas_exact = quokka::EOS<TubeProblem>::ComputeTgasFromEint(rho_exact, Eint_exact);

double Eint = RadSystem<TubeProblem>::ComputeEintFromEgas(rho, x1GasMom, x2GasMom, x3GasMom, Egas);
double Tgas = quokka::EOS<TubeProblem>::ComputeTgasFromEint(rho, Eint);
double const Eint = RadSystem<TubeProblem>::ComputeEintFromEgas(rho, x1GasMom, x2GasMom, x3GasMom, Egas);
double const Tgas = quokka::EOS<TubeProblem>::ComputeTgasFromEint(rho, Eint);

Tgas_arr[i] = Tgas;
Tgas_err[i] = (Tgas - Tgas_exact) / Tgas_exact;
Expand Down
7 changes: 7 additions & 0 deletions src/RadhydroBB/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
add_executable(test_radhydro_bb test_radhydro_bb.cpp ../fextract.cpp ../interpolate.cpp ${QuokkaObjSources})

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

add_test(NAME RadhydroBB COMMAND test_radhydro_bb RadhydroBB.in WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
Loading

0 comments on commit 4686c6a

Please sign in to comment.