From 4c5af29b0b787305744dfebb9212605a7fea1976 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 11 Oct 2024 05:44:52 +0000 Subject: [PATCH] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../RadLineCooling/test_rad_line_cooling.cpp | 22 +++++++------- .../test_rad_line_cooling_MG.cpp | 28 +++++++++-------- src/radiation/radiation_dust_system.hpp | 30 +++++++++++-------- src/radiation/radiation_system.hpp | 24 ++++++++------- src/radiation/source_terms_multi_group.hpp | 7 +++-- src/radiation/source_terms_single_group.hpp | 15 +++++----- 6 files changed, 71 insertions(+), 55 deletions(-) diff --git a/src/problems/RadLineCooling/test_rad_line_cooling.cpp b/src/problems/RadLineCooling/test_rad_line_cooling.cpp index 62ede3d6c..4b8c65951 100644 --- a/src/problems/RadLineCooling/test_rad_line_cooling.cpp +++ b/src/problems/RadLineCooling/test_rad_line_cooling.cpp @@ -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; @@ -81,7 +81,8 @@ template <> struct ISM_Traits { }; template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRate(amrex::Real const temperature, amrex::Real const /*num_density*/) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRate(amrex::Real const temperature, amrex::Real const /*num_density*/) + -> quokka::valarray { quokka::valarray cooling{}; cooling.fillin(0.0); @@ -90,7 +91,8 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRate(amrex:: } template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) + -> quokka::valarray { quokka::valarray cooling{}; cooling.fillin(0.0); @@ -98,13 +100,12 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRateTempDeri return cooling; } -template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double { return CR_heating_rate; } -template <> +template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real { return kappa0; @@ -161,7 +162,7 @@ template <> void QuokkaSimulation::computeAfterTimestep() userData_.t_vec_.push_back(tNew_[0]); const amrex::Real Etot_i = values.at(RadSystem::gasEnergy_index)[0]; - const amrex::Real x1GasMom = values.at(RadSystem::x1GasMomentum_index)[0]; + const amrex::Real x1GasMom = values.at(RadSystem::x1GasMomentum_index)[0]; const amrex::Real x2GasMom = values.at(RadSystem::x2GasMomentum_index)[0]; const amrex::Real x3GasMom = values.at(RadSystem::x3GasMomentum_index)[0]; const amrex::Real rho = values.at(RadSystem::gasDensity_index)[0]; @@ -229,11 +230,12 @@ auto problem_main() -> int std::vector Tgas_exact_vec{}; std::vector 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); } diff --git a/src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp b/src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp index b1bc85963..ce3e8684f 100644 --- a/src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp +++ b/src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp @@ -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; @@ -85,14 +85,15 @@ template <> struct ISM_Traits { }; template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::DefinePhotoelectricHeatingE1Derivative(amrex::Real const /*temperature*/, - amrex::Real const /*num_density*/) -> amrex::Real +AMREX_GPU_HOST_DEVICE auto RadSystem::DefinePhotoelectricHeatingE1Derivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) + -> amrex::Real { return PE_rate / Erad_FUV; } template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRate(amrex::Real const temperature, amrex::Real const /*num_density*/) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRate(amrex::Real const temperature, amrex::Real const /*num_density*/) + -> quokka::valarray { quokka::valarray cooling{}; cooling.fillin(0.0); @@ -101,7 +102,8 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRate(amrex:: } template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) + -> quokka::valarray { quokka::valarray cooling{}; cooling.fillin(0.0); @@ -109,8 +111,7 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRateTempDeri return cooling; } -template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double { return CR_heating_rate; } @@ -232,17 +233,18 @@ auto problem_main() -> int const double heating_rate_ = ISM_Traits::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 t_exact_vec{}; std::vector Tgas_exact_vec{}; std::vector 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; @@ -300,8 +302,10 @@ auto problem_main() -> int std::vector Tgas_interp(t.size()); std::vector Erad_line_interp(t.size()); - interpolate_arrays(t.data(), Tgas_interp.data(), static_cast(t.size()), t_exact_vec.data(), Tgas_exact_vec.data(), static_cast(t_exact_vec.size())); - interpolate_arrays(t.data(), Erad_line_interp.data(), static_cast(t.size()), t_exact_vec.data(), Erad_line_exact_vec.data(), static_cast(t_exact_vec.size())); + interpolate_arrays(t.data(), Tgas_interp.data(), static_cast(t.size()), t_exact_vec.data(), Tgas_exact_vec.data(), + static_cast(t_exact_vec.size())); + interpolate_arrays(t.data(), Erad_line_interp.data(), static_cast(t.size()), t_exact_vec.data(), Erad_line_exact_vec.data(), + static_cast(t_exact_vec.size())); // compute error norm double err_norm = 0.; diff --git a/src/radiation/radiation_dust_system.hpp b/src/radiation/radiation_dust_system.hpp index 3bd9aaa40..fc5774b44 100644 --- a/src/radiation/radiation_dust_system.hpp +++ b/src/radiation/radiation_dust_system.hpp @@ -23,8 +23,8 @@ template AMREX_GPU_DEVICE auto RadSystem::ComputeJacobianForGasAndDust( double T_gas, double T_d, double Egas_diff, quokka::valarray const &Erad_diff, quokka::valarray const &Rvec, quokka::valarray const &Src, double coeff_n, quokka::valarray const &tau, double c_v, double /*lambda_gd_time_dt*/, - quokka::valarray const &kappaPoverE, quokka::valarray const &d_fourpiboverc_d_t, - const double num_den, const double dt) -> JacobianResult + quokka::valarray const &kappaPoverE, quokka::valarray const &d_fourpiboverc_d_t, const double num_den, const double dt) + -> JacobianResult { JacobianResult result; @@ -38,7 +38,7 @@ AMREX_GPU_DEVICE auto RadSystem::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) { @@ -61,7 +61,7 @@ AMREX_GPU_DEVICE auto RadSystem::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) { @@ -146,7 +146,7 @@ AMREX_GPU_DEVICE auto RadSystem::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) { @@ -180,7 +180,7 @@ AMREX_GPU_DEVICE auto RadSystem::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); @@ -541,7 +541,7 @@ AMREX_GPU_DEVICE auto RadSystem::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); @@ -789,8 +789,9 @@ AMREX_GPU_DEVICE auto RadSystem::SolveGasDustRadiationEnergyExchangeW JacobianResult 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); @@ -872,14 +873,17 @@ AMREX_GPU_DEVICE auto RadSystem::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::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::ComputeTgasFromEint(rho, Egas_, massScalars); const auto d_cooling_d_Tgas_ = DefineNetCoolingRateTempDerivative(T_gas_, num_den) * dt; @@ -890,7 +894,7 @@ AMREX_GPU_DEVICE auto RadSystem::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); diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index d924e1a04..8c5b2ef7f 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -349,7 +349,8 @@ template class RadSystem : public HyperbolicSystem quokka::valarray; - AMREX_GPU_HOST_DEVICE static auto DefineNetCoolingRateTempDerivative(amrex::Real temperature, amrex::Real num_density) -> quokka::valarray; + AMREX_GPU_HOST_DEVICE static auto DefineNetCoolingRateTempDerivative(amrex::Real temperature, amrex::Real num_density) + -> quokka::valarray; AMREX_GPU_HOST_DEVICE static auto DefineCosmicRayHeatingRate(amrex::Real num_density) -> double; @@ -368,14 +369,16 @@ template class RadSystem : public HyperbolicSystem const &Rvec, quokka::valarray const &Src, quokka::valarray const &tau, double c_v, quokka::valarray const &kappaPoverE, - quokka::valarray const &d_fourpiboverc_d_t, double num_den, double dt) -> JacobianResult; + quokka::valarray const &d_fourpiboverc_d_t, double num_den, double dt) + -> JacobianResult; AMREX_GPU_DEVICE static auto ComputeJacobianForGasAndDust(double T_gas, double T_d, double Egas_diff, quokka::valarray const &Erad_diff, quokka::valarray const &Rvec, quokka::valarray const &Src, double coeff_n, quokka::valarray const &tau, double c_v, double lambda_gd_time_dt, quokka::valarray const &kappaPoverE, - quokka::valarray const &d_fourpiboverc_d_t, double num_den, double dt) -> JacobianResult; + quokka::valarray const &d_fourpiboverc_d_t, double num_den, double dt) + -> JacobianResult; AMREX_GPU_DEVICE static auto ComputeJacobianForGasAndDustDecoupled( double T_gas, double T_d, double Egas_diff, quokka::valarray const &Erad_diff, quokka::valarray const &Rvec, @@ -507,15 +510,15 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiationTempDeri } // Define the background heating rate for the gas-dust-radiation system. Units in cgs: erg cm^-3 s^-1 -template -AMREX_GPU_HOST_DEVICE auto RadSystem::DefineBackgroundHeatingRate(amrex::Real const /*num_density*/) -> amrex::Real +template AMREX_GPU_HOST_DEVICE auto RadSystem::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 -AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRate(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRate(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) + -> quokka::valarray { quokka::valarray cooling{}; cooling.fillin(0.0); @@ -524,15 +527,15 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::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 -AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto RadSystem::DefineNetCoolingRateTempDerivative(amrex::Real const /*temperature*/, amrex::Real const /*num_density*/) + -> quokka::valarray { quokka::valarray cooling{}; cooling.fillin(0.0); return cooling; } -template -AMREX_GPU_HOST_DEVICE auto RadSystem::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double +template AMREX_GPU_HOST_DEVICE auto RadSystem::DefineCosmicRayHeatingRate(amrex::Real const /*num_density*/) -> double { return 0.0; } @@ -1440,7 +1443,8 @@ AMREX_GPU_DEVICE auto RadSystem::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; diff --git a/src/radiation/source_terms_multi_group.hpp b/src/radiation/source_terms_multi_group.hpp index 1bf42e3ec..379667870 100644 --- a/src/radiation/source_terms_multi_group.hpp +++ b/src/radiation/source_terms_multi_group.hpp @@ -107,8 +107,8 @@ AMREX_GPU_DEVICE auto RadSystem::ComputeJacobianForGas(double /*T_d*/ quokka::valarray const &Rvec, quokka::valarray const &Src, quokka::valarray const &tau, double c_v, quokka::valarray const &kappaPoverE, - quokka::valarray const &d_fourpiboverc_d_t, - double const num_den, double const dt) -> JacobianResult + quokka::valarray const &d_fourpiboverc_d_t, double const num_den, + double const dt) -> JacobianResult { JacobianResult result; @@ -311,7 +311,8 @@ AMREX_GPU_DEVICE auto RadSystem::SolveGasRadiationEnergyExchange( const auto Egas_diff = Egas_guess - Egas0; const auto Erad_diff = EradVec_guess - Erad0Vec; - auto jacobian = ComputeJacobianForGas(T_d, Egas_diff, Erad_diff, Rvec, Src, tau, c_v, opacity_terms.kappaPoverE, d_fourpiboverc_d_t, num_den, dt); + auto jacobian = + ComputeJacobianForGas(T_d, Egas_diff, Erad_diff, Rvec, Src, tau, c_v, opacity_terms.kappaPoverE, d_fourpiboverc_d_t, num_den, dt); if constexpr (use_D_as_base) { jacobian.J0g = jacobian.J0g * tau0; diff --git a/src/radiation/source_terms_single_group.hpp b/src/radiation/source_terms_single_group.hpp index a4a62631a..694c33a8b 100644 --- a/src/radiation/source_terms_single_group.hpp +++ b/src/radiation/source_terms_single_group.hpp @@ -228,8 +228,8 @@ void RadSystem::AddSourceTermsSingleGroup(array_t &consVar, arraycons } } - double cooling = 0.0; - double cooling_derivative = 0.0; + double cooling = 0.0; + double cooling_derivative = 0.0; const double CR_heating = DefineCosmicRayHeatingRate(num_den) * dt; if constexpr (enable_dust_gas_thermal_coupling_model_) { cooling = DefineNetCoolingRate(T_gas, num_den)[0]; @@ -239,7 +239,7 @@ void RadSystem::AddSourceTermsSingleGroup(array_t &consVar, arraycons F_G = Egas_guess - Egas0 + cscale * R + cooling * dt - CR_heating; F_D = Erad_guess - Erad0 - (R + Src); if constexpr (add_line_cooling_to_radiation_in_jac) { - F_D -= (1.0/cscale) * cooling * dt; + F_D -= (1.0 / cscale) * cooling * dt; } double F_D_abs = 0.0; if (tau > 0.0) { @@ -284,7 +284,7 @@ void RadSystem::AddSourceTermsSingleGroup(array_t &consVar, arraycons if constexpr (!enable_dust_gas_thermal_coupling_model_) { J00 = 1.0 + cooling_derivative * dt / c_v; J01 = cscale; - J10 = 1.0 / c_v * dEg_dT - (1/cscale) * cooling_derivative * dt; + J10 = 1.0 / c_v * dEg_dT - (1 / cscale) * cooling_derivative * dt; if (tau <= 0.0) { J11 = -std::numeric_limits::infinity(); } else { @@ -299,7 +299,7 @@ void RadSystem::AddSourceTermsSingleGroup(array_t &consVar, arraycons J01 = cscale; J10 = 1.0 / c_v * dEg_dT; if (tau <= 0.0) { - J11 = - LARGE; + J11 = -LARGE; } else { J11 = kappaPoverE * d_fourpiboverc_d_t * dTd_dRg - kappaPoverE / tau - 1.0; } @@ -350,9 +350,10 @@ void RadSystem::AddSourceTermsSingleGroup(array_t &consVar, arraycons if constexpr (!add_line_cooling_to_radiation_in_jac) { const auto cooling_tend = DefineNetCoolingRate(T_gas, num_den)[0] * dt; - AMREX_ASSERT_WITH_MESSAGE(cooling_tend >= 0., "add_line_cooling_to_radiation has to be enabled when there is negative cooling rate!"); + AMREX_ASSERT_WITH_MESSAGE(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. - Erad_guess += (1/cscale) * cooling_tend; + Erad_guess += (1 / cscale) * cooling_tend; } if (n > 0) {