Skip to content

Commit

Permalink
[pre-commit.ci] auto fixes from pre-commit.com hooks
Browse files Browse the repository at this point in the history 
for more information, see https://pre-commit.ci
  • Loading branch information
@pre-commit-ci
pre-commit-ci[bot] committed Oct 11, 2024
1 parent a8df3b9 commit 4c5af29
Show file tree
Hide file tree
Showing 6 changed files with 71 additions and 55 deletions.
22 changes: 12 additions & 10 deletions src/problems/RadLineCooling/test_rad_line_cooling.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,10 @@

#include "AMReX_Array.H"
#include "AMReX_BC_TYPES.H"
#include "util/fextract.hpp"
#include "AMReX_Print.H"
#include "physics_info.hpp"
#include "test_rad_line_cooling.hpp"
#include "util/fextract.hpp"

static constexpr bool export_csv = true;

Expand Down Expand Up @@ -81,7 +81,8 @@ template <> struct ISM_Traits<PulseProblem> {
};

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineNetCoolingRate(amrex::Real const temperature, amrex::Real const /*num_density*/) -> quokka::valarray<double, nGroups_>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineNetCoolingRate(amrex::Real const temperature, amrex::Real const /*num_density*/)
-> quokka::valarray<double, nGroups_>
{
quokka::valarray<double, nGroups_> cooling{};
cooling.fillin(0.0);
Expand All @@ -90,21 +91,21 @@ AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineNetCoolingRate(amrex::
}

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) -> quokka::valarray<double, nGroups_>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/)
-> quokka::valarray<double, nGroups_>
{
quokka::valarray<double, nGroups_> cooling{};
cooling.fillin(0.0);
cooling[0] = cooling_rate;
return cooling;
}

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double
template <> AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double
{
return CR_heating_rate;
}

template <>
template <>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem<PulseProblem>::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real
{
return kappa0;
Expand Down Expand Up @@ -161,7 +162,7 @@ template <> void QuokkaSimulation<PulseProblem>::computeAfterTimestep()
userData_.t_vec_.push_back(tNew_[0]);

const amrex::Real Etot_i = values.at(RadSystem<PulseProblem>::gasEnergy_index)[0];
const amrex::Real x1GasMom = values.at(RadSystem<PulseProblem>::x1GasMomentum_index)[0];
const amrex::Real x1GasMom = values.at(RadSystem<PulseProblem>::x1GasMomentum_index)[0];
const amrex::Real x2GasMom = values.at(RadSystem<PulseProblem>::x2GasMomentum_index)[0];
const amrex::Real x3GasMom = values.at(RadSystem<PulseProblem>::x3GasMomentum_index)[0];
const amrex::Real rho = values.at(RadSystem<PulseProblem>::gasDensity_index)[0];
Expand Down Expand Up @@ -229,11 +230,12 @@ auto problem_main() -> int

std::vector<double> Tgas_exact_vec{};
std::vector<double> Erad_exact_vec{};
for (const auto& t_exact : t_sim) {
const double Egas_exact_solution = std::exp(-cooling_rate * t_exact) * (cooling_rate * T0 - CR_heating_rate + CR_heating_rate * std::exp(cooling_rate * t_exact)) / cooling_rate;
for (const auto &t_exact : t_sim) {
const double Egas_exact_solution = std::exp(-cooling_rate * t_exact) *
(cooling_rate * T0 - CR_heating_rate + CR_heating_rate * std::exp(cooling_rate * t_exact)) / cooling_rate;
const double T_exact_solution = Egas_exact_solution / C_V;
Tgas_exact_vec.push_back(T_exact_solution);
const double Erad_exact_solution = - (Egas_exact_solution - C_V * T0 - CR_heating_rate * t_exact) * (chat / c);
const double Erad_exact_solution = -(Egas_exact_solution - C_V * T0 - CR_heating_rate * t_exact) * (chat / c);
Erad_exact_vec.push_back(Erad_exact_solution);
}

Expand Down
28 changes: 16 additions & 12 deletions src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,10 @@

#include "AMReX_Array.H"
#include "AMReX_BC_TYPES.H"
#include "util/fextract.hpp"
#include "AMReX_Print.H"
#include "physics_info.hpp"
#include "test_rad_line_cooling_MG.hpp"
#include "util/fextract.hpp"

static constexpr bool export_csv = true;

Expand Down Expand Up @@ -85,14 +85,15 @@ template <> struct ISM_Traits<PulseProblem> {
};

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefinePhotoelectricHeatingE1Derivative(amrex::Real const /*temperature*/,
amrex::Real const /*num_density*/) -> amrex::Real
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefinePhotoelectricHeatingE1Derivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/)
-> amrex::Real
{
return PE_rate / Erad_FUV;
}

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineNetCoolingRate(amrex::Real const temperature, amrex::Real const /*num_density*/) -> quokka::valarray<double, nGroups_>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineNetCoolingRate(amrex::Real const temperature, amrex::Real const /*num_density*/)
-> quokka::valarray<double, nGroups_>
{
quokka::valarray<double, nGroups_> cooling{};
cooling.fillin(0.0);
Expand All @@ -101,16 +102,16 @@ AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineNetCoolingRate(amrex::
}

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) -> quokka::valarray<double, nGroups_>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/)
-> quokka::valarray<double, nGroups_>
{
quokka::valarray<double, nGroups_> cooling{};
cooling.fillin(0.0);
cooling[line_index] = cooling_rate;
return cooling;
}

template <>
AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double
template <> AMREX_GPU_HOST_DEVICE auto RadSystem<PulseProblem>::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double
{
return CR_heating_rate;
}
Expand Down Expand Up @@ -232,17 +233,18 @@ auto problem_main() -> int

const double heating_rate_ = ISM_Traits<PulseProblem>::enable_photoelectric_heating ? PE_rate + CR_heating_rate : CR_heating_rate;

// compute exact solution from t = 0 to t = t_end
// compute exact solution from t = 0 to t = t_end
const double N_dt = 1000.;
double t_exact = 0.0;
std::vector<double> t_exact_vec{};
std::vector<double> Tgas_exact_vec{};
std::vector<double> Erad_line_exact_vec{};
while (true) {
const double Egas_exact_solution = std::exp(-cooling_rate * t_exact) * (cooling_rate * T0 - heating_rate_ + heating_rate_ * std::exp(cooling_rate * t_exact)) / cooling_rate;
const double Egas_exact_solution =
std::exp(-cooling_rate * t_exact) * (cooling_rate * T0 - heating_rate_ + heating_rate_ * std::exp(cooling_rate * t_exact)) / cooling_rate;
const double T_exact_solution = Egas_exact_solution / C_V;
Tgas_exact_vec.push_back(T_exact_solution);
const double Erad_line_exact_solution = - (Egas_exact_solution - C_V * T0 - heating_rate_ * t_exact) * (chat / c);
const double Erad_line_exact_solution = -(Egas_exact_solution - C_V * T0 - heating_rate_ * t_exact) * (chat / c);
Erad_line_exact_vec.push_back(Erad_line_exact_solution);
t_exact_vec.push_back(t_exact);
t_exact += t_end / N_dt;
Expand Down Expand Up @@ -300,8 +302,10 @@ auto problem_main() -> int

std::vector<double> Tgas_interp(t.size());
std::vector<double> Erad_line_interp(t.size());
interpolate_arrays(t.data(), Tgas_interp.data(), static_cast<int>(t.size()), t_exact_vec.data(), Tgas_exact_vec.data(), static_cast<int>(t_exact_vec.size()));
interpolate_arrays(t.data(), Erad_line_interp.data(), static_cast<int>(t.size()), t_exact_vec.data(), Erad_line_exact_vec.data(), static_cast<int>(t_exact_vec.size()));
interpolate_arrays(t.data(), Tgas_interp.data(), static_cast<int>(t.size()), t_exact_vec.data(), Tgas_exact_vec.data(),
static_cast<int>(t_exact_vec.size()));
interpolate_arrays(t.data(), Erad_line_interp.data(), static_cast<int>(t.size()), t_exact_vec.data(), Erad_line_exact_vec.data(),
static_cast<int>(t_exact_vec.size()));

// compute error norm
double err_norm = 0.;
Expand Down
30 changes: 17 additions & 13 deletions src/radiation/radiation_dust_system.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,8 @@ template <typename problem_t>
AMREX_GPU_DEVICE auto RadSystem<problem_t>::ComputeJacobianForGasAndDust(
double T_gas, double T_d, double Egas_diff, quokka::valarray<double, nGroups_> const &Erad_diff, quokka::valarray<double, nGroups_> const &Rvec,
quokka::valarray<double, nGroups_> const &Src, double coeff_n, quokka::valarray<double, nGroups_> const &tau, double c_v, double /*lambda_gd_time_dt*/,
quokka::valarray<double, nGroups_> const &kappaPoverE, quokka::valarray<double, nGroups_> const &d_fourpiboverc_d_t,
const double num_den, const double dt) -> JacobianResult<problem_t>
quokka::valarray<double, nGroups_> const &kappaPoverE, quokka::valarray<double, nGroups_> const &d_fourpiboverc_d_t, const double num_den, const double dt)
-> JacobianResult<problem_t>
{
JacobianResult<problem_t> result;

Expand All @@ -38,7 +38,7 @@ AMREX_GPU_DEVICE auto RadSystem<problem_t>::ComputeJacobianForGasAndDust(
result.F0 = Egas_diff + cscale * sum(Rvec) + sum(cooling) - CR_heating;
result.Fg = Erad_diff - (Rvec + Src);
if constexpr (add_line_cooling_to_radiation_in_jac) {
result.Fg -= (1.0/cscale) * cooling;
result.Fg -= (1.0 / cscale) * cooling;
}
result.Fg_abs_sum = 0.0;
for (int g = 0; g < nGroups_; ++g) {
Expand All @@ -61,7 +61,7 @@ AMREX_GPU_DEVICE auto RadSystem<problem_t>::ComputeJacobianForGasAndDust(
dEg_dT *= d_Td_d_T;
const double dTd_dRg = -1.0 / (coeff_n * std::sqrt(T_gas));
const auto rg = kappaPoverE * d_fourpiboverc_d_t * dTd_dRg;
result.Jg0 = 1.0 / c_v * dEg_dT - (1/cscale) * cooling_derivative - 1.0 / cscale * rg * result.J00;
result.Jg0 = 1.0 / c_v * dEg_dT - (1 / cscale) * cooling_derivative - 1.0 / cscale * rg * result.J00;
// Note that Fg is modified here, but it does not change Fg_abs_sum, which is used to check the convergence.
result.Fg = result.Fg - 1.0 / cscale * rg * result.F0;
for (int g = 0; g < nGroups_; ++g) {
Expand Down Expand Up @@ -146,7 +146,7 @@ AMREX_GPU_DEVICE auto RadSystem<problem_t>::ComputeJacobianForGasAndDustWithPE(
result.F0 = Egas_diff + cscale * sum(Rvec) + sum(cooling) - PE_heating_energy_derivative * Erad[nGroups_ - 1] - CR_heating;
result.Fg = Erad - Erad0 - (Rvec + Src);
if constexpr (add_line_cooling_to_radiation_in_jac) {
result.Fg -= (1.0/cscale) * cooling;
result.Fg -= (1.0 / cscale) * cooling;
}
result.Fg_abs_sum = 0.0;
for (int g = 0; g < nGroups_; ++g) {
Expand Down Expand Up @@ -180,7 +180,7 @@ AMREX_GPU_DEVICE auto RadSystem<problem_t>::ComputeJacobianForGasAndDustWithPE(
const auto dEg_dT = kappaPoverE * d_fourpiboverc_d_t * d_Td_d_T;
const double dTd_dRg = -1.0 / (coeff_n * std::sqrt(T_gas));
const auto rg = kappaPoverE * d_fourpiboverc_d_t * dTd_dRg;
result.Jg0 = 1.0 / c_v * dEg_dT - (1/cscale) * cooling_derivative - 1.0 / cscale * rg * result.J00;
result.Jg0 = 1.0 / c_v * dEg_dT - (1 / cscale) * cooling_derivative - 1.0 / cscale * rg * result.J00;
// Note that Fg is modified here, but it does not change Fg_abs_sum, which is used to check the convergence.
result.Fg = result.Fg - 1.0 / cscale * rg * result.F0;
result.Jgg = d_Eg_d_Rg + (-1.0);
Expand Down Expand Up @@ -541,7 +541,7 @@ AMREX_GPU_DEVICE auto RadSystem<problem_t>::SolveGasDustRadiationEnergyExchange(
if constexpr (!add_line_cooling_to_radiation_in_jac) {
AMREX_ASSERT_WITH_MESSAGE(min(cooling_tend) >= 0., "add_line_cooling_to_radiation has to be enabled when there is negative cooling rate!");
// TODO(CCH): potential GPU-related issue here.
EradVec_guess += (1/cscale) * cooling_tend;
EradVec_guess += (1 / cscale) * cooling_tend;
}

AMREX_ASSERT(Egas_guess > 0.0);
Expand Down Expand Up @@ -789,8 +789,9 @@ AMREX_GPU_DEVICE auto RadSystem<problem_t>::SolveGasDustRadiationEnergyExchangeW
JacobianResult<problem_t> jacobian;

if (dust_model == 1) {
jacobian = ComputeJacobianForGasAndDustWithPE(T_gas, T_d, Egas_diff, EradVec_guess, Erad0Vec, PE_heating_energy_derivative, Rvec, Src,
coeff_n, tau, c_v, lambda_gd_times_dt, opacity_terms.kappaPoverE, d_fourpiboverc_d_t, num_den, dt);
jacobian =
ComputeJacobianForGasAndDustWithPE(T_gas, T_d, Egas_diff, EradVec_guess, Erad0Vec, PE_heating_energy_derivative, Rvec, Src, coeff_n,
tau, c_v, lambda_gd_times_dt, opacity_terms.kappaPoverE, d_fourpiboverc_d_t, num_den, dt);
} else {
jacobian = ComputeJacobianForGasAndDustDecoupled(T_gas, T_d, Egas_diff, Erad_diff, Rvec, Src, coeff_n, tau, c_v, lambda_gd_times_dt,
opacity_terms.kappaPoverE, d_fourpiboverc_d_t);
Expand Down Expand Up @@ -872,14 +873,17 @@ AMREX_GPU_DEVICE auto RadSystem<problem_t>::SolveGasDustRadiationEnergyExchangeW
const double CR_heating = DefineCosmicRayHeatingRate(num_den) * dt;
const double compare = Egas_guess + cscale * lambda_gd_times_dt + sum(abs(cooling_tend)) + CR_heating;

// RHS of the equation 0 = Egas - Egas0 + cscale * lambda_gd_times_dt + sum(cooling) - PE_heating_energy_derivative * EradVec_guess[nGroups_ - 1];
// RHS of the equation 0 = Egas - Egas0 + cscale * lambda_gd_times_dt + sum(cooling) - PE_heating_energy_derivative * EradVec_guess[nGroups_ -
// 1];
auto rhs = [=](double Egas_) -> double {
const double T_gas_ = quokka::EOS<problem_t>::ComputeTgasFromEint(rho, Egas_, massScalars);
const auto cooling_ = DefineNetCoolingRate(T_gas_, num_den) * dt;
return Egas_ - Egas0 + cscale * lambda_gd_times_dt + sum(cooling_) - PE_heating_energy_derivative * EradVec_guess[nGroups_ - 1] - CR_heating;
return Egas_ - Egas0 + cscale * lambda_gd_times_dt + sum(cooling_) - PE_heating_energy_derivative * EradVec_guess[nGroups_ - 1] -
CR_heating;
};

// Jacobian of the RHS of the equation 0 = Egas - Egas0 + cscale * lambda_gd_times_dt + sum(cooling) + PE_heating_energy_derivative * EradVec_guess[nGroups_ - 1];
// Jacobian of the RHS of the equation 0 = Egas - Egas0 + cscale * lambda_gd_times_dt + sum(cooling) + PE_heating_energy_derivative *
// EradVec_guess[nGroups_ - 1];
auto jac = [=](double Egas_) -> double {
const double T_gas_ = quokka::EOS<problem_t>::ComputeTgasFromEint(rho, Egas_, massScalars);
const auto d_cooling_d_Tgas_ = DefineNetCoolingRateTempDerivative(T_gas_, num_den) * dt;
Expand All @@ -890,7 +894,7 @@ AMREX_GPU_DEVICE auto RadSystem<problem_t>::SolveGasDustRadiationEnergyExchangeW
}

if constexpr (!add_line_cooling_to_radiation_in_jac) {
EradVec_guess += (1/cscale) * cooling_tend;
EradVec_guess += (1 / cscale) * cooling_tend;
}

AMREX_ASSERT(Egas_guess > 0.0);
Expand Down
24 changes: 14 additions & 10 deletions src/radiation/radiation_system.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -349,7 +349,8 @@ template <typename problem_t> class RadSystem : public HyperbolicSystem<problem_

AMREX_GPU_HOST_DEVICE static auto DefineNetCoolingRate(amrex::Real temperature, amrex::Real num_density) -> quokka::valarray<double, nGroups_>;

AMREX_GPU_HOST_DEVICE static auto DefineNetCoolingRateTempDerivative(amrex::Real temperature, amrex::Real num_density) -> quokka::valarray<double, nGroups_>;
AMREX_GPU_HOST_DEVICE static auto DefineNetCoolingRateTempDerivative(amrex::Real temperature, amrex::Real num_density)
-> quokka::valarray<double, nGroups_>;

AMREX_GPU_HOST_DEVICE static auto DefineCosmicRayHeatingRate(amrex::Real num_density) -> double;

Expand All @@ -368,14 +369,16 @@ template <typename problem_t> class RadSystem : public HyperbolicSystem<problem_
quokka::valarray<double, nGroups_> const &Rvec, quokka::valarray<double, nGroups_> const &Src,
quokka::valarray<double, nGroups_> const &tau, double c_v,
quokka::valarray<double, nGroups_> const &kappaPoverE,
quokka::valarray<double, nGroups_> const &d_fourpiboverc_d_t, double num_den, double dt) -> JacobianResult<problem_t>;
quokka::valarray<double, nGroups_> const &d_fourpiboverc_d_t, double num_den, double dt)
-> JacobianResult<problem_t>;

AMREX_GPU_DEVICE static auto ComputeJacobianForGasAndDust(double T_gas, double T_d, double Egas_diff,
quokka::valarray<double, nGroups_> const &Erad_diff,
quokka::valarray<double, nGroups_> const &Rvec, quokka::valarray<double, nGroups_> const &Src,
double coeff_n, quokka::valarray<double, nGroups_> const &tau, double c_v,
double lambda_gd_time_dt, quokka::valarray<double, nGroups_> const &kappaPoverE,
quokka::valarray<double, nGroups_> const &d_fourpiboverc_d_t, double num_den, double dt) -> JacobianResult<problem_t>;
quokka::valarray<double, nGroups_> const &d_fourpiboverc_d_t, double num_den, double dt)
-> JacobianResult<problem_t>;

AMREX_GPU_DEVICE static auto ComputeJacobianForGasAndDustDecoupled(
double T_gas, double T_d, double Egas_diff, quokka::valarray<double, nGroups_> const &Erad_diff, quokka::valarray<double, nGroups_> const &Rvec,
Expand Down Expand Up @@ -507,15 +510,15 @@ AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::ComputeThermalRadiationTempDeri
}

// Define the background heating rate for the gas-dust-radiation system. Units in cgs: erg cm^-3 s^-1
template <typename problem_t>
AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::DefineBackgroundHeatingRate(amrex::Real const /*num_density*/) -> amrex::Real
template <typename problem_t> AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::DefineBackgroundHeatingRate(amrex::Real const /*num_density*/) -> amrex::Real
{
return 0.0;
}

// Define the net cooling rate (line cooling + heating) for the gas-dust-radiation system. Units in cgs: erg cm^-3 s^-1
template <typename problem_t>
AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::DefineNetCoolingRate(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) -> quokka::valarray<double, nGroups_>
AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::DefineNetCoolingRate(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/)
-> quokka::valarray<double, nGroups_>
{
quokka::valarray<double, nGroups_> cooling{};
cooling.fillin(0.0);
Expand All @@ -524,15 +527,15 @@ AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::DefineNetCoolingRate(amrex::Rea

// Define the derivative of the net cooling rate with respect to temperature for the gas-dust-radiation system. Units in cgs: erg cm^-3 s^-1 K^-1
template <typename problem_t>
AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) -> quokka::valarray<double, nGroups_>
AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/)
-> quokka::valarray<double, nGroups_>
{
quokka::valarray<double, nGroups_> cooling{};
cooling.fillin(0.0);
return cooling;
}

template <typename problem_t>
AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double
template <typename problem_t> AMREX_GPU_HOST_DEVICE auto RadSystem<problem_t>::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double
{
return 0.0;
}
Expand Down Expand Up @@ -1440,7 +1443,8 @@ AMREX_GPU_DEVICE auto RadSystem<problem_t>::ComputeDustTemperatureBateKeto(doubl
} else {
const auto fourPiBoverC = ComputeThermalRadiationMultiGroup(T_d, rad_boundaries);
const auto opacity_terms = ComputeModelDependentKappaEAndKappaP(T_d, rho, rad_boundaries, rad_boundary_ratios, fourPiBoverC, Erad, 0);
LHS = c_hat_ * dt * rho * sum(opacity_terms.kappaE * Erad - opacity_terms.kappaP * fourPiBoverC) + N_d * std::sqrt(T_gas) * (T_gas - T_d);
LHS =
c_hat_ * dt * rho * sum(opacity_terms.kappaE * Erad - opacity_terms.kappaP * fourPiBoverC) + N_d * std::sqrt(T_gas) * (T_gas - T_d);
}

return LHS;
Expand Down
Loading

0 comments on commit 4c5af29

Please sign in to comment.