From c4f575dee97cad0f0e2bfc7b6a8b36496c6c6d07 Mon Sep 17 00:00:00 2001 From: chongchonghe Date: Thu, 8 Aug 2024 12:14:38 +1000 Subject: [PATCH 01/25] mod gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index bf8e1ceb9..9efbe8521 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,8 @@ __pycache__ *.code-workspace *.sarif +/docs/site/ + *.aux *.dbl *.bbl From a332f4aea60ee4aef6c22910c8233227da1ae955 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 15 Aug 2024 14:55:29 +1000 Subject: [PATCH 02/25] mod gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 9efbe8521..9de410af7 100644 --- a/.gitignore +++ b/.gitignore @@ -17,6 +17,7 @@ sims/ __pycache__ *.code-workspace *.sarif +Backtrace.0 /docs/site/ From eb0ba648a5570085fabac24b3154983a5f6dddd6 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 15 Aug 2024 18:23:05 +1000 Subject: [PATCH 03/25] cleanup: remove repeating calculation of kappa --- src/radiation/radiation_system.hpp | 470 ++++++++++++----------------- 1 file changed, 201 insertions(+), 269 deletions(-) diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 32d154172..4e1e63226 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -1344,131 +1344,161 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // dF_{D,i} / dE_g = 1 / (chat * C_v) * (kappa_{P,i} / kappa_{E,i}) * d/dT (4 \pi B_i) // dF_{D,i} / dD_i = - (1 / (chat * dt * rho * kappa_{E,i}) + 1) * tau0_i = - ((1 / tau_i)(kappa_Pi / kappa_Ei) + 1) * tau0_i - T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas0, massScalars); - AMREX_ASSERT(T_gas >= 0.); - fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); + double F_G = NAN; + double dFG_dEgas = NAN; + double deltaEgas = NAN; + quokka::valarray dFG_dD{}; + quokka::valarray dFR_dEgas{}; + quokka::valarray dFR_i_dD_i{}; + quokka::valarray deltaD{}; + quokka::valarray F_D{}; - if constexpr (opacity_model_ == OpacityModel::single_group) { - kappaPVec[0] = ComputePlanckOpacity(rho, T_gas); - kappaEVec[0] = ComputeEnergyMeanOpacity(rho, T_gas); - kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_gas); - } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { - kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); - for (int g = 0; g < nGroups_; ++g) { - kappaPVec[g] = kappa_expo_and_lower_value[1][g]; - kappaEVec[g] = kappa_expo_and_lower_value[1][g]; - kappaFVec[g] = kappa_expo_and_lower_value[1][g]; + const double resid_tol = 1.0e-11; // 1.0e-15; + const int maxIter = 400; + int n = 0; + for (; n < maxIter; ++n) { + T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); + AMREX_ASSERT(T_gas >= 0.); + fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); + + if constexpr (opacity_model_ == OpacityModel::single_group) { + kappaPVec[0] = ComputePlanckOpacity(rho, T_gas); + kappaEVec[0] = ComputeEnergyMeanOpacity(rho, T_gas); + } else { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + for (int g = 0; g < nGroups_; ++g) { + auto const nu_L = radBoundaries_g_copy[g]; + auto const nu_R = radBoundaries_g_copy[g + 1]; + auto const B_L = PlanckFunction(nu_L, T_gas); // 4 pi B(nu) / c + auto const B_R = PlanckFunction(nu_R, T_gas); // 4 pi B(nu) / c + auto const kappa_L = kappa_expo_and_lower_value[1][g]; + auto const kappa_R = kappa_L * std::pow(nu_R / nu_L, kappa_expo_and_lower_value[0][g]); + delta_nu_kappa_B_at_edge[g] = nu_R * kappa_R * B_R - nu_L * kappa_L * B_L; + delta_nu_B_at_edge[g] = nu_R * B_R - nu_L * B_L; + } + if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + for (int g = 0; g < nGroups_; ++g) { + kappaPVec[g] = kappa_expo_and_lower_value[1][g]; + kappaEVec[g] = kappa_expo_and_lower_value[1][g]; + } + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); + kappaEVec = kappaPVec; + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + if (n < max_ite_to_update_alpha_E) { + alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); + alpha_E = ComputeRadQuantityExponents(EradVec_guess, radBoundaries_g_copy); + } + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); + kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); + } } - } else { - kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + AMREX_ASSERT(!kappaPVec.hasnan()); + AMREX_ASSERT(!kappaEVec.hasnan()); for (int g = 0; g < nGroups_; ++g) { - auto const nu_L = radBoundaries_g_copy[g]; - auto const nu_R = radBoundaries_g_copy[g + 1]; - auto const B_L = PlanckFunction(nu_L, T_gas); // 4 pi B(nu) / c - auto const B_R = PlanckFunction(nu_R, T_gas); // 4 pi B(nu) / c - auto const kappa_L = kappa_expo_and_lower_value[1][g]; - auto const kappa_R = kappa_L * std::pow(nu_R / nu_L, kappa_expo_and_lower_value[0][g]); - delta_nu_kappa_B_at_edge[g] = nu_R * kappa_R * B_R - nu_L * kappa_L * B_L; - delta_nu_B_at_edge[g] = nu_R * B_R - nu_L * B_L; - } - if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { - kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); - kappaEVec = kappaPVec; - if constexpr (use_diffuse_flux_mean_opacity) { - kappaFVec = - ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, delta_nu_kappa_B_at_edge, - delta_nu_B_at_edge, kappa_expo_and_lower_value[0]); + if (kappaEVec[g] > 0.0) { + kappaPoverE[g] = kappaPVec[g] / kappaEVec[g]; } else { - kappaFVec = kappaPVec; - } - } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { - alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); - alpha_E = ComputeRadQuantityExponents(Erad0Vec, radBoundaries_g_copy); - kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); - kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); - if constexpr (use_diffuse_flux_mean_opacity) { - kappaFVec = - ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, delta_nu_kappa_B_at_edge, - delta_nu_B_at_edge, kappa_expo_and_lower_value[0]); - } else { // fall back to bin-center opacity. For simplicity we don't do PPL on the flux mean opacity because it - // requires fitting to each components of the flux. - kappaFVec = ComputeBinCenterOpacity(radBoundaries_g_copy, kappa_expo_and_lower_value); + kappaPoverE[g] = 1.0; } } - } - AMREX_ASSERT(!kappaPVec.hasnan()); - AMREX_ASSERT(!kappaEVec.hasnan()); - AMREX_ASSERT(!kappaFVec.hasnan()); - - for (int g = 0; g < nGroups_; ++g) { - if (kappaEVec[g] > 0.0) { - kappaPoverE[g] = kappaPVec[g] / kappaEVec[g]; - } else { - kappaPoverE[g] = 1.0; - } - } - if constexpr ((beta_order_ != 0) && (include_work_term_in_source)) { - // compute the work term at the old state - // const double gamma = 1.0 / sqrt(1.0 - vsqr / (c * c)); - if (ite == 0) { + // In the first loop, calculate kappaF, work, tau0, R + if (n == 0) { if constexpr (opacity_model_ == OpacityModel::single_group) { - const double frad0 = consPrev(i, j, k, x1RadFlux_index); - const double frad1 = consPrev(i, j, k, x2RadFlux_index); - const double frad2 = consPrev(i, j, k, x3RadFlux_index); - // work = v * F * chi - work[0] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * - (2.0 * kappaEVec[0] - kappaFVec[0]) * chat / (c * c) * lorentz_factor_v * dt; - } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_gas); + } else { for (int g = 0; g < nGroups_; ++g) { - const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); - const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); - const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); - // work = v * F * chi - work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * kappaFVec[g] * chat / - (c * c) * dt; + auto const nu_L = radBoundaries_g_copy[g]; + auto const nu_R = radBoundaries_g_copy[g + 1]; + auto const B_L = PlanckFunction(nu_L, T_gas); // 4 pi B(nu) / c + auto const B_R = PlanckFunction(nu_R, T_gas); // 4 pi B(nu) / c + auto const kappa_L = kappa_expo_and_lower_value[1][g]; + auto const kappa_R = kappa_L * std::pow(nu_R / nu_L, kappa_expo_and_lower_value[0][g]); + delta_nu_kappa_B_at_edge[g] = nu_R * kappa_R * B_R - nu_L * kappa_L * B_L; + delta_nu_B_at_edge[g] = nu_R * B_R - nu_L * B_L; } - } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum || - opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { - for (int g = 0; g < nGroups_; ++g) { - const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); - const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); - const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); - // work = v * F * chi - work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * - (1.0 + kappa_expo_and_lower_value[0][g]) * kappaFVec[g] * chat / (c * c) * dt; + if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + kappaFVec = kappaPVec; + } else { + if constexpr (use_diffuse_flux_mean_opacity) { + kappaFVec = + ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, delta_nu_kappa_B_at_edge, + delta_nu_B_at_edge, kappa_expo_and_lower_value[0]); + } else { + // for simplicity, I assume kappaF = kappaE when opacity_model_ == OpacityModel::PPL_opacity_full_spectrum, if !use_diffuse_flux_mean_opacity. + // We won't use this option anyway. + kappaFVec = kappaEVec; + } + } + } + AMREX_ASSERT(!kappaFVec.hasnan()); + + if constexpr ((beta_order_ != 0) && (include_work_term_in_source)) { + // compute the work term at the old state + // const double gamma = 1.0 / sqrt(1.0 - vsqr / (c * c)); + if (ite == 0) { + if constexpr (opacity_model_ == OpacityModel::single_group) { + const double frad0 = consPrev(i, j, k, x1RadFlux_index); + const double frad1 = consPrev(i, j, k, x2RadFlux_index); + const double frad2 = consPrev(i, j, k, x3RadFlux_index); + // work = v * F * chi + work[0] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * + (2.0 * kappaEVec[0] - kappaFVec[0]) * chat / (c * c) * lorentz_factor_v * dt; + } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + for (int g = 0; g < nGroups_; ++g) { + const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); + const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); + const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); + // work = v * F * chi + work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * kappaFVec[g] * chat / + (c * c) * dt; + } + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum || + opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + for (int g = 0; g < nGroups_; ++g) { + const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); + const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); + const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); + // work = v * F * chi + work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * + (1.0 + kappa_expo_and_lower_value[0][g]) * kappaFVec[g] * chat / (c * c) * dt; + } + } } } - } - } - tau0 = dt * rho * kappaPVec * chat * lorentz_factor; - Rvec = (fourPiBoverC - Erad0Vec / kappaPoverE) * tau0 + work; - // tau0 is used as a scaling factor for Rvec - if constexpr (use_D_as_base) { - for (int g = 0; g < nGroups_; ++g) { - if (tau0[g] <= 1.0) { - tau0[g] = 1.0; + tau0 = dt * rho * kappaPVec * chat * lorentz_factor; + tau = tau0; + Rvec = (fourPiBoverC - EradVec_guess / kappaPoverE) * tau0 + work; + if constexpr (use_D_as_base) { + // tau0 is used as a scaling factor for Rvec + for (int g = 0; g < nGroups_; ++g) { + if (tau0[g] <= 1.0) { + tau0[g] = 1.0; + } + } + D = Rvec / tau0; + } + } else { // in the second and later loops, calculate tau and E (given R) + tau = dt * rho * kappaPVec * chat * lorentz_factor; + if constexpr (use_D_as_base) { + Rvec = tau0 * D; + } + for (int g = 0; g < nGroups_; ++g) { + // If tau = 0.0, Erad_guess shouldn't change + if (tau[g] > 0.0) { + EradVec_guess[g] = kappaPoverE[g] * (fourPiBoverC[g] - (Rvec[g] - work[g]) / tau[g]); + if constexpr (force_rad_floor_in_iteration) { + if (EradVec_guess[g] < 0.0) { + Egas_guess -= (c_light_ / c_hat_) * (Erad_floor_ - EradVec_guess[g]); + EradVec_guess[g] = Erad_floor_; + } + } + } } } - D = Rvec / tau0; - } - - double F_G = NAN; - double dFG_dEgas = NAN; - double deltaEgas = NAN; - quokka::valarray dFG_dD{}; - quokka::valarray dFR_dEgas{}; - quokka::valarray dFR_i_dD_i{}; - quokka::valarray deltaD{}; - quokka::valarray F_D{}; - tau = tau0; - const double resid_tol = 1.0e-11; // 1.0e-15; - const int maxIter = 400; - int n = 0; - for (; n < maxIter; ++n) { - // F_G = Egas_guess - Egas0 + (c / chat) * sum(Rvec); F_G = Egas_guess - Egas0; F_D = EradVec_guess - Erad0Vec - (Rvec + Src); double F_D_abs_sum = 0.0; @@ -1527,63 +1557,6 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne Rvec += deltaD; } - // compute material temperature - T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); - AMREX_ASSERT(T_gas >= 0.); - // compute opacity, emissivity - fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); - - if constexpr (opacity_model_ == OpacityModel::single_group) { - kappaPVec[0] = ComputePlanckOpacity(rho, T_gas); - kappaEVec[0] = ComputeEnergyMeanOpacity(rho, T_gas); - } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { - kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); - for (int g = 0; g < nGroups_; ++g) { - kappaPVec[g] = kappa_expo_and_lower_value[1][g]; - kappaEVec[g] = kappa_expo_and_lower_value[1][g]; - } - } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { - kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); - // kappaPVec = ComputeGroupMeanOpacityWithMinusOneSlope(kappa_expo_and_lower_value, radBoundaryRatios_copy); - kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); - kappaEVec = kappaPVec; - } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { - kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); - if (n < max_ite_to_update_alpha_E) { - alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); - alpha_E = ComputeRadQuantityExponents(EradVec_guess, radBoundaries_g_copy); - } - kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); - kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); - } - AMREX_ASSERT(!kappaPVec.hasnan()); - AMREX_ASSERT(!kappaEVec.hasnan()); - - for (int g = 0; g < nGroups_; ++g) { - if (kappaEVec[g] > 0.0) { - kappaPoverE[g] = kappaPVec[g] / kappaEVec[g]; - } else { - kappaPoverE[g] = 1.0; - } - } - - tau = dt * rho * kappaEVec * chat * lorentz_factor; - if constexpr (use_D_as_base) { - Rvec = tau0 * D; - } - for (int g = 0; g < nGroups_; ++g) { - // If tau = 0.0, Erad_guess shouldn't change - if (tau[g] > 0.0) { - EradVec_guess[g] = kappaPoverE[g] * (fourPiBoverC[g] - (Rvec[g] - work[g]) / tau[g]); - if constexpr (force_rad_floor_in_iteration) { - if (EradVec_guess[g] < 0.0) { - Egas_guess -= (c_light_ / c_hat_) * (Erad_floor_ - EradVec_guess[g]); - EradVec_guess[g] = Erad_floor_; - } - } - } - } - // check relative and absolute convergence of E_r // if (std::abs(deltaEgas / Egas_guess) < 1e-7) { // break; @@ -1594,100 +1567,59 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // std::cout << "Newton-Raphson converged after " << n << " it." << std::endl; AMREX_ALWAYS_ASSERT(Egas_guess > 0.0); AMREX_ALWAYS_ASSERT(min(EradVec_guess) >= 0.0); - } // endif gamma != 1.0 - - // Erad_guess is the new radiation energy (excluding work term) - // Egas_guess is the new gas internal energy - - // 2. Compute radiation flux update - - amrex::GpuArray Frad_t0{}; - - T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); - if constexpr (gamma_ != 1.0) { - fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); - } - if constexpr (opacity_model_ == OpacityModel::single_group) { - kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_gas); // note that kappaFVec is used no matter what the value of gamma is - if constexpr (gamma_ != 1.0) { - kappaPVec[0] = ComputePlanckOpacity(rho, T_gas); - kappaEVec[0] = ComputeEnergyMeanOpacity(rho, T_gas); - AMREX_ASSERT(!std::isnan(kappaPVec[0])); - AMREX_ASSERT(!std::isnan(kappaEVec[0])); - } - } else { - kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); - if constexpr (gamma_ != 1.0) { - for (int g = 0; g < nGroups_; ++g) { - auto const nu_L = radBoundaries_g_copy[g]; - auto const nu_R = radBoundaries_g_copy[g + 1]; - auto const B_L = PlanckFunction(nu_L, T_gas); // 4 pi B(nu) / c - auto const B_R = PlanckFunction(nu_R, T_gas); // 4 pi B(nu) / c - if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { - delta_nu_kappa_B_at_edge[g] = kappa_expo_and_lower_value[1][g] * (nu_R * B_R - nu_L * B_L); - } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum || - opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + if (n > 0) { + // calculate kappaF since the temperature has changed + if constexpr (opacity_model_ == OpacityModel::single_group) { + kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_gas); + } else { + for (int g = 0; g < nGroups_; ++g) { + auto const nu_L = radBoundaries_g_copy[g]; + auto const nu_R = radBoundaries_g_copy[g + 1]; + auto const B_L = PlanckFunction(nu_L, T_gas); // 4 pi B(nu) / c + auto const B_R = PlanckFunction(nu_R, T_gas); // 4 pi B(nu) / c auto const kappa_L = kappa_expo_and_lower_value[1][g]; auto const kappa_R = kappa_L * std::pow(nu_R / nu_L, kappa_expo_and_lower_value[0][g]); delta_nu_kappa_B_at_edge[g] = nu_R * kappa_R * B_R - nu_L * kappa_L * B_L; delta_nu_B_at_edge[g] = nu_R * B_R - nu_L * B_L; } + if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + kappaFVec = kappaPVec; + } else { + if constexpr (use_diffuse_flux_mean_opacity) { + kappaFVec = ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, delta_nu_kappa_B_at_edge, + delta_nu_B_at_edge, kappa_expo_and_lower_value[0]); + } else { + // for simplicity, I assume kappaF = kappaE when opacity_model_ == OpacityModel::PPL_opacity_full_spectrum, if !use_diffuse_flux_mean_opacity. + // We won't use this option anyway. + kappaFVec = kappaEVec; + } + } } } - if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + } else { // if constexpr gamma_ == 1.0 + if constexpr (opacity_model_ == OpacityModel::single_group) { + kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_gas); + } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); for (int g = 0; g < nGroups_; ++g) { kappaFVec[g] = kappa_expo_and_lower_value[1][g]; - if constexpr (gamma_ != 1.0) { - kappaPVec[g] = kappaFVec[g]; - kappaEVec[g] = kappaFVec[g]; - } } } else { - if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { - if constexpr (gamma_ != 1.0) { - kappaPVec = - ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); - kappaEVec = kappaPVec; - if constexpr (use_diffuse_flux_mean_opacity) { - kappaFVec = ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, - delta_nu_kappa_B_at_edge, delta_nu_B_at_edge, - kappa_expo_and_lower_value[0]); - } else { - kappaFVec = kappaPVec; - } - } else { - kappaFVec = - ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); - } - } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { - // Note that alpha_F has not been changed in the Newton iteration - if constexpr (gamma_ != 1.0) { - alpha_B = ComputeRadQuantityExponents(fourPiBoverC, radBoundaries_g_copy); - alpha_E = ComputeRadQuantityExponents(EradVec_guess, radBoundaries_g_copy); - kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_B); - kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_E); - AMREX_ASSERT(!kappaPVec.hasnan()); - AMREX_ASSERT(!kappaEVec.hasnan()); - if constexpr (use_diffuse_flux_mean_opacity) { - kappaFVec = ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, - delta_nu_kappa_B_at_edge, delta_nu_B_at_edge, - kappa_expo_and_lower_value[0]); - } else { // fall back to bin-center opacity - kappaFVec = ComputeBinCenterOpacity(radBoundaries_g_copy, kappa_expo_and_lower_value); - } - } else { - kappaFVec = - ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); - } - } + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + kappaFVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); } } + // Erad_guess is the new radiation energy (excluding work term) + // Egas_guess is the new gas internal energy + + // 2. Compute radiation flux update + + amrex::GpuArray Frad_t0{}; dMomentum = {0., 0., 0.}; for (int g = 0; g < nGroups_; ++g) { - Frad_t0[0] = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); Frad_t0[1] = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); Frad_t0[2] = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); @@ -1792,7 +1724,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } } } - } else { + } else { // if constexpr (gamma_ == 1.0 || beta_order_ == 0) for (int n = 0; n < 3; ++n) { Frad_t1[n][g] = Frad_t0[n] / (1.0 + rho * kappaFVec[g] * chat * dt); // Compute conservative gas momentum update @@ -1854,34 +1786,34 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne if constexpr ((beta_order_ == 0) || (gamma_ == 1.0) || (!include_work_term_in_source)) { break; - } - - // If you are here, then you are using the new scheme. Step 3 is skipped. The work term is included in the source term, but it is - // lagged. The work term is updated in the next step. - for (int g = 0; g < nGroups_; ++g) { - // copy work to work_prev - work_prev[g] = work[g]; - // compute new work term from the updated radiation flux and velocity - // work = v * F * chi - if constexpr (opacity_model_ == OpacityModel::single_group) { - work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * chat / (c * c) * - lorentz_factor_v * (2.0 * kappaEVec[g] - kappaFVec[g]) * dt; - } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { - work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * kappaFVec[g] * chat / - (c * c) * dt; - } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum || - opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { - work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * - (1.0 + kappa_expo_and_lower_value[0][g]) * kappaFVec[g] * chat / (c * c) * dt; + } else { + // If you are here, then you are using the new scheme. Step 3 is skipped. The work term is included in the source term, but it is + // lagged. The work term is updated in the next step. + for (int g = 0; g < nGroups_; ++g) { + // copy work to work_prev + work_prev[g] = work[g]; + // compute new work term from the updated radiation flux and velocity + // work = v * F * chi + if constexpr (opacity_model_ == OpacityModel::single_group) { + work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * chat / (c * c) * + lorentz_factor_v * (2.0 * kappaEVec[g] - kappaFVec[g]) * dt; + } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * kappaFVec[g] * chat / + (c * c) * dt; + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum || + opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * + (1.0 + kappa_expo_and_lower_value[0][g]) * kappaFVec[g] * chat / (c * c) * dt; + } } - } - // Check for convergence of the work term: if the relative change in the work term is less than 1e-13, then break the loop - const double lag_tol = 1.0e-13; - if ((sum(abs(work)) == 0.0) || ((c / chat) * sum(abs(work - work_prev)) / Etot0 < lag_tol) || - (sum(abs(work - work_prev)) <= lag_tol * sum(Rvec)) || - (sum(abs(work)) > 0.0 && sum(abs(work - work_prev)) <= 1.0e-8 * sum(abs(work)))) { - break; + // Check for convergence of the work term: if the relative change in the work term is less than 1e-13, then break the loop + const double lag_tol = 1.0e-13; + if ((sum(abs(work)) == 0.0) || ((c / chat) * sum(abs(work - work_prev)) / Etot0 < lag_tol) || + (sum(abs(work - work_prev)) <= lag_tol * sum(Rvec)) || + (sum(abs(work)) > 0.0 && sum(abs(work - work_prev)) <= 1.0e-8 * sum(abs(work)))) { + break; + } } } // end full-step iteration From f0962c7faab6840037959b59de137b81a14174c3 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 15 Aug 2024 18:25:44 +1000 Subject: [PATCH 04/25] revert gitignore --- .gitignore | 3 --- 1 file changed, 3 deletions(-) diff --git a/.gitignore b/.gitignore index 9de410af7..bf8e1ceb9 100644 --- a/.gitignore +++ b/.gitignore @@ -17,9 +17,6 @@ sims/ __pycache__ *.code-workspace *.sarif -Backtrace.0 - -/docs/site/ *.aux *.dbl From 7dae30ffff62ccec514d592bd556547c416ec22f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 15 Aug 2024 08:29:13 +0000 Subject: [PATCH 05/25] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/radiation/radiation_system.hpp | 53 ++++++++++++++++-------------- 1 file changed, 29 insertions(+), 24 deletions(-) diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 4e1e63226..202c4213f 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -1382,7 +1382,8 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne kappaEVec[g] = kappa_expo_and_lower_value[1][g]; } } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { - kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); + kappaPVec = + ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); kappaEVec = kappaPVec; } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { if (n < max_ite_to_update_alpha_E) { @@ -1422,12 +1423,13 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne kappaFVec = kappaPVec; } else { if constexpr (use_diffuse_flux_mean_opacity) { - kappaFVec = - ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, delta_nu_kappa_B_at_edge, - delta_nu_B_at_edge, kappa_expo_and_lower_value[0]); + kappaFVec = ComputeDiffusionFluxMeanOpacity( + kappaPVec, kappaEVec, fourPiBoverC, delta_nu_kappa_B_at_edge, delta_nu_B_at_edge, + kappa_expo_and_lower_value[0]); } else { - // for simplicity, I assume kappaF = kappaE when opacity_model_ == OpacityModel::PPL_opacity_full_spectrum, if !use_diffuse_flux_mean_opacity. - // We won't use this option anyway. + // for simplicity, I assume kappaF = kappaE when opacity_model_ == + // OpacityModel::PPL_opacity_full_spectrum, if !use_diffuse_flux_mean_opacity. We won't + // use this option anyway. kappaFVec = kappaEVec; } } @@ -1444,25 +1446,26 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne const double frad2 = consPrev(i, j, k, x3RadFlux_index); // work = v * F * chi work[0] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * - (2.0 * kappaEVec[0] - kappaFVec[0]) * chat / (c * c) * lorentz_factor_v * dt; + (2.0 * kappaEVec[0] - kappaFVec[0]) * chat / (c * c) * lorentz_factor_v * dt; } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { for (int g = 0; g < nGroups_; ++g) { const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); // work = v * F * chi - work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * kappaFVec[g] * chat / - (c * c) * dt; + work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * + kappaFVec[g] * chat / (c * c) * dt; } } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum || - opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { for (int g = 0; g < nGroups_; ++g) { const double frad0 = consPrev(i, j, k, x1RadFlux_index + numRadVars_ * g); const double frad1 = consPrev(i, j, k, x2RadFlux_index + numRadVars_ * g); const double frad2 = consPrev(i, j, k, x3RadFlux_index + numRadVars_ * g); // work = v * F * chi work[g] = (x1GasMom0 * frad0 + x2GasMom0 * frad1 + x3GasMom0 * frad2) * - (1.0 + kappa_expo_and_lower_value[0][g]) * kappaFVec[g] * chat / (c * c) * dt; + (1.0 + kappa_expo_and_lower_value[0][g]) * kappaFVec[g] * chat / + (c * c) * dt; } } } @@ -1587,11 +1590,13 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne kappaFVec = kappaPVec; } else { if constexpr (use_diffuse_flux_mean_opacity) { - kappaFVec = ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, delta_nu_kappa_B_at_edge, - delta_nu_B_at_edge, kappa_expo_and_lower_value[0]); + kappaFVec = ComputeDiffusionFluxMeanOpacity(kappaPVec, kappaEVec, fourPiBoverC, + delta_nu_kappa_B_at_edge, delta_nu_B_at_edge, + kappa_expo_and_lower_value[0]); } else { - // for simplicity, I assume kappaF = kappaE when opacity_model_ == OpacityModel::PPL_opacity_full_spectrum, if !use_diffuse_flux_mean_opacity. - // We won't use this option anyway. + // for simplicity, I assume kappaF = kappaE when opacity_model_ == + // OpacityModel::PPL_opacity_full_spectrum, if !use_diffuse_flux_mean_opacity. We won't use this + // option anyway. kappaFVec = kappaEVec; } } @@ -1787,8 +1792,8 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne if constexpr ((beta_order_ == 0) || (gamma_ == 1.0) || (!include_work_term_in_source)) { break; } else { - // If you are here, then you are using the new scheme. Step 3 is skipped. The work term is included in the source term, but it is - // lagged. The work term is updated in the next step. + // If you are here, then you are using the new scheme. Step 3 is skipped. The work term is included in the source term, but it + // is lagged. The work term is updated in the next step. for (int g = 0; g < nGroups_; ++g) { // copy work to work_prev work_prev[g] = work[g]; @@ -1796,22 +1801,22 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // work = v * F * chi if constexpr (opacity_model_ == OpacityModel::single_group) { work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * chat / (c * c) * - lorentz_factor_v * (2.0 * kappaEVec[g] - kappaFVec[g]) * dt; + lorentz_factor_v * (2.0 * kappaEVec[g] - kappaFVec[g]) * dt; } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { - work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * kappaFVec[g] * chat / - (c * c) * dt; + work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * kappaFVec[g] * + chat / (c * c) * dt; } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum || - opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { work[g] = (x1GasMom1 * Frad_t1[0][g] + x2GasMom1 * Frad_t1[1][g] + x3GasMom1 * Frad_t1[2][g]) * - (1.0 + kappa_expo_and_lower_value[0][g]) * kappaFVec[g] * chat / (c * c) * dt; + (1.0 + kappa_expo_and_lower_value[0][g]) * kappaFVec[g] * chat / (c * c) * dt; } } // Check for convergence of the work term: if the relative change in the work term is less than 1e-13, then break the loop const double lag_tol = 1.0e-13; if ((sum(abs(work)) == 0.0) || ((c / chat) * sum(abs(work - work_prev)) / Etot0 < lag_tol) || - (sum(abs(work - work_prev)) <= lag_tol * sum(Rvec)) || - (sum(abs(work)) > 0.0 && sum(abs(work - work_prev)) <= 1.0e-8 * sum(abs(work)))) { + (sum(abs(work - work_prev)) <= lag_tol * sum(Rvec)) || + (sum(abs(work)) > 0.0 && sum(abs(work - work_prev)) <= 1.0e-8 * sum(abs(work)))) { break; } } From ccef32b435b0cd0d2f366b8a8fa2745e6c7a7061 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 15 Aug 2024 19:59:55 +1000 Subject: [PATCH 06/25] merge in RadDust folder from T_d branch --- src/problems/CMakeLists.txt | 1 + src/problems/RadDust/CMakeLists.txt | 7 + src/problems/RadDust/test_rad_dust.cpp | 273 +++++++++++++++++++++++++ src/problems/RadDust/test_rad_dust.hpp | 26 +++ tests/RadDust.in | 23 +++ 5 files changed, 330 insertions(+) create mode 100644 src/problems/RadDust/CMakeLists.txt create mode 100644 src/problems/RadDust/test_rad_dust.cpp create mode 100644 src/problems/RadDust/test_rad_dust.hpp create mode 100644 tests/RadDust.in diff --git a/src/problems/CMakeLists.txt b/src/problems/CMakeLists.txt index d710eca96..0ca223b15 100644 --- a/src/problems/CMakeLists.txt +++ b/src/problems/CMakeLists.txt @@ -33,6 +33,7 @@ add_subdirectory(RadStreamingY) add_subdirectory(RadSuOlson) add_subdirectory(RadTophat) add_subdirectory(RadTube) +add_subdirectory(RadDust) add_subdirectory(RadhydroShell) add_subdirectory(RadhydroShock) diff --git a/src/problems/RadDust/CMakeLists.txt b/src/problems/RadDust/CMakeLists.txt new file mode 100644 index 000000000..15b8ac221 --- /dev/null +++ b/src/problems/RadDust/CMakeLists.txt @@ -0,0 +1,7 @@ +add_executable(test_rad_dust test_rad_dust.cpp ../../util/fextract.cpp ../../math/interpolate.cpp ${QuokkaObjSources}) + +if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_rad_dust) +endif(AMReX_GPU_BACKEND MATCHES "CUDA") + +add_test(NAME RadhydroUniformAdvecting COMMAND test_rad_dust RadDust.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/problems/RadDust/test_rad_dust.cpp b/src/problems/RadDust/test_rad_dust.cpp new file mode 100644 index 000000000..c3dad98c5 --- /dev/null +++ b/src/problems/RadDust/test_rad_dust.cpp @@ -0,0 +1,273 @@ +/// \file test_rad_dust.cpp +/// \brief Defines a test problem for radiation advection in a uniform medium with grey radiation. +/// + +#include "test_rad_dust.hpp" +#include "AMReX_BC_TYPES.H" +#include "AMReX_Print.H" +#include "QuokkaSimulation.hpp" +#include "fundamental_constants.H" +#include "physics_info.hpp" +#include "util/fextract.hpp" +#include + +struct DustProblem { +}; // dummy type to allow compile-type polymorphism via template specialization + +// CGS + +constexpr int beta_order_ = 1; // order of beta in the radiation four-force +constexpr double c = 1.0e8; +constexpr double chat = c; +constexpr double v0 = 0.0; +constexpr double chi0 = 10000.0; + +constexpr double T0 = 1.0; +constexpr double T_equi = 0.7680325; +constexpr double rho0 = 1.0; +constexpr double a_rad = 1.0; +constexpr double mu = 1.0; +constexpr double k_B = 1.0; + +constexpr double max_time = 1.0e-5; +constexpr double delta_time = 1.0e-8; + +constexpr double Erad0 = a_rad * T0 * T0 * T0 * T0; +constexpr double erad_floor = 1.0e-20 * Erad0; + +template <> struct SimulationData { + std::vector t_vec_; + std::vector Trad_vec_; + std::vector Tgas_vec_; +}; + +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = mu; + static constexpr double boltzmann_constant = k_B; + static constexpr double gamma = 5. / 3.; +}; + +template <> struct RadSystem_Traits { + 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 = beta_order_; +}; + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + 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; +}; + +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double rho, const double /*Tgas*/) -> amrex::Real +{ + return chi0 / rho; +} + +template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanOpacity(const double rho, const double Tgas) -> amrex::Real +{ + return ComputePlanckOpacity(rho, Tgas); +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiation(amrex::Real temperature, amrex::GpuArray const &boundaries) + -> quokka::valarray +{ + auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); + const double power = radiation_constant_ * temperature; + auto Erad_g = power * radEnergyFractions; + return Erad_g; +} + +template <> +AMREX_GPU_HOST_DEVICE auto +RadSystem::ComputeThermalRadiationTempDerivative(amrex::Real temperature, + amrex::GpuArray const &boundaries) -> quokka::valarray +{ + auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); + const double d_power_dt = radiation_constant_; + return d_power_dt * radEnergyFractions; +} + +template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) +{ + // extract variables required from the geom object + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + const double Egas = quokka::EOS::ComputeEintFromTgas(rho0, T0); + + double erad = NAN; + double frad = NAN; + erad = Erad0; + frad = 0.0; + + // test + erad = erad_floor; + frad = 0.0; + + // loop over the grid and set the initial condition + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + state_cc(i, j, k, RadSystem::radEnergy_index) = erad; + state_cc(i, j, k, RadSystem::x1RadFlux_index) = frad; + state_cc(i, j, k, RadSystem::x2RadFlux_index) = 0; + state_cc(i, j, k, RadSystem::x3RadFlux_index) = 0; + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho0 * v0 * v0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho0; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = v0 * rho0; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} + +template <> void QuokkaSimulation::computeAfterTimestep() +{ + auto [position, values] = fextract(state_new_cc_[0], Geom(0), 0, 0.5); + + if (amrex::ParallelDescriptor::IOProcessor()) { + 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 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]; + const amrex::Real Egas_i = RadSystem::ComputeEintFromEgas(rho, x1GasMom, x2GasMom, x3GasMom, Etot_i); + const amrex::Real Erad_i = values.at(RadSystem::radEnergy_index)[0]; + // userData_.Trad_vec_.push_back(std::pow(Erad_i / a_rad, 1. / 4.)); + userData_.Trad_vec_.push_back(Erad_i / a_rad); + userData_.Tgas_vec_.push_back(quokka::EOS::ComputeTgasFromEint(rho, Egas_i)); + } +} + +auto problem_main() -> int +{ + // Problem parameters + const int max_timesteps = 1e6; + const double CFL_number_gas = 0.8; + const double CFL_number_rad = 8.0; + + // Boundary conditions + constexpr int nvars = RadSystem::nvar_; + amrex::Vector BCs_cc(nvars); + for (int n = 0; n < nvars; ++n) { + for (int i = 0; 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 sim(BCs_cc); + + sim.radiationReconstructionOrder_ = 3; // PPM + sim.stopTime_ = max_time; + sim.radiationCflNumber_ = CFL_number_rad; + sim.cflNumber_ = CFL_number_gas; + sim.maxTimesteps_ = max_timesteps; + sim.plotfileInterval_ = -1; + sim.initDt_ = delta_time; + sim.maxDt_ = delta_time; + + // initialize + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + // read in exact solution + std::vector ts_exact{}; + std::vector Trad_exact{}; + std::vector Tgas_exact{}; + + std::ifstream fstream("../extern/data/dust/rad_dust_exact.csv", std::ios::in); + AMREX_ALWAYS_ASSERT(fstream.is_open()); + std::string header; + std::getline(fstream, header); + + for (std::string line; std::getline(fstream, line);) { + std::istringstream iss(line); + std::vector values; + std::string value; + + while (std::getline(iss, value, ',')) { + values.push_back(std::stod(value)); + } + auto t_val = values.at(0); + auto Tmat_val = values.at(1); + auto Trad_val = values.at(2); + if (t_val <= 0.0) { + continue; + } + ts_exact.push_back(t_val); + Tgas_exact.push_back(Tmat_val); + Trad_exact.push_back(Trad_val); + } + + std::vector &Tgas = sim.userData_.Tgas_vec_; + std::vector &Trad = sim.userData_.Trad_vec_; + std::vector &t = sim.userData_.t_vec_; + + std::vector Tgas_interp(t.size()); + std::vector Trad_interp(t.size()); + interpolate_arrays(t.data(), Tgas_interp.data(), static_cast(t.size()), ts_exact.data(), Tgas_exact.data(), static_cast(ts_exact.size())); + interpolate_arrays(t.data(), Trad_interp.data(), static_cast(t.size()), ts_exact.data(), Trad_exact.data(), static_cast(ts_exact.size())); + + // compute error norm + double err_norm = 0.; + double sol_norm = 0.; + for (size_t i = 0; i < t.size(); ++i) { + err_norm += std::abs(Tgas[i] - Tgas_interp[i]); + err_norm += std::abs(Trad[i] - Trad_interp[i]); + sol_norm += std::abs(Tgas_interp[i]) + std::abs(Trad_interp[i]); + } + const double error_tol = 0.0008; + const double rel_error = err_norm / sol_norm; + amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; + +#ifdef HAVE_PYTHON + // plot temperature + matplotlibcpp::clf(); + matplotlibcpp::xscale("log"); + std::map Trad_args; + std::map Tgas_args; + std::map Texact_args; + std::map Tradexact_args; + Trad_args["label"] = "radiation (numerical)"; + Trad_args["linestyle"] = "--"; + Trad_args["color"] = "C1"; + Tradexact_args["label"] = "radiation (exact)"; + Tradexact_args["linestyle"] = "-"; + Tradexact_args["color"] = "k"; + Tgas_args["label"] = "gas (numerical)"; + Tgas_args["linestyle"] = "--"; + Tgas_args["color"] = "C2"; + Texact_args["label"] = "gas (exact)"; + Texact_args["linestyle"] = "-"; + Texact_args["color"] = "k"; + matplotlibcpp::plot(ts_exact, Tgas_exact, Texact_args); + matplotlibcpp::plot(ts_exact, Trad_exact, Tradexact_args); + matplotlibcpp::plot(t, Tgas, Tgas_args); + matplotlibcpp::plot(t, Trad, Trad_args); + matplotlibcpp::xlabel("t (dimensionless)"); + matplotlibcpp::ylabel("T (dimensionless)"); + matplotlibcpp::legend(); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./rad_dust_T.pdf"); +#endif + + // Cleanup and exit + int status = 0; + if ((rel_error > error_tol) || std::isnan(rel_error)) { + status = 1; + } + return status; +} diff --git a/src/problems/RadDust/test_rad_dust.hpp b/src/problems/RadDust/test_rad_dust.hpp new file mode 100644 index 000000000..20ddcfd67 --- /dev/null +++ b/src/problems/RadDust/test_rad_dust.hpp @@ -0,0 +1,26 @@ +#ifndef TEST_RAD_DUST_HPP_ // NOLINT +#define TEST_RAD_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_rad_dust.hpp +/// \brief Defines a test problem for radiation in the static diffusion regime. +/// + +// external headers +#ifdef HAVE_PYTHON +#include "util/matplotlibcpp.h" +#endif +#include +#include + +// internal headers + +#include "math/interpolate.hpp" +#include "radiation/radiation_system.hpp" + +// function definitions + +#endif // TEST_RAD_DUST_HPP_ diff --git a/tests/RadDust.in b/tests/RadDust.in new file mode 100644 index 000000000..82d70ab3f --- /dev/null +++ b/tests/RadDust.in @@ -0,0 +1,23 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 1.0 1.0 1.0 +geometry.is_periodic = 1 1 1 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 0 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 8 4 4 +amr.max_level = 0 # number of levels = max_level + 1 +amr.blocking_factor = 4 # grid size must be divisible by this +amr.max_grid_size = 16 + +do_reflux = 0 +do_subcycle = 0 +suppress_output = 0 \ No newline at end of file From f4ae7ea7b8743df4a7b890255c4b00a2b1d4b718 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 15 Aug 2024 20:52:47 +1000 Subject: [PATCH 07/25] add exact solution to the dust problem --- extern/data/dust/rad_dust_exact.csv | 1026 +++++++++++++++++++++++++++ 1 file changed, 1026 insertions(+) create mode 100644 extern/data/dust/rad_dust_exact.csv diff --git a/extern/data/dust/rad_dust_exact.csv b/extern/data/dust/rad_dust_exact.csv new file mode 100644 index 000000000..99d48b6f9 --- /dev/null +++ b/extern/data/dust/rad_dust_exact.csv @@ -0,0 +1,1026 @@ +"t","T_gas","T_rad" +0.,1.,0. +9.765625e-9,0.9935526673232676,0.00967099901509845 +1.953125e-8,0.9872294252476815,0.019155862128477617 +2.9296875e-8,0.9810273280234489,0.028459007964826463 +3.90625e-8,0.9749434809259232,0.03758477861111502 +4.8828125e-8,0.9689751052086036,0.04653734218709428 +5.859375e-8,0.9631194765946542,0.05532078510801854 +6.8359375e-8,0.9573739704351772,0.063939044347234 +7.8125e-8,0.9517360182313913,0.07239597265291281 +8.7890625e-8,0.9462031444578127,0.08069528331328064 +9.765625e-8,0.9407729262647306,0.08884061060290394 +1.07421875e-7,0.9354430252098653,0.09683546218520178 +1.171875e-7,0.9302111467014769,0.10468327994778455 +1.26953125e-7,0.9250751030953223,0.11238734535701654 +1.3671875e-7,0.9200327252100756,0.1199509121848865 +1.46484375e-7,0.9150819074757611,0.12737713878635823 +1.5625e-7,0.9102206397248788,0.13466904041268174 +1.66015625e-7,0.9054469210448486,0.141829618432727 +1.7578125e-7,0.9007588163700884,0.14886177544486734 +1.85546875e-7,0.8961544778423224,0.15576828323651626 +1.953125e-7,0.89163206471411,0.16255190292883487 +2.05078125e-7,0.8871897942034853,0.16921530869477192 +2.1484375e-7,0.8828259614535108,0.1757610578197337 +2.24609375e-7,0.878538868494201,0.18219169725869833 +2.34375e-7,0.874326871792332,0.18850969231150191 +2.44140625e-7,0.8701883968831077,0.1947174046753384 +2.5390625e-7,0.8661218752973174,0.2008171870540237 +2.63671875e-7,0.8621257899451216,0.2068113150823175 +2.734375e-7,0.8581986849297515,0.2127019726053727 +2.83203125e-7,0.8543391075105059,0.21849133873424118 +2.9296875e-7,0.8505456461009504,0.22418153084857445 +3.02734375e-7,0.8468169693360215,0.22977454599596786 +3.125e-7,0.8431517534676893,0.23527236979846614 +3.22265625e-7,0.8395486749375626,0.24067698759365622 +3.3203125e-7,0.8360064492033011,0.2459903261950484 +3.41796875e-7,0.8325238769703621,0.251214184544457 +3.515625e-7,0.8290997700971222,0.2563503448543168 +3.61328125e-7,0.8257329404419582,0.2614005893370627 +3.7109375e-7,0.8224222160776555,0.26636667588351676 +3.80859375e-7,0.8191665033957773,0.27125024490633404 +3.90625e-7,0.8159647323253035,0.27605290151204487 +4.00390625e-7,0.8128158327980433,0.28077625080293517 +4.1015625e-7,0.8097187403937732,0.2854218894093403 +4.19921875e-7,0.8066724517098149,0.28999132243527775 +4.296875e-7,0.8036760004310484,0.29448599935342745 +4.39453125e-7,0.8007284207532316,0.2989073688701526 +4.4921875e-7,0.7978287481208854,0.30325687781867183 +4.58984375e-7,0.7949760584178817,0.30753591237317757 +4.6875e-7,0.7921694754700492,0.3117457867949262 +4.78515625e-7,0.7894081258536431,0.3158878112195355 +4.8828125e-7,0.7866911362156422,0.3199632956765367 +4.98046875e-7,0.784017655188405,0.3239735172173926 +5.078125e-7,0.7813868837103285,0.32791967443450737 +5.17578125e-7,0.7787980302707683,0.33180295459384773 +5.2734375e-7,0.77625030335908,0.3356245449613803 +5.37109375e-7,0.773742920909645,0.3393856186355326 +5.46875e-7,0.771275149033554,0.3430872764496693 +5.56640625e-7,0.768846269054321,0.3467305964185188 +5.6640625e-7,0.766455562298091,0.35031665655286387 +5.76171875e-7,0.7641023092359293,0.35384653614610634 +5.859375e-7,0.7617858227221189,0.35732126591682195 +5.95703125e-7,0.759505459989334,0.36074181001599914 +6.0546875e-7,0.7572605803043606,0.36410912954345936 +6.15234375e-7,0.755050542933984,0.3674241855990243 +6.25e-7,0.7528747071449898,0.37068793928251553 +6.34765625e-7,0.7507324348251997,0.37390134776220074 +6.4453125e-7,0.748623130076254,0.3770653048856191 +6.54296875e-7,0.7465462304979351,0.38018065425309755 +6.640625e-7,0.7445011745585666,0.3832482381621505 +6.73828125e-7,0.7424874007264719,0.38626889891029237 +6.8359375e-7,0.7405043474699751,0.38924347879503773 +6.93359375e-7,0.7385514545195476,0.3921728182206789 +7.03125e-7,0.7366281935472386,0.3950577096791424 +7.12890625e-7,0.7347340692679728,0.39789889609804113 +7.2265625e-7,0.7328685878640794,0.40069711820388126 +7.32421875e-7,0.7310312555178877,0.4034531167231688 +7.421875e-7,0.729221578411727,0.4061676323824098 +7.51953125e-7,0.7274390632589113,0.40884140511163336 +7.6171875e-7,0.7256832403599346,0.41147513946009845 +7.71484375e-7,0.7239536718387849,0.41406949224182293 +7.8125e-7,0.7222499220049802,0.4166251169925301 +7.91015625e-7,0.7205715551680382,0.41914266724794297 +8.0078125e-7,0.718918135637477,0.4216227965437849 +8.10546875e-7,0.7172892279011063,0.42406615814834087 +8.203125e-7,0.7156844134050308,0.4264733798924541 +8.30078125e-7,0.7141033035705255,0.4288450446442119 +8.3984375e-7,0.7125455128094647,0.4311817307858033 +8.49609375e-7,0.7110106555337221,0.4334840166994171 +8.59375e-7,0.709498346155172,0.43575248076724216 +8.69140625e-7,0.7080081991231668,0.43798770131525 +8.7890625e-7,0.7065398407058581,0.4401902389412131 +8.88671875e-7,0.7050929248114342,0.4423606127828489 +8.984375e-7,0.7036671092078198,0.4444993361882704 +9.08203125e-7,0.7022620516629395,0.44660692250559086 +9.1796875e-7,0.7008774099447178,0.4486838850829234 +9.27734375e-7,0.6995128418224691,0.4507307372662966 +9.375e-7,0.6981680130172093,0.4527479804741862 +9.47265625e-7,0.6968426141800251,0.4547360787299625 +9.5703125e-7,0.6955363407309759,0.45669548890353634 +9.66796875e-7,0.6942488880901212,0.4586266678648184 +9.765625e-7,0.6929799516775206,0.46053007248371935 +9.86328125e-7,0.6917292269132335,0.46240615963015 +9.9609375e-7,0.6904964113282019,0.4642553830076972 +1.005859375e-6,0.6892812257439994,0.46607816138400104 +1.015625e-6,0.6880834019334213,0.4678748970998682 +1.025390625e-6,0.6869026716847273,0.46964599247290917 +1.03515625e-6,0.685738766786177,0.4713918498207347 +1.044921875e-6,0.6845914190260298,0.4731128714609554 +1.0546875e-6,0.6834603601925455,0.47480945971118177 +1.064453125e-6,0.682345322645923,0.4764820160311157 +1.07421875e-6,0.6812460547319646,0.47813091790205303 +1.083984375e-6,0.6801623218049705,0.47975651729254426 +1.09375e-6,0.6790938899659903,0.48135916505101456 +1.103515625e-6,0.6780405253160741,0.482939212025889 +1.11328125e-6,0.6770019939562716,0.48449700906559273 +1.123046875e-6,0.6759780619876328,0.48603290701855084 +1.1328125e-6,0.6749684955112079,0.48754725673318833 +1.142578125e-6,0.6739730646471207,0.48904040302931895 +1.15234375e-6,0.6729915591850558,0.4905126612224164 +1.162109375e-6,0.6720237745607723,0.49196433815884155 +1.171875e-6,0.6710695061962452,0.49339574070563236 +1.181640625e-6,0.6701285495134487,0.494807175729827 +1.19140625e-6,0.6692006999343576,0.49619895009846365 +1.201171875e-6,0.6682857528809466,0.4975713706785802 +1.2109375e-6,0.6673835040078937,0.49892474398815956 +1.220703125e-6,0.6664937593543191,0.5002593609685215 +1.23046875e-6,0.6656163385912978,0.5015754921130535 +1.240234375e-6,0.6647510622382601,0.5028734066426099 +1.25e-6,0.6638977508146366,0.5041533737780451 +1.259765625e-6,0.6630562248398578,0.5054156627402133 +1.26953125e-6,0.6622263048333541,0.5066605427499689 +1.279296875e-6,0.661407811314556,0.507888283028166 +1.2890625e-6,0.6606005669576923,0.5090991495634615 +1.298828125e-6,0.6598044084500235,0.5102933873249648 +1.30859375e-6,0.6590191778588054,0.511471233211792 +1.318359375e-6,0.6582447172495286,0.512632924125707 +1.328125e-6,0.6574808686876841,0.5137786969684739 +1.337890625e-6,0.6567274742387625,0.5149087886418562 +1.34765625e-6,0.6559843759682548,0.5160234360476179 +1.357421875e-6,0.6552514160164322,0.5171228759753517 +1.3671875e-6,0.6545284430485814,0.5182073354271279 +1.376953125e-6,0.6538153167155442,0.5192770249266837 +1.38671875e-6,0.6531118976667172,0.5203321534999242 +1.396484375e-6,0.6524180465514969,0.5213729301727547 +1.40625e-6,0.65173362401928,0.52239956397108 +1.416015625e-6,0.651058490719463,0.5234122639208055 +1.42578125e-6,0.6503925073014425,0.5244112390478363 +1.435546875e-6,0.6497355355706362,0.5253966966440458 +1.4453125e-6,0.6490874473026599,0.5263688290460101 +1.455078125e-6,0.6484481192082455,0.5273278211876316 +1.46484375e-6,0.6478174280219424,0.5282738579670863 +1.474609375e-6,0.6471952504783,0.52920712428255 +1.484375e-6,0.6465814633118674,0.5301278050321988 +1.494140625e-6,0.6459759432571942,0.5310360851142086 +1.50390625e-6,0.645378567067376,0.5319321493989358 +1.513671875e-6,0.644789215606922,0.5328161765896168 +1.5234375e-6,0.6442077785131411,0.533688332230288 +1.533203125e-6,0.6436341465009007,0.5345487802486488 +1.54296875e-6,0.6430682102850677,0.5353976845723983 +1.552734375e-6,0.6425098605805094,0.5362352091292357 +1.5625e-6,0.6419589881020931,0.5370615178468602 +1.572265625e-6,0.6414154835646859,0.5378767746529709 +1.58203125e-6,0.6408792369162177,0.5386811446256732 +1.591796875e-6,0.6403501421601343,0.5394747867597983 +1.6015625e-6,0.639828104086679,0.5402578438699812 +1.611328125e-6,0.6393130281871416,0.5410304577192873 +1.62109375e-6,0.638804819952812,0.5417927700707817 +1.630859375e-6,0.63830338487498,0.5425449226875299 +1.640625e-6,0.6378086284449351,0.5432870573325971 +1.650390625e-6,0.6373204561539675,0.5440193157690485 +1.66015625e-6,0.6368387734933668,0.5447418397599496 +1.669921875e-6,0.6363634859544228,0.5454547710683656 +1.6796875e-6,0.6358944990284253,0.5461582514573619 +1.689453125e-6,0.6354317188259678,0.546852421761048 +1.69921875e-6,0.6349750600195935,0.5475374099706095 +1.708984375e-6,0.634524442847803,0.5482133357282953 +1.71875e-6,0.6340797875688184,0.5488803186467721 +1.728515625e-6,0.6336410144408621,0.5495384783387066 +1.73828125e-6,0.6332080437221561,0.5501879344167656 +1.748046875e-6,0.6327807956709226,0.5508288064936159 +1.7578125e-6,0.6323591905453837,0.5514612141819242 +1.767578125e-6,0.6319431486037617,0.5520852770943572 +1.77734375e-6,0.6315325901042786,0.5527011148435818 +1.787109375e-6,0.6311274353051568,0.5533088470422646 +1.796875e-6,0.6307276047623348,0.5539085928564976 +1.806640625e-6,0.6303330253881951,0.554500461917707 +1.81640625e-6,0.6299436295535721,0.5550845556696415 +1.826171875e-6,0.6295593497299579,0.5556609754050631 +1.8359375e-6,0.6291801183888445,0.5562298224167331 +1.845703125e-6,0.6288058680017241,0.5567911979974137 +1.85546875e-6,0.628436531040089,0.5573452034398664 +1.865234375e-6,0.6280720399754314,0.5578919400368528 +1.875e-6,0.6277123272792432,0.5584315090811349 +1.884765625e-6,0.6273573254230169,0.5589640118654744 +1.89453125e-6,0.6270069668782445,0.559489549682633 +1.904296875e-6,0.6266611842425095,0.5600082236362357 +1.9140625e-6,0.6263199147253893,0.5605201279119159 +1.923828125e-6,0.6259831007521037,0.5610253488718442 +1.93359375e-6,0.6256506849500094,0.5615239725749857 +1.943359375e-6,0.6253226099464632,0.5620160850803051 +1.953125e-6,0.6249988183688218,0.5625017724467672 +1.962890625e-6,0.6246792528444418,0.5629811207333371 +1.97265625e-6,0.6243638560006803,0.5634542159989794 +1.982421875e-6,0.6240525704648938,0.563921144302659 +1.9921875e-6,0.6237453388644392,0.564381991703341 +2.001953125e-6,0.6234421038266731,0.5648368442599901 +2.01171875e-6,0.6231428080220461,0.5652857879669306 +2.021484375e-6,0.6228473973840281,0.5657289039239577 +2.03125e-6,0.6225558227133627,0.5661662659299557 +2.041015625e-6,0.6222680351273524,0.5665979473089712 +2.05078125e-6,0.6219839857432993,0.5670240213850509 +2.060546875e-6,0.6217036256785055,0.5674445614822416 +2.0703125e-6,0.6214269060502733,0.56785964092459 +2.080078125e-6,0.6211537779759048,0.5682693330361427 +2.08984375e-6,0.6208841925727023,0.5686737111409464 +2.099609375e-6,0.620618100957968,0.569072848563048 +2.109375e-6,0.620355454249004,0.5694668186264938 +2.119140625e-6,0.6200962035724198,0.5698556946413702 +2.12890625e-6,0.6198403022964484,0.5702395465553274 +2.138671875e-6,0.6195877082286402,0.5706184376570396 +2.1484375e-6,0.6193383796129764,0.5709924305805352 +2.158203125e-6,0.619092274693438,0.5713615879598429 +2.16796875e-6,0.6188493517140057,0.5717259724289913 +2.177734375e-6,0.6186095689186607,0.5720856466220088 +2.1875e-6,0.6183728845513838,0.572440673172924 +2.197265625e-6,0.6181392568561562,0.5727911147157655 +2.20703125e-6,0.6179086440769587,0.5731370338845618 +2.216796875e-6,0.6176810044577722,0.5734784933133414 +2.2265625e-6,0.617456296242848,0.5738155556357277 +2.236328125e-6,0.6172344791668299,0.5741482812497549 +2.24609375e-6,0.6170155169274565,0.5744767246088149 +2.255859375e-6,0.6167993737777587,0.5748009393333615 +2.265625e-6,0.6165860139707675,0.5751209790438484 +2.275390625e-6,0.6163754017595138,0.5754368973607291 +2.28515625e-6,0.6161675013970282,0.5757487479044573 +2.294921875e-6,0.615962277136342,0.5760565842954867 +2.3046875e-6,0.6157596932304858,0.576360460154271 +2.314453125e-6,0.6155597139324906,0.5766604291012637 +2.32421875e-6,0.6153623034953873,0.5769565447569186 +2.333984375e-6,0.615167426172207,0.5772488607416892 +2.34375e-6,0.6149750471722468,0.5775374292416295 +2.353515625e-6,0.6147851351565901,0.5778222972651146 +2.36328125e-6,0.6145976594640024,0.5781035108039961 +2.373046875e-6,0.6144125894319113,0.5783811158521326 +2.3828125e-6,0.6142298943977446,0.5786551584033828 +2.392578125e-6,0.6140495436989297,0.5789256844516051 +2.40234375e-6,0.6138715066728943,0.5791927399906582 +2.412109375e-6,0.6136957526570658,0.579456371014401 +2.421875e-6,0.613522250988872,0.5797166235166916 +2.431640625e-6,0.6133509710057404,0.579973543491389 +2.44140625e-6,0.6131818820450987,0.5802271769323515 +2.451171875e-6,0.613014953092231,0.5804775703616532 +2.4609375e-6,0.6128501556254347,0.5807247665618477 +2.470703125e-6,0.612687463939222,0.5809688040911665 +2.48046875e-6,0.6125268523661768,0.5812097214507345 +2.490234375e-6,0.612368295238882,0.5814475571416766 +2.5e-6,0.6122117668899212,0.5816823496651178 +2.509765625e-6,0.6120572416518777,0.581914137522183 +2.51953125e-6,0.6119046938573348,0.5821429592139973 +2.529296875e-6,0.611754097838876,0.5823688532416855 +2.5390625e-6,0.6116054279290846,0.5825918581063727 +2.548828125e-6,0.6114586584605439,0.5828120123091838 +2.55859375e-6,0.6113137637658373,0.5830293543512437 +2.568359375e-6,0.6111707181775482,0.5832439227336773 +2.578125e-6,0.6110294960282598,0.5834557559576098 +2.587890625e-6,0.6108900716505558,0.583664892524166 +2.59765625e-6,0.6107524200909932,0.58387136986351 +2.607421875e-6,0.6106165197451393,0.5840752203822906 +2.6171875e-6,0.610482349682067,0.5842764754768991 +2.626953125e-6,0.610349888952715,0.5844751665709271 +2.63671875e-6,0.6102191166080219,0.5846713250879669 +2.646484375e-6,0.6100900116989261,0.5848649824516106 +2.65625e-6,0.6099625532763664,0.5850561700854501 +2.666015625e-6,0.6098367203912813,0.5852449194130777 +2.67578125e-6,0.6097124920946095,0.5854312618580854 +2.685546875e-6,0.6095898474372896,0.5856152288440652 +2.6953125e-6,0.6094687654702602,0.5857968517946095 +2.705078125e-6,0.6093492252444598,0.5859761621333099 +2.71484375e-6,0.6092312058108271,0.5861531912837589 +2.724609375e-6,0.6091146862203008,0.5863279706695483 +2.734375e-6,0.6089996455238195,0.5865005317142704 +2.744140625e-6,0.6088860632549539,0.5866709051175688 +2.75390625e-6,0.6087739216237654,0.5868391175643515 +2.763671875e-6,0.6086632035116374,0.5870051947325435 +2.7734375e-6,0.608553891780983,0.5871691623285252 +2.783203125e-6,0.6084459692942152,0.5873310460586769 +2.79296875e-6,0.6083394189137474,0.5874908716293785 +2.802734375e-6,0.6082342235019929,0.5876486647470104 +2.8125e-6,0.6081303659213647,0.5878044511179527 +2.822265625e-6,0.6080278290342761,0.5879582564485856 +2.83203125e-6,0.6079265957031403,0.5881101064452893 +2.841796875e-6,0.6078266487903705,0.588260026814444 +2.8515625e-6,0.60772797115838,0.5884080432624299 +2.861328125e-6,0.6076305456695817,0.5885541814956271 +2.87109375e-6,0.6075343551863892,0.588698467220416 +2.880859375e-6,0.6074393825712154,0.5888409261431766 +2.890625e-6,0.6073456110113941,0.5889815834829086 +2.900390625e-6,0.6072530258273621,0.5891204612589567 +2.91015625e-6,0.60716161299117,0.5892575805132447 +2.919921875e-6,0.6070713584565482,0.5893929623151773 +2.9296875e-6,0.6069822481772272,0.5895266277341589 +2.939453125e-6,0.6068942681069374,0.5896585978395936 +2.94921875e-6,0.6068074041994089,0.5897888937008863 +2.958984375e-6,0.6067216424083723,0.5899175363874412 +2.96875e-6,0.6066369686875579,0.5900445469686628 +2.978515625e-6,0.6065533689906961,0.5901699465139555 +2.98828125e-6,0.6064708292715173,0.5902937560927237 +2.998046875e-6,0.6063893354837518,0.590415996774372 +3.0078125e-6,0.6063088735811301,0.5905366896283046 +3.017578125e-6,0.6062294295173823,0.5906558557239262 +3.02734375e-6,0.606150989246239,0.5907735161306411 +3.037109375e-6,0.6060735389375688,0.5908896915936465 +3.046875e-6,0.6059970664537311,0.591004400319403 +3.056640625e-6,0.6059215602815826,0.5911176595776257 +3.06640625e-6,0.6058470088918474,0.5912294866622286 +3.076171875e-6,0.6057734007552495,0.5913398988671253 +3.0859375e-6,0.6057007243425132,0.5914489134862297 +3.095703125e-6,0.6056289681243625,0.5915565478134558 +3.10546875e-6,0.6055581205715217,0.5916628191427171 +3.115234375e-6,0.6054881701547146,0.5917677447679276 +3.125e-6,0.6054191053446656,0.5918713419830012 +3.134765625e-6,0.6053509146120987,0.5919736280818516 +3.14453125e-6,0.6052835864277379,0.5920746203583926 +3.154296875e-6,0.6052171092623075,0.5921743361065381 +3.1640625e-6,0.6051514715865317,0.592272792620202 +3.173828125e-6,0.6050866618711344,0.592370007193298 +3.18359375e-6,0.605022668729118,0.5924659969063225 +3.193359375e-6,0.604959482109843,0.5925607768352351 +3.203125e-6,0.6048970925492276,0.5926543611761581 +3.212890625e-6,0.604835490570007,0.5927467641449891 +3.22265625e-6,0.604774666694916,0.5928379999576256 +3.232421875e-6,0.6047146114466895,0.5929280828299652 +3.2421875e-6,0.6046553153480625,0.5930170269779057 +3.251953125e-6,0.6045967689217701,0.5931048466173444 +3.26171875e-6,0.6045389626905469,0.5931915559641792 +3.271484375e-6,0.604481887177128,0.5932771692343074 +3.28125e-6,0.6044255329042485,0.5933617006436267 +3.291015625e-6,0.6043698903946432,0.5934451644080347 +3.30078125e-6,0.6043149501710469,0.5935275747434292 +3.310546875e-6,0.6042607027561947,0.5936089458657074 +3.3203125e-6,0.6042071386728215,0.5936892919907673 +3.330078125e-6,0.6041542485367823,0.593768627194826 +3.33984375e-6,0.6041020240186976,0.593846963971953 +3.349609375e-6,0.6040504573308034,0.5939243140037944 +3.359375e-6,0.6039995406751419,0.5940006889872865 +3.369140625e-6,0.6039492662537557,0.5940761006193661 +3.37890625e-6,0.6038996262686867,0.5941505605969695 +3.388671875e-6,0.6038506129219775,0.5942240806170332 +3.3984375e-6,0.6038022184156704,0.594296672376494 +3.408203125e-6,0.6037544349518076,0.5943683475722882 +3.41796875e-6,0.6037072547324315,0.5944391179013524 +3.427734375e-6,0.6036606699595843,0.594508995060623 +3.4375e-6,0.6036146728353086,0.5945779907470368 +3.447265625e-6,0.6035692555616463,0.59464611665753 +3.45703125e-6,0.6035244103406402,0.5947133844890393 +3.466796875e-6,0.6034801293743322,0.5947798059385012 +3.4765625e-6,0.6034364049241694,0.5948453926137455 +3.486328125e-6,0.603393230077246,0.5949101548841306 +3.49609375e-6,0.6033505984177787,0.5949741023733315 +3.505859375e-6,0.6033085035235729,0.5950372447146403 +3.515625e-6,0.6032669389724339,0.5950995915413487 +3.525390625e-6,0.6032258983421669,0.5951611524867492 +3.53515625e-6,0.6031853752105775,0.5952219371841333 +3.544921875e-6,0.6031453631554708,0.5952819552667933 +3.5546875e-6,0.6031058557546523,0.5953412163680212 +3.564453125e-6,0.6030668465859271,0.5953997301211088 +3.57421875e-6,0.6030283292271009,0.5954575061593482 +3.583984375e-6,0.6029902972559787,0.5955145541160315 +3.59375e-6,0.602952744250366,0.5955708836244505 +3.603515625e-6,0.602915663788068,0.5956265043178974 +3.61328125e-6,0.6028790494468903,0.5956814258296641 +3.623046875e-6,0.6028428948165476,0.5957356577751781 +3.6328125e-6,0.6028071939742906,0.5957892090385637 +3.642578125e-6,0.6027719416370088,0.5958420875444865 +3.65234375e-6,0.6027371325644075,0.5958943011533883 +3.662109375e-6,0.6027027615161925,0.5959458577257108 +3.671875e-6,0.602668823252069,0.595996765121896 +3.681640625e-6,0.6026353125317427,0.5960470312023854 +3.69140625e-6,0.602602224114919,0.596096663827621 +3.701171875e-6,0.6025695527613033,0.5961456708580445 +3.7109375e-6,0.6025372932306012,0.5961940601540976 +3.720703125e-6,0.6025054402825182,0.5962418395762222 +3.73046875e-6,0.6024739886767596,0.5962890169848601 +3.740234375e-6,0.602442933173031,0.596335600240453 +3.75e-6,0.6024122685310379,0.5963815972034426 +3.759765625e-6,0.6023819895104857,0.5964270157342709 +3.76953125e-6,0.6023520908710799,0.5964718636933796 +3.779296875e-6,0.6023225673725261,0.5965161489412104 +3.7890625e-6,0.6022934139132782,0.5965598791300822 +3.798828125e-6,0.602264626115817,0.596603060826274 +3.80859375e-6,0.6022361997557296,0.596645700366405 +3.818359375e-6,0.6022081306024952,0.5966878040962567 +3.828125e-6,0.6021804144255926,0.5967293783616108 +3.837890625e-6,0.6021530469945005,0.5967704295082488 +3.84765625e-6,0.602126024078698,0.5968109638819525 +3.857421875e-6,0.6020993414476639,0.5968509878285037 +3.8671875e-6,0.6020729948708772,0.5968905076936838 +3.876953125e-6,0.6020469801178168,0.5969295298232744 +3.88671875e-6,0.6020212929579615,0.5969680605630574 +3.896484375e-6,0.6019959291607901,0.5970061062588142 +3.90625e-6,0.6019708844957818,0.5970436732563267 +3.916015625e-6,0.6019461547324154,0.5970807679013764 +3.92578125e-6,0.6019217356401697,0.5971173965397449 +3.935546875e-6,0.6018976229885237,0.597153565517214 +3.9453125e-6,0.6018738125470963,0.5971892811793551 +3.955078125e-6,0.6018503003370466,0.5972245494944298 +3.96484375e-6,0.6018270829074609,0.5972593756388083 +3.974609375e-6,0.6018041568419368,0.5972937647370945 +3.984375e-6,0.6017815187240715,0.5973277219138924 +3.994140625e-6,0.6017591651374624,0.5973612522938059 +4.00390625e-6,0.6017370926657071,0.597394361001439 +4.013671875e-6,0.6017152978924026,0.5974270531613955 +4.0234375e-6,0.6016937774011466,0.5974593338982797 +4.033203125e-6,0.6016725277755364,0.5974912083366951 +4.04296875e-6,0.6016515455991691,0.597522681601246 +4.052734375e-6,0.6016308274556423,0.5975537588165362 +4.0625e-6,0.6016103699285533,0.5975844451071696 +4.072265625e-6,0.6015901696014996,0.5976147455977502 +4.08203125e-6,0.6015702230580784,0.5976446654128821 +4.091796875e-6,0.6015505268818871,0.597674209677169 +4.1015625e-6,0.6015310776565231,0.5977033835152149 +4.111328125e-6,0.6015118719773678,0.597732192033948 +4.12109375e-6,0.6014929067719365,0.5977606398420948 +4.130859375e-6,0.6014741792754092,0.5977887310868858 +4.140625e-6,0.6014556867259054,0.5978164699111415 +4.150390625e-6,0.6014374263615446,0.5978438604576827 +4.16015625e-6,0.6014193954204461,0.5978709068693304 +4.169921875e-6,0.6014015911407294,0.5978976132889053 +4.1796875e-6,0.6013840107605142,0.5979239838592283 +4.189453125e-6,0.6013666515179196,0.5979500227231201 +4.19921875e-6,0.6013495106510653,0.5979757340234017 +4.208984375e-6,0.6013325853980707,0.5980011219028936 +4.21875e-6,0.6013158729970552,0.5980261905044167 +4.228515625e-6,0.6012993706861384,0.598050943970792 +4.23828125e-6,0.6012830757034396,0.5980753864448402 +4.248046875e-6,0.6012669852870783,0.5980995220693819 +4.2578125e-6,0.6012510966751742,0.5981233549872382 +4.267578125e-6,0.6012354071058464,0.5981468893412298 +4.27734375e-6,0.6012199138614217,0.5981701292078669 +4.287109375e-6,0.6012046145742815,0.5981930781385771 +4.296875e-6,0.601189507013407,0.5982157394798889 +4.306640625e-6,0.6011745889452552,0.5982381165821167 +4.31640625e-6,0.6011598581362833,0.5982602127955746 +4.326171875e-6,0.6011453123529482,0.5982820314705772 +4.3359375e-6,0.6011309493617071,0.5983035759574389 +4.345703125e-6,0.6011167669290172,0.5983248496064738 +4.35546875e-6,0.6011027628213353,0.5983458557679965 +4.365234375e-6,0.6010889348051189,0.5983665977923212 +4.375e-6,0.6010752806468248,0.5983870790297623 +4.384765625e-6,0.6010617981129103,0.5984073028306341 +4.39453125e-6,0.6010484849698323,0.5984272725452511 +4.404296875e-6,0.601035338984048,0.5984469915239276 +4.4140625e-6,0.6010223579220145,0.5984664631169777 +4.423828125e-6,0.601009539550189,0.5984856906747161 +4.43359375e-6,0.6009968816350284,0.5985046775474571 +4.443359375e-6,0.600984382038781,0.5985234269418281 +4.453125e-6,0.6009720389156844,0.598541941626473 +4.462890625e-6,0.6009598504588702,0.5985602243116943 +4.47265625e-6,0.6009478148613481,0.5985782777079774 +4.482421875e-6,0.600935930316128,0.5985961045258076 +4.4921875e-6,0.6009241950162197,0.5986137074756701 +4.501953125e-6,0.600912607154633,0.5986310892680502 +4.51171875e-6,0.6009011649243777,0.5986482526134331 +4.521484375e-6,0.6008898665184637,0.5986652002223042 +4.53125e-6,0.6008787101299006,0.5986819348051488 +4.541015625e-6,0.6008676939516985,0.5986984590724519 +4.55078125e-6,0.600856816176867,0.598714775734699 +4.560546875e-6,0.6008460749984161,0.5987308875023754 +4.5703125e-6,0.6008354686093557,0.5987467970859662 +4.580078125e-6,0.6008249952026953,0.5987625071959568 +4.58984375e-6,0.6008146529714448,0.5987780205428325 +4.599609375e-6,0.600804440111435,0.5987933398328471 +4.609375e-6,0.6007943549693945,0.598808467545908 +4.619140625e-6,0.6007843960771644,0.598823405884253 +4.62890625e-6,0.6007745619723416,0.5988381570414872 +4.638671875e-6,0.6007648511925228,0.5988527232112155 +4.6484375e-6,0.6007552622753046,0.5988671065870427 +4.658203125e-6,0.6007457937582841,0.5988813093625736 +4.66796875e-6,0.6007364441790577,0.598895333731413 +4.677734375e-6,0.6007272120752224,0.598909181887166 +4.6875e-6,0.6007180959843748,0.5989228560234374 +4.697265625e-6,0.6007090944441117,0.5989363583338321 +4.70703125e-6,0.6007002059920299,0.5989496910119548 +4.716796875e-6,0.600691429165726,0.5989628562514105 +4.7265625e-6,0.600682762502797,0.5989758562458041 +4.736328125e-6,0.6006742045408394,0.5989886931887404 +4.74609375e-6,0.6006657538174502,0.5990013692738243 +4.755859375e-6,0.600657408870226,0.5990138866946606 +4.765625e-6,0.6006491682316545,0.5990262476525179 +4.775390625e-6,0.6006410305338233,0.5990384541992646 +4.78515625e-6,0.6006329946021889,0.5990505080967162 +4.794921875e-6,0.6006250592761341,0.5990624110857984 +4.8046875e-6,0.600617223395042,0.5990741649074366 +4.814453125e-6,0.6006094857982954,0.5990857713025565 +4.82421875e-6,0.6006018453252774,0.5990972320120835 +4.833984375e-6,0.6005943008153708,0.5991085487769433 +4.84375e-6,0.6005868511079586,0.5991197233380616 +4.853515625e-6,0.6005794950424239,0.5991307574363638 +4.86328125e-6,0.6005722314581494,0.5991416528127754 +4.873046875e-6,0.6005650591945182,0.5991524112082223 +4.8828125e-6,0.6005579770909132,0.5991630343636298 +4.892578125e-6,0.6005509839867174,0.5991735240199236 +4.90234375e-6,0.6005440787213135,0.5991838819180292 +4.912109375e-6,0.6005372601340849,0.5991941097988722 +4.921875e-6,0.6005305270644142,0.5992042094033783 +4.931640625e-6,0.6005238783516844,0.5992141824724729 +4.94140625e-6,0.6005173128352786,0.5992240307470816 +4.951171875e-6,0.6005108293545796,0.5992337559681301 +4.9609375e-6,0.6005044267904684,0.5992433598142968 +4.970703125e-6,0.6004981042060347,0.5992528436909476 +4.98046875e-6,0.6004918606976344,0.5992622089535481 +4.990234375e-6,0.6004856953607576,0.5992714569588632 +5.e-6,0.6004796072908944,0.599280589063658 +5.009765625e-6,0.600473595583535,0.5992896066246971 +5.01953125e-6,0.6004676593341692,0.5992985109987458 +5.029296875e-6,0.6004617976382873,0.5993073035425686 +5.0390625e-6,0.6004560095913793,0.5993159856129308 +5.048828125e-6,0.6004502942889351,0.5993245585665969 +5.05859375e-6,0.600444650826445,0.5993330237603322 +5.068359375e-6,0.6004390782993989,0.5993413825509014 +5.078125e-6,0.6004335758032868,0.5993496362950694 +5.087890625e-6,0.600428142433599,0.5993577863496012 +5.09765625e-6,0.6004227772858253,0.5993658340712619 +5.107421875e-6,0.6004174794554559,0.5993737808168159 +5.1171875e-6,0.6004122480379808,0.5993816279430286 +5.126953125e-6,0.60040708212889,0.5993893768066646 +5.13671875e-6,0.6004019808236738,0.599397028764489 +5.146484375e-6,0.6003969432179087,0.5994045851731368 +5.15625e-6,0.6003919684738296,0.5994120472892553 +5.166015625e-6,0.6003870558808969,0.5994194161786545 +5.17578125e-6,0.6003822047344906,0.599426692898264 +5.185546875e-6,0.6003774143299906,0.5994338785050138 +5.1953125e-6,0.6003726839627773,0.5994409740558339 +5.205078125e-6,0.6003680129282307,0.5994479806076538 +5.21484375e-6,0.6003634005217308,0.5994548992174036 +5.224609375e-6,0.6003588460386577,0.5994617309420132 +5.234375e-6,0.6003543487743916,0.5994684768384123 +5.244140625e-6,0.6003499080243125,0.5994751379635309 +5.25390625e-6,0.6003455230838006,0.5994817153742988 +5.263671875e-6,0.6003411932482359,0.599488210127646 +5.2734375e-6,0.6003369178129985,0.5994946232805022 +5.283203125e-6,0.6003326960734684,0.5995009558897972 +5.29296875e-6,0.6003285273250258,0.599507209012461 +5.302734375e-6,0.6003244108630509,0.5995133837054235 +5.3125e-6,0.6003203459829236,0.5995194810256144 +5.322265625e-6,0.600316331980024,0.5995255020299637 +5.33203125e-6,0.6003123681497323,0.5995314477754012 +5.341796875e-6,0.6003084537900534,0.5995373193149197 +5.3515625e-6,0.6003045882794622,0.5995431175808065 +5.361328125e-6,0.600300771069444,0.5995488433958338 +5.37109375e-6,0.600297001611059,0.5995544975834112 +5.380859375e-6,0.6002932793553678,0.5995600809669481 +5.390625e-6,0.6002896037534304,0.599565594369854 +5.400390625e-6,0.6002859742563075,0.5995710386155385 +5.41015625e-6,0.6002823903150591,0.599576414527411 +5.419921875e-6,0.6002788513807459,0.5995817229288809 +5.4296875e-6,0.6002753569044279,0.5995869646433579 +5.439453125e-6,0.6002719063371655,0.5995921404942514 +5.44921875e-6,0.6002684991300193,0.5995972513049708 +5.458984375e-6,0.6002651347340494,0.5996022978989257 +5.46875e-6,0.6002618126003161,0.5996072810995255 +5.478515625e-6,0.6002585321798799,0.5996122017301798 +5.48828125e-6,0.6002552929238011,0.5996170606142981 +5.498046875e-6,0.60025209428314,0.5996218585752897 +5.5078125e-6,0.6002489357089569,0.5996265964365644 +5.517578125e-6,0.6002458166523122,0.5996312750215314 +5.52734375e-6,0.6002427365642663,0.5996358951536003 +5.537109375e-6,0.6002396949043242,0.5996404576435135 +5.546875e-6,0.6002366912120876,0.5996449631818683 +5.556640625e-6,0.6002337250606555,0.5996494124090165 +5.56640625e-6,0.6002307960220113,0.5996538059669826 +5.576171875e-6,0.6002279036681387,0.5996581444977916 +5.5859375e-6,0.6002250475710208,0.5996624286434685 +5.595703125e-6,0.6002222273026411,0.5996666590460381 +5.60546875e-6,0.6002194424349829,0.5996708363475252 +5.615234375e-6,0.60021669254003,0.5996749611899547 +5.625e-6,0.6002139771897654,0.5996790342153515 +5.634765625e-6,0.6002112959561728,0.5996830560657403 +5.64453125e-6,0.6002086484112356,0.5996870273831463 +5.654296875e-6,0.6002060341269371,0.5996909488095941 +5.6640625e-6,0.6002034526752608,0.5996948209871086 +5.673828125e-6,0.60020090362819,0.5996986445577147 +5.68359375e-6,0.6001983865577083,0.5997024201634372 +5.693359375e-6,0.600195901035799,0.5997061484463012 +5.703125e-6,0.6001934466344456,0.5997098300483312 +5.712890625e-6,0.6001910229256315,0.5997134656115525 +5.72265625e-6,0.6001886294813401,0.5997170557779895 +5.732421875e-6,0.6001862658904241,0.5997206011643635 +5.7421875e-6,0.6001839318093617,0.5997241022859572 +5.751953125e-6,0.6001816269050002,0.5997275596424992 +5.76171875e-6,0.6001793508439062,0.5997309737341404 +5.771484375e-6,0.6001771032926457,0.5997343450610311 +5.78125e-6,0.6001748839177854,0.5997376741233215 +5.791015625e-6,0.6001726923858915,0.5997409614211624 +5.80078125e-6,0.6001705283635304,0.599744207454704 +5.810546875e-6,0.6001683915172684,0.599747412724097 +5.8203125e-6,0.6001662815136719,0.5997505777294918 +5.830078125e-6,0.6001641980193073,0.5997537029710387 +5.83984375e-6,0.6001621407007409,0.5997567889488883 +5.849609375e-6,0.6001601092245391,0.5997598361631911 +5.859375e-6,0.6001581032572681,0.5997628451140973 +5.869140625e-6,0.6001561224654947,0.5997658163017576 +5.87890625e-6,0.6001541665157847,0.5997687502263226 +5.888671875e-6,0.6001522350747048,0.5997716473879423 +5.8984375e-6,0.6001503278088214,0.5997745082867676 +5.908203125e-6,0.6001484443847006,0.5997773334229487 +5.91796875e-6,0.6001465844689712,0.5997801232965428 +5.927734375e-6,0.6001447477537022,0.5997828783694463 +5.9375e-6,0.6001429339768508,0.5997855990347234 +5.947265625e-6,0.6001411428784529,0.5997882856823202 +5.95703125e-6,0.6001393741985442,0.5997909387021831 +5.966796875e-6,0.6001376276771608,0.5997935584842584 +5.9765625e-6,0.6001359030543383,0.5997961454184922 +5.986328125e-6,0.6001342000701125,0.5997986998948308 +5.99609375e-6,0.6001325184645193,0.5998012223032206 +6.005859375e-6,0.6001308579775946,0.5998037130336077 +6.015625e-6,0.600129218349374,0.5998061724759385 +6.025390625e-6,0.6001275993198936,0.5998086010201591 +6.03515625e-6,0.6001260006291891,0.5998109990562159 +6.044921875e-6,0.6001244220172962,0.5998133669740552 +6.0546875e-6,0.6001228632242509,0.599815705163623 +6.064453125e-6,0.600121323990089,0.5998180140148659 +6.07421875e-6,0.6001198040548463,0.59982029391773 +6.083984375e-6,0.6001183031585586,0.5998225452621615 +6.09375e-6,0.6001168210412618,0.5998247684381067 +6.103515625e-6,0.6001153574429917,0.5998269638355119 +6.11328125e-6,0.6001139121027919,0.5998291318458117 +6.123046875e-6,0.6001124847665015,0.5998312728502473 +6.1328125e-6,0.6001110752280886,0.5998333871578665 +6.142578125e-6,0.6001096832900119,0.5998354750649817 +6.15234375e-6,0.6001083087547294,0.5998375368679053 +6.162109375e-6,0.6001069514246997,0.59983957286295 +6.171875e-6,0.6001056111023809,0.599841583346428 +6.181640625e-6,0.6001042875902316,0.5998435686146519 +6.19140625e-6,0.6001029806907101,0.5998455289639344 +6.201171875e-6,0.6001016902062746,0.5998474646905876 +6.2109375e-6,0.6001004159393835,0.5998493760909243 +6.220703125e-6,0.6000991576924951,0.5998512634612567 +6.23046875e-6,0.600097915268068,0.5998531270978975 +6.240234375e-6,0.6000966884685602,0.5998549672971591 +6.25e-6,0.6000954770964303,0.5998567843553541 +6.259765625e-6,0.6000942809541365,0.5998585785687948 +6.26953125e-6,0.6000930998441372,0.5998603502337937 +6.279296875e-6,0.6000919335688908,0.5998620996466633 +6.2890625e-6,0.6000907819308555,0.5998638271037162 +6.298828125e-6,0.6000896447324898,0.5998655329012648 +6.30859375e-6,0.6000885217762519,0.5998672173356215 +6.318359375e-6,0.6000874128646003,0.599868880703099 +6.328125e-6,0.6000863177999933,0.5998705233000095 +6.337890625e-6,0.6000852363848892,0.5998721454226656 +6.34765625e-6,0.6000841684340893,0.5998737473488656 +6.357421875e-6,0.6000831138004324,0.5998753292993508 +6.3671875e-6,0.6000820723362431,0.5998768914956348 +6.376953125e-6,0.6000810438935925,0.5998784341596107 +6.38671875e-6,0.6000800283245521,0.5998799575131714 +6.396484375e-6,0.6000790254811932,0.5998814617782098 +6.40625e-6,0.6000780352155871,0.5998829471766188 +6.416015625e-6,0.6000770573798052,0.5998844139302917 +6.42578125e-6,0.6000760918259189,0.5998858622611212 +6.435546875e-6,0.6000751384059995,0.5998872923910002 +6.4453125e-6,0.6000741969721184,0.599888704541822 +6.455078125e-6,0.6000732673763468,0.5998900989354793 +6.46484375e-6,0.6000723494707563,0.5998914757938651 +6.474609375e-6,0.600071443107418,0.5998928353388725 +6.484375e-6,0.6000705481384034,0.5998941777923944 +6.494140625e-6,0.6000696644157838,0.5998955033763237 +6.50390625e-6,0.6000687917916306,0.5998968123125535 +6.513671875e-6,0.6000679301180152,0.5998981048229768 +6.5234375e-6,0.6000670792470088,0.5998993811294864 +6.533203125e-6,0.6000662390306828,0.5999006414539754 +6.54296875e-6,0.6000654093211085,0.5999018860183366 +6.552734375e-6,0.6000645899703575,0.5999031150444633 +6.5625e-6,0.6000637808305009,0.5999043287542482 +6.572265625e-6,0.6000629817557825,0.5999055273663257 +6.58203125e-6,0.6000621926371683,0.599906711044247 +6.591796875e-6,0.6000614133698277,0.5999078799452581 +6.6015625e-6,0.6000606438432212,0.5999090342351676 +6.611328125e-6,0.60005988394681,0.5999101740797845 +6.62109375e-6,0.6000591335700547,0.5999112996449174 +6.630859375e-6,0.6000583926024164,0.5999124110963748 +6.640625e-6,0.6000576609333559,0.5999135085999657 +6.650390625e-6,0.6000569384523339,0.5999145923214986 +6.66015625e-6,0.6000562250488115,0.5999156624267822 +6.669921875e-6,0.6000555206122494,0.5999167190816252 +6.6796875e-6,0.6000548250321086,0.5999177624518365 +6.689453125e-6,0.6000541381978498,0.5999187927032247 +6.69921875e-6,0.6000534599989341,0.5999198100015983 +6.708984375e-6,0.6000527903248221,0.5999208145127662 +6.71875e-6,0.6000521290649748,0.599921806402537 +6.728515625e-6,0.6000514761088531,0.5999227858367195 +6.73828125e-6,0.6000508313459179,0.5999237529811224 +6.748046875e-6,0.6000501946656299,0.5999247080015544 +6.7578125e-6,0.6000495659574501,0.5999256510638241 +6.767578125e-6,0.6000489451108393,0.5999265823337403 +6.77734375e-6,0.6000483320152584,0.5999275019771116 +6.787109375e-6,0.6000477265601682,0.5999284101597468 +6.796875e-6,0.6000471286350297,0.5999293070474546 +6.806640625e-6,0.6000465381433459,0.5999301927849805 +6.81640625e-6,0.6000459550161507,0.5999310674757732 +6.826171875e-6,0.6000453791718565,0.5999319312422146 +6.8359375e-6,0.6000448105284687,0.5999327842072961 +6.845703125e-6,0.6000442490039928,0.59993362649401 +6.85546875e-6,0.6000436945164342,0.599934458225348 +6.865234375e-6,0.6000431469837985,0.5999352795243016 +6.875e-6,0.6000426063240909,0.599936090513863 +6.884765625e-6,0.600042072455317,0.5999368913170238 +6.89453125e-6,0.6000415452954824,0.5999376820567758 +6.904296875e-6,0.6000410247625922,0.5999384628561109 +6.9140625e-6,0.6000405107746523,0.5999392338380208 +6.923828125e-6,0.6000400032496679,0.5999399951254974 +6.93359375e-6,0.6000395021056445,0.5999407468415325 +6.943359375e-6,0.6000390072605875,0.599941489109118 +6.953125e-6,0.6000385186325025,0.5999422220512456 +6.962890625e-6,0.6000380361393948,0.5999429457909071 +6.97265625e-6,0.60003755969927,0.5999436604510944 +6.982421875e-6,0.6000370892301334,0.5999443661547991 +6.9921875e-6,0.6000366246499906,0.5999450630250133 +7.001953125e-6,0.6000361658768472,0.5999457511847286 +7.01171875e-6,0.6000357128287083,0.5999464307569369 +7.021484375e-6,0.6000352654235795,0.59994710186463 +7.03125e-6,0.6000348235809146,0.5999477646286272 +7.041015625e-6,0.6000343872433481,0.599948419134977 +7.05078125e-6,0.6000339563538089,0.5999490654692861 +7.060546875e-6,0.6000335308511676,0.5999497037232479 +7.0703125e-6,0.6000331106742953,0.5999503339885562 +7.080078125e-6,0.6000326957620629,0.599950956356905 +7.08984375e-6,0.600032286053341,0.5999515709199876 +7.099609375e-6,0.6000318814870008,0.5999521777694979 +7.109375e-6,0.600031482001913,0.5999527769971297 +7.119140625e-6,0.6000310875369484,0.5999533686945764 +7.12890625e-6,0.6000306980309781,0.599953952953532 +7.138671875e-6,0.6000303134228727,0.59995452986569 +7.1484375e-6,0.6000299336515033,0.5999550995227443 +7.158203125e-6,0.6000295586557406,0.5999556620163883 +7.16796875e-6,0.6000291883744554,0.599956217438316 +7.177734375e-6,0.6000288227465188,0.599956765880221 +7.1875e-6,0.6000284617108015,0.5999573074337968 +7.197265625e-6,0.6000281052061744,0.5999578421907374 +7.20703125e-6,0.6000277531715085,0.5999583702427365 +7.216796875e-6,0.6000274055456744,0.5999588916814875 +7.2265625e-6,0.6000270622675431,0.5999594065986844 +7.236328125e-6,0.6000267232759856,0.5999599150860208 +7.24609375e-6,0.6000263885098726,0.5999604172351902 +7.255859375e-6,0.600026057908075,0.5999609131378867 +7.265625e-6,0.6000257314174525,0.5999614028738204 +7.275390625e-6,0.6000254090001699,0.5999618864997442 +7.28515625e-6,0.6000250906110025,0.5999623640834953 +7.294921875e-6,0.6000247762044956,0.5999628356932558 +7.3046875e-6,0.6000244657351941,0.5999633013972079 +7.314453125e-6,0.6000241591576434,0.599963761263534 +7.32421875e-6,0.6000238564263884,0.5999642153604163 +7.333984375e-6,0.6000235574959746,0.5999646637560372 +7.34375e-6,0.6000232623209468,0.5999651065185789 +7.353515625e-6,0.6000229708558502,0.5999655437162237 +7.36328125e-6,0.6000226830552302,0.5999659754171538 +7.373046875e-6,0.6000223988736316,0.5999664016895516 +7.3828125e-6,0.6000221182655998,0.5999668226015994 +7.392578125e-6,0.6000218411856798,0.5999672382214793 +7.40234375e-6,0.6000215675884167,0.5999676486173738 +7.412109375e-6,0.6000212974283559,0.5999680538574652 +7.421875e-6,0.6000210306600423,0.5999684540099355 +7.431640625e-6,0.6000207672380212,0.5999688491429672 +7.44140625e-6,0.6000205071168376,0.5999692393247427 +7.451171875e-6,0.6000202502510368,0.5999696246234439 +7.4609375e-6,0.6000199965951637,0.5999700051072534 +7.470703125e-6,0.6000197461037637,0.5999703808443534 +7.48046875e-6,0.6000194987313819,0.5999707519029261 +7.490234375e-6,0.600019254433257,0.5999711183501135 +7.5e-6,0.600019013176347,0.5999714802354785 +7.509765625e-6,0.6000187749291618,0.5999718376062563 +7.51953125e-6,0.6000185396584508,0.5999721905123229 +7.529296875e-6,0.600018307330963,0.5999725390035545 +7.5390625e-6,0.6000180779134476,0.5999728831298275 +7.548828125e-6,0.600017851372654,0.5999732229410181 +7.55859375e-6,0.6000176276753312,0.5999735584870022 +7.568359375e-6,0.6000174067882285,0.5999738898176563 +7.578125e-6,0.600017188678095,0.5999742169828566 +7.587890625e-6,0.6000169733116799,0.5999745400324792 +7.59765625e-6,0.6000167606557325,0.5999748590164002 +7.607421875e-6,0.6000165506770019,0.5999751739844962 +7.6171875e-6,0.6000163433422374,0.599975484986643 +7.626953125e-6,0.600016138618188,0.5999757920727169 +7.63671875e-6,0.6000159364716031,0.5999760952925943 +7.646484375e-6,0.6000157368692318,0.5999763946961513 +7.65625e-6,0.6000155397778234,0.599976690333264 +7.666015625e-6,0.6000153451641269,0.5999769822538087 +7.67578125e-6,0.6000151529948916,0.5999772705076616 +7.685546875e-6,0.6000149632368668,0.599977555144699 +7.6953125e-6,0.6000147758568014,0.5999778362147968 +7.705078125e-6,0.6000145908214449,0.5999781137678316 +7.71484375e-6,0.6000144080975464,0.5999783878536795 +7.724609375e-6,0.600014227651855,0.5999786585222165 +7.734375e-6,0.60001404945112,0.599978925823319 +7.744140625e-6,0.6000138734620906,0.5999791898068632 +7.75390625e-6,0.6000136996515159,0.5999794505227252 +7.763671875e-6,0.6000135279925198,0.5999797080112191 +7.7734375e-6,0.6000133584710301,0.5999799622934539 +7.783203125e-6,0.6000131910641018,0.5999802134038463 +7.79296875e-6,0.6000130257483776,0.5999804613774326 +7.802734375e-6,0.6000128625005002,0.5999807062492487 +7.8125e-6,0.6000127012971123,0.5999809480543307 +7.822265625e-6,0.6000125421148563,0.5999811868277145 +7.83203125e-6,0.6000123849303751,0.5999814226044363 +7.841796875e-6,0.6000122297203113,0.599981655419532 +7.8515625e-6,0.6000120764613076,0.5999818853080376 +7.861328125e-6,0.6000119251300066,0.5999821123049892 +7.87109375e-6,0.6000117757030508,0.5999823364454228 +7.880859375e-6,0.6000116281570831,0.5999825577643744 +7.890625e-6,0.6000114824687459,0.59998277629688 +7.900390625e-6,0.6000113386146821,0.5999829920779757 +7.91015625e-6,0.6000111965715343,0.5999832051426974 +7.919921875e-6,0.6000110563159451,0.5999834155260814 +7.9296875e-6,0.6000109178245571,0.5999836232631633 +7.939453125e-6,0.600010781074013,0.5999838283889795 +7.94921875e-6,0.6000106460409554,0.5999840309385658 +7.958984375e-6,0.6000105127020271,0.5999842309469583 +7.96875e-6,0.6000103810338706,0.599984428449193 +7.978515625e-6,0.6000102510131287,0.5999846234803059 +7.98828125e-6,0.6000101226164439,0.599984816075333 +7.998046875e-6,0.6000099958204589,0.5999850062693105 +8.0078125e-6,0.6000098706018164,0.5999851940972742 +8.017578125e-6,0.600009746937159,0.5999853795942603 +8.02734375e-6,0.6000096248034474,0.5999855627948278 +8.037109375e-6,0.6000095041873517,0.5999857437189713 +8.046875e-6,0.6000093850753333,0.5999859223869989 +8.056640625e-6,0.6000092674508595,0.5999860988237097 +8.06640625e-6,0.6000091512973974,0.5999862730539027 +8.076171875e-6,0.6000090365984145,0.5999864451023772 +8.0859375e-6,0.6000089233373779,0.5999866149939321 +8.095703125e-6,0.6000088114977548,0.5999867827533667 +8.10546875e-6,0.6000087010630126,0.5999869484054801 +8.115234375e-6,0.6000085920166185,0.5999871119750713 +8.125e-6,0.6000084843420397,0.5999872734869394 +8.134765625e-6,0.6000083780227435,0.5999874329658836 +8.14453125e-6,0.6000082730421972,0.5999875904367031 +8.154296875e-6,0.6000081693838679,0.599987745924197 +8.1640625e-6,0.6000080670312231,0.5999878994531642 +8.173828125e-6,0.6000079659677299,0.599988051048404 +8.18359375e-6,0.6000078661768555,0.5999882007347156 +8.193359375e-6,0.6000077676420673,0.5999883485368979 +8.203125e-6,0.6000076703468326,0.5999884944797501 +8.212890625e-6,0.6000075742746184,0.5999886385880714 +8.22265625e-6,0.6000074794088921,0.5999887808866609 +8.232421875e-6,0.600007385733121,0.5999889214003175 +8.2421875e-6,0.6000072932307723,0.5999890601538406 +8.251953125e-6,0.6000072018853132,0.5999891971720291 +8.26171875e-6,0.600007111680211,0.5999893324796823 +8.271484375e-6,0.6000070225989331,0.5999894661015992 +8.28125e-6,0.6000069346249466,0.599989598062579 +8.291015625e-6,0.6000068477417188,0.5999897283874207 +8.30078125e-6,0.6000067619349286,0.5999898570976059 +8.310546875e-6,0.6000066771973394,0.59998998420399 +8.3203125e-6,0.6000065935176556,0.5999901097235155 +8.330078125e-6,0.6000065108841406,0.599990233673788 +8.33984375e-6,0.6000064292850575,0.5999903560724127 +8.349609375e-6,0.6000063487086693,0.599990476936995 +8.359375e-6,0.6000062691432392,0.5999905962851402 +8.369140625e-6,0.6000061905770303,0.5999907141344535 +8.37890625e-6,0.6000061129983056,0.5999908305025405 +8.388671875e-6,0.6000060363953285,0.5999909454070064 +8.3984375e-6,0.6000059607563616,0.5999910588654566 +8.408203125e-6,0.6000058860696685,0.5999911708954963 +8.41796875e-6,0.600005812323512,0.599991281514731 +8.427734375e-6,0.6000057395061554,0.599991390740766 +8.4375e-6,0.6000056676058617,0.5999914985912066 +8.447265625e-6,0.600005596610894,0.5999916050836581 +8.45703125e-6,0.6000055265095154,0.599991710235726 +8.466796875e-6,0.6000054572899891,0.5999918140650154 +8.4765625e-6,0.6000053889405782,0.5999919165891319 +8.486328125e-6,0.6000053214495457,0.5999920178256807 +8.49609375e-6,0.6000052548051547,0.5999921177922671 +8.505859375e-6,0.6000051889956685,0.5999922165064965 +8.515625e-6,0.6000051240093499,0.5999923139859743 +8.525390625e-6,0.6000050598344623,0.5999924102483057 +8.53515625e-6,0.6000049964592686,0.5999925053110962 +8.544921875e-6,0.6000049338720321,0.599992599191951 +8.5546875e-6,0.6000048720610158,0.5999926919084755 +8.564453125e-6,0.6000048110145273,0.5999927834782082 +8.57421875e-6,0.6000047507255325,0.5999928739117004 +8.583984375e-6,0.6000046911879144,0.5999929632181275 +8.59375e-6,0.6000046323934365,0.5999930514098445 +8.603515625e-6,0.6000045743338615,0.5999931384992069 +8.61328125e-6,0.6000045170009528,0.5999932244985701 +8.623046875e-6,0.6000044603864734,0.5999933094202892 +8.6328125e-6,0.6000044044821864,0.5999933932767196 +8.642578125e-6,0.6000043492798549,0.5999934760802167 +8.65234375e-6,0.6000042947712422,0.5999935578431358 +8.662109375e-6,0.6000042409481112,0.5999936385778323 +8.671875e-6,0.6000041878022252,0.5999937182966614 +8.681640625e-6,0.6000041353253471,0.5999937970119784 +8.69140625e-6,0.6000040835092403,0.5999938747361387 +8.701171875e-6,0.6000040323456677,0.5999939514814977 +8.7109375e-6,0.6000039818263924,0.5999940272604105 +8.720703125e-6,0.6000039319431776,0.5999941020852327 +8.73046875e-6,0.6000038826877865,0.5999941759683194 +8.740234375e-6,0.600003834051982,0.5999942489220261 +8.75e-6,0.6000037860275275,0.599994320958708 +8.759765625e-6,0.6000037386061858,0.5999943920907205 +8.76953125e-6,0.6000036917797202,0.5999944623304188 +8.779296875e-6,0.6000036455398938,0.5999945316901584 +8.7890625e-6,0.6000035998784697,0.5999946001822944 +8.798828125e-6,0.6000035547872111,0.5999946678191824 +8.80859375e-6,0.600003510257881,0.5999947346131775 +8.818359375e-6,0.6000034662822425,0.5999948005766352 +8.828125e-6,0.6000034228520589,0.5999948657219107 +8.837890625e-6,0.600003379959807,0.5999949300602886 +8.84765625e-6,0.6000033376016136,0.5999949935975787 +8.857421875e-6,0.6000032957720034,0.599995056341994 +8.8671875e-6,0.6000032544651398,0.5999951183022894 +8.876953125e-6,0.6000032136751867,0.599995179487219 +8.88671875e-6,0.6000031733963076,0.5999952399055376 +8.896484375e-6,0.6000031336226662,0.5999952995659997 +8.90625e-6,0.6000030943484261,0.59999535847736 +8.916015625e-6,0.6000030555677509,0.5999954166483726 +8.92578125e-6,0.6000030172748044,0.5999954740877925 +8.935546875e-6,0.60000297946375,0.599995530804374 +8.9453125e-6,0.6000029421287516,0.5999955868068716 +8.955078125e-6,0.6000029052639727,0.59999564210404 +8.96484375e-6,0.6000028688635769,0.5999956967046337 +8.974609375e-6,0.600002832921728,0.5999957506174072 +8.984375e-6,0.6000027974325894,0.599995803851115 +8.994140625e-6,0.600002762390325,0.5999958564145116 +9.00390625e-6,0.6000027277890982,0.5999959083163516 +9.013671875e-6,0.6000026936230729,0.5999959595653896 +9.0234375e-6,0.6000026598864125,0.5999960101703802 +9.033203125e-6,0.6000026265732809,0.5999960601400778 +9.04296875e-6,0.6000025936778415,0.5999961094832369 +9.052734375e-6,0.600002561194258,0.5999961582086121 +9.0625e-6,0.600002529116694,0.599996206324958 +9.072265625e-6,0.6000024974393134,0.5999962538410291 +9.08203125e-6,0.6000024661562795,0.5999963007655797 +9.091796875e-6,0.6000024352617562,0.5999963471073648 +9.1015625e-6,0.600002404749907,0.5999963928751386 +9.111328125e-6,0.6000023746162494,0.5999964380756251 +9.12109375e-6,0.600002344857503,0.5999964827137446 +9.130859375e-6,0.6000023154696673,0.5999965267954982 +9.140625e-6,0.6000022864487415,0.5999965703268868 +9.150390625e-6,0.600002257790725,0.5999966133139115 +9.16015625e-6,0.6000022294916172,0.5999966557625733 +9.169921875e-6,0.6000022015474172,0.5999966976788733 +9.1796875e-6,0.6000021739541245,0.5999967390688125 +9.189453125e-6,0.6000021467077382,0.5999967799383917 +9.19921875e-6,0.6000021198042579,0.5999968202936122 +9.208984375e-6,0.6000020932396828,0.5999968601404749 +9.21875e-6,0.6000020670100121,0.5999968994849809 +9.228515625e-6,0.6000020411112453,0.5999969383331311 +9.23828125e-6,0.6000020155393817,0.5999969766909266 +9.248046875e-6,0.6000019902904205,0.5999970145643684 +9.2578125e-6,0.6000019653603611,0.5999970519594575 +9.267578125e-6,0.6000019407452027,0.5999970888821949 +9.27734375e-6,0.6000019164409449,0.5999971253385817 +9.287109375e-6,0.6000018924435867,0.5999971613346189 +9.296875e-6,0.6000018687491276,0.5999971968763075 +9.306640625e-6,0.600001845353567,0.5999972319696486 +9.31640625e-6,0.600001822252904,0.599997266620643 +9.326171875e-6,0.600001799443138,0.5999973008352919 +9.3359375e-6,0.6000017769202685,0.5999973346195964 +9.345703125e-6,0.6000017546802945,0.5999973679795573 +9.35546875e-6,0.6000017327192155,0.5999974009211757 +9.365234375e-6,0.6000017110330308,0.5999974334504528 +9.375e-6,0.6000016896177398,0.5999974655733893 +9.384765625e-6,0.6000016684693417,0.5999974972959865 +9.39453125e-6,0.6000016475838358,0.5999975286242453 +9.404296875e-6,0.6000016269572216,0.5999975595641667 +9.4140625e-6,0.6000016065854982,0.5999975901217518 +9.423828125e-6,0.6000015864646651,0.5999976203030015 +9.43359375e-6,0.6000015665907213,0.5999976501139169 +9.443359375e-6,0.6000015469596667,0.5999976795604991 +9.453125e-6,0.6000015275677082,0.5999977086484367 +9.462890625e-6,0.6000015084134519,0.5999977373798211 +9.47265625e-6,0.6000014894949817,0.5999977657575265 +9.482421875e-6,0.6000014708099507,0.599997793785073 +9.4921875e-6,0.6000014523560121,0.5999978214659809 +9.501953125e-6,0.6000014341308192,0.5999978488037703 +9.51171875e-6,0.6000014161320251,0.5999978758019614 +9.521484375e-6,0.6000013983572832,0.5999979024640743 +9.53125e-6,0.6000013808042465,0.5999979287936292 +9.541015625e-6,0.6000013634705684,0.5999979547941464 +9.55078125e-6,0.6000013463539021,0.5999979804691459 +9.560546875e-6,0.6000013294519008,0.5999980058221479 +9.5703125e-6,0.6000013127622176,0.5999980308566727 +9.580078125e-6,0.6000012962825059,0.5999980555762403 +9.58984375e-6,0.6000012800104187,0.5999980799843709 +9.599609375e-6,0.6000012639436094,0.5999981040845849 +9.609375e-6,0.6000012480797312,0.5999981278804022 +9.619140625e-6,0.6000012324164373,0.5999981513753431 +9.62890625e-6,0.600001216951381,0.5999981745729277 +9.638671875e-6,0.6000012016822153,0.5999981974766762 +9.6484375e-6,0.6000011866065935,0.5999982200901087 +9.658203125e-6,0.6000011717221689,0.5999982424167456 +9.66796875e-6,0.6000011570265948,0.5999982644601068 +9.677734375e-6,0.6000011425175242,0.5999982862237128 +9.6875e-6,0.6000011281926105,0.5999983077110834 +9.697265625e-6,0.6000011140495067,0.5999983289257389 +9.70703125e-6,0.6000011000858664,0.5999983498711996 +9.716796875e-6,0.6000010862993423,0.5999983705509856 +9.7265625e-6,0.6000010726875165,0.5999983909687244 +9.736328125e-6,0.6000010592456387,0.599998411131541 +9.74609375e-6,0.6000010459706229,0.5999984310440648 +9.755859375e-6,0.6000010328607138,0.5999984507089283 +9.765625e-6,0.6000010199141569,0.5999984701287638 +9.775390625e-6,0.6000010071291969,0.5999984893062037 +9.78515625e-6,0.6000009945040792,0.5999985082438805 +9.794921875e-6,0.6000009820370485,0.5999985269444262 +9.8046875e-6,0.6000009697263504,0.5999985454104735 +9.814453125e-6,0.6000009575702296,0.5999985636446546 +9.82421875e-6,0.6000009455669313,0.5999985816496021 +9.833984375e-6,0.6000009337147006,0.5999985994279481 +9.84375e-6,0.6000009220117827,0.5999986169823252 +9.853515625e-6,0.6000009104564223,0.5999986343153655 +9.86328125e-6,0.6000008990468649,0.5999986514297017 +9.873046875e-6,0.6000008877813555,0.5999986683279658 +9.8828125e-6,0.6000008766581391,0.5999986850127905 +9.892578125e-6,0.6000008656754607,0.599998701486808 +9.90234375e-6,0.6000008548315656,0.5999987177526507 +9.912109375e-6,0.6000008441246987,0.599998733812951 +9.921875e-6,0.6000008335531052,0.5999987496703413 +9.931640625e-6,0.6000008231150301,0.5999987653274538 +9.94140625e-6,0.6000008128087186,0.599998780786921 +9.951171875e-6,0.6000008026324157,0.5999987960513754 +9.9609375e-6,0.6000007925843666,0.5999988111234492 +9.970703125e-6,0.6000007826628162,0.5999988260057747 +9.98046875e-6,0.6000007728660097,0.5999988407009845 +9.990234375e-6,0.6000007631921922,0.5999988552117108 +0.00001,0.6000007536396087,0.599998869540586 From 063a46ae7b2828f339659ef0503a767d3fe75035 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 15 Aug 2024 20:53:14 +1000 Subject: [PATCH 08/25] add dust temperature coupled with gas and radiation --- src/radiation/radiation_system.hpp | 133 +++++++++++++++++++++++------ 1 file changed, 107 insertions(+), 26 deletions(-) diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 202c4213f..9f055598f 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -45,7 +45,7 @@ static const int max_ite_to_update_alpha_E = 5; // Apply to the PPL_opacity_full // iterations of the Newton iteration static constexpr bool include_work_term_in_source = true; -static constexpr bool use_D_as_base = true; +static constexpr bool use_D_as_base = false; static const bool PPL_free_slope_st_total = false; // PPL with free slopes for all, but subject to the constraint sum_g alpha_g B_g = - sum_g B_g. Not working // well -- Newton iteration convergence issue. @@ -62,6 +62,10 @@ static constexpr double c_light_cgs_ = C::c_light; // cgs static constexpr double radiation_constant_cgs_ = C::a_rad; // cgs static constexpr double inf = std::numeric_limits::max(); +static constexpr bool use_dust = true; +static const double Lambda_gd_0 = 1e6; +// static const double Lambda_gd_0 = 2.5e-34; // erg cm^3 s^−1 K^−3/2 + // enum for opacity_model enum class OpacityModel { single_group = 0, // user-defined opacity for each group, given as a function of density and temperature. @@ -174,7 +178,7 @@ template class RadSystem : public HyperbolicSystem::mean_molecular_mass; + static constexpr double mean_molecular_mass_ = quokka::EOS_Traits::mean_molecular_weight; static constexpr double boltzmann_constant_ = quokka::EOS_Traits::boltzmann_constant; static constexpr double gamma_ = quokka::EOS_Traits::gamma; @@ -262,6 +266,10 @@ template class RadSystem : public HyperbolicSystem const &boundaries) -> quokka::valarray; + AMREX_GPU_HOST_DEVICE static auto + ComputeDustTemperature(double T_gas, double T_d_init, double rho, quokka::valarray const &Erad, + amrex::GpuArray const &rad_boundaries) -> double; + template AMREX_GPU_DEVICE static auto ComputeCellOpticalDepth(const quokka::Array4View &consVar, amrex::GpuArray dx, int i, int j, int k, @@ -309,8 +317,8 @@ template AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiation(amrex::Real temperature, amrex::GpuArray const &boundaries) -> quokka::valarray { - auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); - double power = radiation_constant_ * std::pow(temperature, 4); + const auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); + const double power = radiation_constant_ * std::pow(temperature, 4); auto Erad_g = power * radEnergyFractions; // set floor for (int g = 0; g < nGroups_; ++g) { @@ -327,8 +335,9 @@ RadSystem::ComputeThermalRadiationTempDerivative(amrex::Real temperat amrex::GpuArray const &boundaries) -> quokka::valarray { // by default, d emission/dT = 4 emission / T - auto erad = ComputeThermalRadiation(temperature, boundaries); - return 4. * erad / temperature; + auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); + double d_power_dt = 4. * radiation_constant_ * std::pow(temperature, 3); + return d_power_dt * radEnergyFractions; } // Linear equation solver for matrix with non-zeros at the first row, first column, and diagonal only. @@ -1169,6 +1178,51 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxInDiffusionLimit(con return flux; } +template +AMREX_GPU_HOST_DEVICE auto +RadSystem::ComputeDustTemperature(double const T_gas, double const T_d_init, double const rho, quokka::valarray const &Erad, + amrex::GpuArray const &rad_boundaries) -> double +{ + quokka::valarray fourPiBoverC{}; + quokka::valarray chiPVec{}; + quokka::valarray chiEVec{}; + + const double num_density = rho / mean_molecular_mass_; + const double Lambda_compare = Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas) * T_gas; + + // solve for dust temperature: Newton iteration on the dust temperature T_d + double T_d = T_d_init; + const double lambda_rel_tol = 1.0e-8; + const int max_ite_td = 100; + int ite_td = 0; + double delta_T_d = NAN; + for (; ite_td < max_ite_td; ++ite_td) { + chiPVec[0] = ComputePlanckOpacity(rho, T_d) * rho; + chiEVec[0] = ComputeEnergyMeanOpacity(rho, T_d) * rho; + + fourPiBoverC = ComputeThermalRadiation(T_d, rad_boundaries); + const double LHS = c_light_ * (chiEVec[0] * Erad[0] - chiPVec[0] * fourPiBoverC[0]) + Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas) * (T_gas - T_d); + + if (std::abs(LHS) < lambda_rel_tol * std::abs(Lambda_compare)) { + break; + } + + const auto d_fourpib_over_c_d_t = ComputeThermalRadiationTempDerivative(T_d, rad_boundaries); + const double dLHS_dTd = - c_light_ * chiPVec[0] * d_fourpib_over_c_d_t[0] - Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas); + delta_T_d = LHS / dLHS_dTd; + T_d -= delta_T_d; + + if (ite_td > 0) { + if (std::abs(delta_T_d) < lambda_rel_tol * std::abs(T_d)) { + break; + } + } + } + + AMREX_ASSERT_WITH_MESSAGE(ite_td < max_ite_td, "Newton iteration for dust temperature did not converge."); + return T_d; +} + template void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEnergySource, amrex::Box const &indexRange, amrex::Real dt_radiation, const int stage) @@ -1229,6 +1283,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne double Etot0 = NAN; double Egas_guess = NAN; double T_gas = NAN; + double T_d = NAN; double lorentz_factor = NAN; double lorentz_factor_v = NAN; double lorentz_factor_v_v = NAN; @@ -1291,6 +1346,8 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne gas_update_factor = IMEX_a32; } + const double num_den = rho / mean_molecular_mass_; + const int max_ite = 5; int ite = 0; for (; ite < max_ite; ++ite) { @@ -1359,18 +1416,32 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne for (; n < maxIter; ++n) { T_gas = quokka::EOS::ComputeTgasFromEint(rho, Egas_guess, massScalars); AMREX_ASSERT(T_gas >= 0.); - fourPiBoverC = ComputeThermalRadiation(T_gas, radBoundaries_g_copy); + + // dust temperature + if constexpr (!use_dust) { + T_d = T_gas; + } else { + if (n == 0) { + T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, radBoundaries_g_copy); + } else { + const auto Lambda_gd = Rvec[0] / (dt * chat / c); + T_d = T_gas - Lambda_gd / (Lambda_gd_0 * num_den * num_den * std::sqrt(T_gas)); + } + } + AMREX_ASSERT(T_d >= 0.); + + fourPiBoverC = ComputeThermalRadiation(T_d, radBoundaries_g_copy); if constexpr (opacity_model_ == OpacityModel::single_group) { - kappaPVec[0] = ComputePlanckOpacity(rho, T_gas); - kappaEVec[0] = ComputeEnergyMeanOpacity(rho, T_gas); + kappaPVec[0] = ComputePlanckOpacity(rho, T_d); + kappaEVec[0] = ComputeEnergyMeanOpacity(rho, T_d); } else { - kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_d); for (int g = 0; g < nGroups_; ++g) { auto const nu_L = radBoundaries_g_copy[g]; auto const nu_R = radBoundaries_g_copy[g + 1]; - auto const B_L = PlanckFunction(nu_L, T_gas); // 4 pi B(nu) / c - auto const B_R = PlanckFunction(nu_R, T_gas); // 4 pi B(nu) / c + auto const B_L = PlanckFunction(nu_L, T_d); // 4 pi B(nu) / c + auto const B_R = PlanckFunction(nu_R, T_d); // 4 pi B(nu) / c auto const kappa_L = kappa_expo_and_lower_value[1][g]; auto const kappa_R = kappa_L * std::pow(nu_R / nu_L, kappa_expo_and_lower_value[0][g]); delta_nu_kappa_B_at_edge[g] = nu_R * kappa_R * B_R - nu_L * kappa_L * B_L; @@ -1407,13 +1478,13 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // In the first loop, calculate kappaF, work, tau0, R if (n == 0) { if constexpr (opacity_model_ == OpacityModel::single_group) { - kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_gas); + kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_d); } else { for (int g = 0; g < nGroups_; ++g) { auto const nu_L = radBoundaries_g_copy[g]; auto const nu_R = radBoundaries_g_copy[g + 1]; - auto const B_L = PlanckFunction(nu_L, T_gas); // 4 pi B(nu) / c - auto const B_R = PlanckFunction(nu_R, T_gas); // 4 pi B(nu) / c + auto const B_L = PlanckFunction(nu_L, T_d); // 4 pi B(nu) / c + auto const B_R = PlanckFunction(nu_R, T_d); // 4 pi B(nu) / c auto const kappa_L = kappa_expo_and_lower_value[1][g]; auto const kappa_R = kappa_L * std::pow(nu_R / nu_L, kappa_expo_and_lower_value[0][g]); delta_nu_kappa_B_at_edge[g] = nu_R * kappa_R * B_R - nu_L * kappa_L * B_L; @@ -1526,18 +1597,27 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne const double c_v = quokka::EOS::ComputeEintTempDerivative(rho, T_gas, massScalars); // Egas = c_v * T - const auto dfourPiB_dTgas = chat * ComputeThermalRadiationTempDerivative(T_gas, radBoundaries_g_copy); - AMREX_ASSERT(!dfourPiB_dTgas.hasnan()); + const auto d_fourpiboverc_d_t = ComputeThermalRadiationTempDerivative(T_d, radBoundaries_g_copy); + AMREX_ASSERT(!d_fourpiboverc_d_t.hasnan()); + + const auto d_Td_d_T = 3. / 2. - T_d / (2. * T_gas); // compute Jacobian elements // I assume (kappaPVec / kappaEVec) is constant here. This is usually a reasonable assumption. Note that this assumption // only affects the convergence rate of the Newton-Raphson iteration and does not affect the converged solution at all. dFG_dEgas = 1.0; + const double coeff_n = (chat / c) * dt * Lambda_gd_0 * num_den * num_den; for (int g = 0; g < nGroups_; ++g) { if (tau[g] <= 0.0) { dFR_i_dD_i[g] = -std::numeric_limits::infinity(); } else { - dFR_i_dD_i[g] = -1.0 * (1.0 / tau[g] * kappaPoverE[g] + 1.0); + double dBg_dRg = NAN; + if (use_dust) { + dBg_dRg = d_fourpiboverc_d_t[g] * (-1.0 / (coeff_n * std::sqrt(T_gas))); + } else { + dBg_dRg = 0.0; + } + dFR_i_dD_i[g] = kappaPoverE[g] * (dBg_dRg - 1.0 / tau[g]) - 1.0; } } if constexpr (use_D_as_base) { @@ -1546,7 +1626,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } else { dFG_dD.fillin(c / chat); } - dFR_dEgas = 1.0 / c_v * kappaPoverE * (dfourPiB_dTgas / chat); + dFR_dEgas = 1.0 / c_v * kappaPoverE * d_fourpiboverc_d_t * d_Td_d_T; // update variables RadSystem::SolveLinearEqs(dFG_dEgas, dFG_dD, dFR_dEgas, dFR_i_dD_i, -F_G, -1. * F_D, deltaEgas, deltaD); @@ -1567,20 +1647,20 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } // END NEWTON-RAPHSON LOOP AMREX_ALWAYS_ASSERT_WITH_MESSAGE(n < maxIter, "Newton-Raphson iteration failed to converge!"); - // std::cout << "Newton-Raphson converged after " << n << " it." << std::endl; + std::cout << "Newton-Raphson converged after " << n << " it." << std::endl; AMREX_ALWAYS_ASSERT(Egas_guess > 0.0); AMREX_ALWAYS_ASSERT(min(EradVec_guess) >= 0.0); if (n > 0) { // calculate kappaF since the temperature has changed if constexpr (opacity_model_ == OpacityModel::single_group) { - kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_gas); + kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_d); } else { for (int g = 0; g < nGroups_; ++g) { auto const nu_L = radBoundaries_g_copy[g]; auto const nu_R = radBoundaries_g_copy[g + 1]; - auto const B_L = PlanckFunction(nu_L, T_gas); // 4 pi B(nu) / c - auto const B_R = PlanckFunction(nu_R, T_gas); // 4 pi B(nu) / c + auto const B_L = PlanckFunction(nu_L, T_d); // 4 pi B(nu) / c + auto const B_R = PlanckFunction(nu_R, T_d); // 4 pi B(nu) / c auto const kappa_L = kappa_expo_and_lower_value[1][g]; auto const kappa_R = kappa_L * std::pow(nu_R / nu_L, kappa_expo_and_lower_value[0][g]); delta_nu_kappa_B_at_edge[g] = nu_R * kappa_R * B_R - nu_L * kappa_L * B_L; @@ -1603,15 +1683,16 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } } } else { // if constexpr gamma_ == 1.0 + T_d = T_gas; if constexpr (opacity_model_ == OpacityModel::single_group) { - kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_gas); + kappaFVec[0] = ComputeFluxMeanOpacity(rho, T_d); } else if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { - kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_d); for (int g = 0; g < nGroups_; ++g) { kappaFVec[g] = kappa_expo_and_lower_value[1][g]; } } else { - kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_gas); + kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_d); kappaFVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, radBoundaryRatios_copy, alpha_quant_minus_one); } } From db5badb80d3f8fa379079a30f1026debf76cf155 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 15 Aug 2024 21:11:00 +1000 Subject: [PATCH 09/25] make Rvec vector --- src/radiation/radiation_system.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 9f055598f..063091256 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -1201,6 +1201,9 @@ RadSystem::ComputeDustTemperature(double const T_gas, double const T_ chiEVec[0] = ComputeEnergyMeanOpacity(rho, T_d) * rho; fourPiBoverC = ComputeThermalRadiation(T_d, rad_boundaries); + + + const double LHS = c_light_ * (chiEVec[0] * Erad[0] - chiPVec[0] * fourPiBoverC[0]) + Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas) * (T_gas - T_d); if (std::abs(LHS) < lambda_rel_tol * std::abs(Lambda_compare)) { @@ -1422,9 +1425,9 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne T_d = T_gas; } else { if (n == 0) { - T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, radBoundaries_g_copy); + T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, radBoundaries_g_copy); // TODO } else { - const auto Lambda_gd = Rvec[0] / (dt * chat / c); + const auto Lambda_gd = sum(Rvec) / (dt * chat / c); T_d = T_gas - Lambda_gd / (Lambda_gd_0 * num_den * num_den * std::sqrt(T_gas)); } } From 870deb9244e79f79b0e8d710dd3933a62a39c2b4 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 15 Aug 2024 12:47:51 +0000 Subject: [PATCH 10/25] add Backtrace to gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 7f29f70b8..bfc387a09 100644 --- a/.gitignore +++ b/.gitignore @@ -17,6 +17,7 @@ sims/ __pycache__ *.code-workspace *.sarif +Backtrace.* *.aux *.dbl From 6af5e610295e34175939ab4b7a7766a1543f1bc8 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 15 Aug 2024 12:48:43 +0000 Subject: [PATCH 11/25] add RadDustMG --- src/problems/CMakeLists.txt | 1 + src/problems/RadDust/CMakeLists.txt | 2 +- src/problems/RadDust/test_rad_dust.hpp | 5 - src/problems/RadDustMG/CMakeLists.txt | 7 + src/problems/RadDustMG/test_rad_dust_MG.cpp | 273 ++++++++++++++++++++ src/problems/RadDustMG/test_rad_dust_MG.hpp | 21 ++ 6 files changed, 303 insertions(+), 6 deletions(-) create mode 100644 src/problems/RadDustMG/CMakeLists.txt create mode 100644 src/problems/RadDustMG/test_rad_dust_MG.cpp create mode 100644 src/problems/RadDustMG/test_rad_dust_MG.hpp diff --git a/src/problems/CMakeLists.txt b/src/problems/CMakeLists.txt index 0ca223b15..24c827be7 100644 --- a/src/problems/CMakeLists.txt +++ b/src/problems/CMakeLists.txt @@ -34,6 +34,7 @@ add_subdirectory(RadSuOlson) add_subdirectory(RadTophat) add_subdirectory(RadTube) add_subdirectory(RadDust) +add_subdirectory(RadDustMG) add_subdirectory(RadhydroShell) add_subdirectory(RadhydroShock) diff --git a/src/problems/RadDust/CMakeLists.txt b/src/problems/RadDust/CMakeLists.txt index 15b8ac221..51079803b 100644 --- a/src/problems/RadDust/CMakeLists.txt +++ b/src/problems/RadDust/CMakeLists.txt @@ -4,4 +4,4 @@ if(AMReX_GPU_BACKEND MATCHES "CUDA") setup_target_for_cuda_compilation(test_rad_dust) endif(AMReX_GPU_BACKEND MATCHES "CUDA") -add_test(NAME RadhydroUniformAdvecting COMMAND test_rad_dust RadDust.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) +add_test(NAME RadDust COMMAND test_rad_dust RadDust.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/problems/RadDust/test_rad_dust.hpp b/src/problems/RadDust/test_rad_dust.hpp index 20ddcfd67..7abf33616 100644 --- a/src/problems/RadDust/test_rad_dust.hpp +++ b/src/problems/RadDust/test_rad_dust.hpp @@ -1,10 +1,5 @@ #ifndef TEST_RAD_DUST_HPP_ // NOLINT #define TEST_RAD_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_rad_dust.hpp /// \brief Defines a test problem for radiation in the static diffusion regime. /// diff --git a/src/problems/RadDustMG/CMakeLists.txt b/src/problems/RadDustMG/CMakeLists.txt new file mode 100644 index 000000000..7154a874a --- /dev/null +++ b/src/problems/RadDustMG/CMakeLists.txt @@ -0,0 +1,7 @@ +add_executable(test_rad_dust_MG test_rad_dust_MG.cpp ../../util/fextract.cpp ../../math/interpolate.cpp ${QuokkaObjSources}) + +if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_rad_dust_MG) +endif(AMReX_GPU_BACKEND MATCHES "CUDA") + +add_test(NAME RadDustMG COMMAND test_rad_dust_MG RadDust.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp new file mode 100644 index 000000000..966909a19 --- /dev/null +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -0,0 +1,273 @@ +/// \file test_rad_dust_MG.cpp +/// \brief Defines a test problem for radiation advection in a uniform medium with grey radiation. +/// + +#include "test_rad_dust_MG.hpp" +#include "AMReX_BC_TYPES.H" +#include "AMReX_Print.H" +#include "QuokkaSimulation.hpp" +#include "fundamental_constants.H" +#include "physics_info.hpp" +#include "util/fextract.hpp" +#include + +struct DustProblem { +}; // dummy type to allow compile-type polymorphism via template specialization + +// CGS + +constexpr int beta_order_ = 1; // order of beta in the radiation four-force +constexpr double c = 1.0e8; +constexpr double chat = c; +constexpr double v0 = 0.0; +constexpr double chi0 = 10000.0; + +constexpr double T0 = 1.0; +constexpr double T_equi = 0.7680325; +constexpr double rho0 = 1.0; +constexpr double a_rad = 1.0; +constexpr double mu = 1.0; +constexpr double k_B = 1.0; + +constexpr double max_time = 1.0e-5; +constexpr double delta_time = 1.0e-8; + +constexpr double Erad0 = a_rad * T0 * T0 * T0 * T0; +constexpr double erad_floor = 1.0e-20 * Erad0; + +template <> struct SimulationData { + std::vector t_vec_; + std::vector Trad_vec_; + std::vector Tgas_vec_; +}; + +template <> struct quokka::EOS_Traits { + static constexpr double mean_molecular_weight = mu; + static constexpr double boltzmann_constant = k_B; + static constexpr double gamma = 5. / 3.; +}; + +template <> struct Physics_Traits { + // cell-centred + static constexpr bool is_hydro_enabled = true; + 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 = 2; +}; + +template <> struct RadSystem_Traits { + 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 = beta_order_; + static constexpr double energy_unit = 1.; + static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{1.0e-3, 1.0, 1.0e3}; + // 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 AMREX_FORCE_INLINE auto +RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, const double rho, + const double /*Tgas*/) -> amrex::GpuArray, 2> +{ + amrex::GpuArray, 2> exponents_and_values{}; + for (int i = 0; i < nGroups_ + 1; ++i) { + exponents_and_values[0][i] = 0.0; + exponents_and_values[1][i] = chi0 / rho; + } + return exponents_and_values; +} + +template <> +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiation(amrex::Real temperature, amrex::GpuArray const &boundaries) + -> quokka::valarray +{ + auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); + const double power = radiation_constant_ * temperature; + auto Erad_g = power * radEnergyFractions; + return Erad_g; +} + +template <> +AMREX_GPU_HOST_DEVICE auto +RadSystem::ComputeThermalRadiationTempDerivative(amrex::Real temperature, + amrex::GpuArray const &boundaries) -> quokka::valarray +{ + auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); + const double d_power_dt = radiation_constant_; + return d_power_dt * radEnergyFractions; +} + +template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) +{ + // extract variables required from the geom object + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + const double Egas = quokka::EOS::ComputeEintFromTgas(rho0, T0); + + // 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::nGroups; ++g) { + state_cc(i, j, k, RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g) = erad_floor; + state_cc(i, j, k, RadSystem::x1RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + state_cc(i, j, k, RadSystem::x2RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + state_cc(i, j, k, RadSystem::x3RadFlux_index + Physics_NumVars::numRadVars * g) = 0; + } + state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho0 * v0 * v0; + state_cc(i, j, k, RadSystem::gasDensity_index) = rho0; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Egas; + state_cc(i, j, k, RadSystem::x1GasMomentum_index) = v0 * rho0; + state_cc(i, j, k, RadSystem::x2GasMomentum_index) = 0.; + state_cc(i, j, k, RadSystem::x3GasMomentum_index) = 0.; + }); +} + +template <> void QuokkaSimulation::computeAfterTimestep() +{ + auto [position, values] = fextract(state_new_cc_[0], Geom(0), 0, 0.5); + + if (amrex::ParallelDescriptor::IOProcessor()) { + 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 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]; + const amrex::Real Egas_i = RadSystem::ComputeEintFromEgas(rho, x1GasMom, x2GasMom, x3GasMom, Etot_i); + const amrex::Real Erad_i = values.at(RadSystem::radEnergy_index)[0]; + // userData_.Trad_vec_.push_back(std::pow(Erad_i / a_rad, 1. / 4.)); + userData_.Trad_vec_.push_back(Erad_i / a_rad); + userData_.Tgas_vec_.push_back(quokka::EOS::ComputeTgasFromEint(rho, Egas_i)); + } +} + +auto problem_main() -> int +{ + // Problem parameters + const int max_timesteps = 1e6; + const double CFL_number_gas = 0.8; + const double CFL_number_rad = 8.0; + + // Boundary conditions + constexpr int nvars = RadSystem::nvar_; + amrex::Vector BCs_cc(nvars); + for (int n = 0; n < nvars; ++n) { + for (int i = 0; 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 sim(BCs_cc); + + sim.radiationReconstructionOrder_ = 3; // PPM + sim.stopTime_ = max_time; + sim.radiationCflNumber_ = CFL_number_rad; + sim.cflNumber_ = CFL_number_gas; + sim.maxTimesteps_ = max_timesteps; + sim.plotfileInterval_ = -1; + sim.initDt_ = delta_time; + sim.maxDt_ = delta_time; + + // initialize + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + // read in exact solution + std::vector ts_exact{}; + std::vector Trad_exact{}; + std::vector Tgas_exact{}; + + std::ifstream fstream("../extern/data/dust/rad_dust_exact.csv", std::ios::in); + AMREX_ALWAYS_ASSERT(fstream.is_open()); + std::string header; + std::getline(fstream, header); + + for (std::string line; std::getline(fstream, line);) { + std::istringstream iss(line); + std::vector values; + std::string value; + + while (std::getline(iss, value, ',')) { + values.push_back(std::stod(value)); + } + auto t_val = values.at(0); + auto Tmat_val = values.at(1); + auto Trad_val = values.at(2); + if (t_val <= 0.0) { + continue; + } + ts_exact.push_back(t_val); + Tgas_exact.push_back(Tmat_val); + Trad_exact.push_back(Trad_val); + } + + std::vector &Tgas = sim.userData_.Tgas_vec_; + std::vector &Trad = sim.userData_.Trad_vec_; + std::vector &t = sim.userData_.t_vec_; + + std::vector Tgas_interp(t.size()); + std::vector Trad_interp(t.size()); + interpolate_arrays(t.data(), Tgas_interp.data(), static_cast(t.size()), ts_exact.data(), Tgas_exact.data(), static_cast(ts_exact.size())); + interpolate_arrays(t.data(), Trad_interp.data(), static_cast(t.size()), ts_exact.data(), Trad_exact.data(), static_cast(ts_exact.size())); + + // compute error norm + double err_norm = 0.; + double sol_norm = 0.; + for (size_t i = 0; i < t.size(); ++i) { + err_norm += std::abs(Tgas[i] - Tgas_interp[i]); + err_norm += std::abs(Trad[i] - Trad_interp[i]); + sol_norm += std::abs(Tgas_interp[i]) + std::abs(Trad_interp[i]); + } + const double error_tol = 0.0008; + const double rel_error = err_norm / sol_norm; + amrex::Print() << "Relative L1 error norm = " << rel_error << std::endl; + +#ifdef HAVE_PYTHON + // plot temperature + matplotlibcpp::clf(); + matplotlibcpp::xscale("log"); + std::map Trad_args; + std::map Tgas_args; + std::map Texact_args; + std::map Tradexact_args; + Trad_args["label"] = "radiation (numerical)"; + Trad_args["linestyle"] = "--"; + Trad_args["color"] = "C1"; + Tradexact_args["label"] = "radiation (exact)"; + Tradexact_args["linestyle"] = "-"; + Tradexact_args["color"] = "k"; + Tgas_args["label"] = "gas (numerical)"; + Tgas_args["linestyle"] = "--"; + Tgas_args["color"] = "C2"; + Texact_args["label"] = "gas (exact)"; + Texact_args["linestyle"] = "-"; + Texact_args["color"] = "k"; + matplotlibcpp::plot(ts_exact, Tgas_exact, Texact_args); + matplotlibcpp::plot(ts_exact, Trad_exact, Tradexact_args); + matplotlibcpp::plot(t, Tgas, Tgas_args); + matplotlibcpp::plot(t, Trad, Trad_args); + matplotlibcpp::xlabel("t (dimensionless)"); + matplotlibcpp::ylabel("T (dimensionless)"); + matplotlibcpp::legend(); + matplotlibcpp::tight_layout(); + matplotlibcpp::save("./rad_dust_MG_T.pdf"); +#endif + + // Cleanup and exit + int status = 0; + if ((rel_error > error_tol) || std::isnan(rel_error)) { + status = 1; + } + return status; +} diff --git a/src/problems/RadDustMG/test_rad_dust_MG.hpp b/src/problems/RadDustMG/test_rad_dust_MG.hpp new file mode 100644 index 000000000..6187b8b7f --- /dev/null +++ b/src/problems/RadDustMG/test_rad_dust_MG.hpp @@ -0,0 +1,21 @@ +#ifndef TEST_RAD_DUST_MG_HPP_ // NOLINT +#define TEST_RAD_DUST_MG_HPP_ +/// \file test_rad_dust_MG.hpp +/// \brief Defines a test problem for radiation in the static diffusion regime. +/// + +// external headers +#ifdef HAVE_PYTHON +#include "util/matplotlibcpp.h" +#endif +#include +#include + +// internal headers + +#include "math/interpolate.hpp" +#include "radiation/radiation_system.hpp" + +// function definitions + +#endif // TEST_RAD_DUST_MG_HPP_ From 7a81c19ea2544d916fa559e3caf235ca5f14b6ca Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Thu, 15 Aug 2024 12:50:06 +0000 Subject: [PATCH 12/25] make T_d work in MG case; got negative T_d error --- src/radiation/radiation_system.hpp | 103 +++++++++++++++++++---------- 1 file changed, 68 insertions(+), 35 deletions(-) diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 063091256..fdab8b305 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -268,7 +268,7 @@ template class RadSystem : public HyperbolicSystem const &Erad, - amrex::GpuArray const &rad_boundaries) -> double; + amrex::GpuArray const &rad_boundaries, amrex::GpuArray const &rad_boundary_ratios) -> double; template AMREX_GPU_DEVICE static auto @@ -287,27 +287,45 @@ template class RadSystem : public HyperbolicSystem AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckEnergyFractions(amrex::GpuArray const &boundaries, amrex::Real temperature) -> quokka::valarray { quokka::valarray radEnergyFractions{}; if constexpr (nGroups_ == 1) { - // TODO(CCH): allow the total radEnergyFraction to be smaller than 1. One usage case is to allow, say, a single group representing IR radiation. radEnergyFractions[0] = 1.0; return radEnergyFractions; } else { amrex::Real const energy_unit_over_kT = RadSystem_Traits::energy_unit / (boltzmann_constant_ * temperature); - amrex::Real previous = integrate_planck_from_0_to_x(boundaries[0] * energy_unit_over_kT); - for (int g = 0; g < nGroups_; ++g) { - amrex::Real y = integrate_planck_from_0_to_x(boundaries[g + 1] * energy_unit_over_kT); + // amrex::Real x = boundaries[0] * energy_unit_over_kT; + // amrex::Real previous = NAN; + // if (x >= 100.) { + // previous = 1.0; + // } else { + // previous = integrate_planck_from_0_to_x(x); + // } + amrex::Real y = NAN; + amrex::Real previous = 0.0; + for (int g = 0; g < nGroups_ - 1; ++g) { + const amrex::Real x = boundaries[g + 1] * energy_unit_over_kT; + if (x >= 100.) { + y = 1.0; + } else { + y = integrate_planck_from_0_to_x(x); + } radEnergyFractions[g] = y - previous; previous = y; } - auto tote = sum(radEnergyFractions); + // last group + y = 1.0; + radEnergyFractions[nGroups_ - 1] = y - previous; + + const auto tote = sum(radEnergyFractions); + AMREX_ALWAYS_ASSERT(std::abs(tote - 1.0) < 1.0e-10); // AMREX_ALWAYS_ASSERT(tote <= 1.0); // AMREX_ALWAYS_ASSERT(tote > 0.9999); - radEnergyFractions /= tote; + // radEnergyFractions /= tote; return radEnergyFractions; } } @@ -1138,10 +1156,13 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeDiffusionFluxMeanOpacity // delta_nu_kappa_B_at_edge[g]; kappaF[g] = (kappaPVec[g] + 1. / 3. * kappaEVec[g]) * fourPiBoverC[g] + 1. / 3. * (kappa_slope[g] * kappaEVec[g] * fourPiBoverC[g] - delta_nu_kappa_B_at_edge[g]); - AMREX_ASSERT(kappaF[g] > 0.0); auto const denom = 4. / 3. * fourPiBoverC[g] - 1. / 3. * delta_nu_B_at_edge[g]; - AMREX_ASSERT(denom > 0.0); - kappaF[g] /= denom; + if (denom <= 0.0) { + AMREX_ASSERT(kappaF[g] == 0.0); + kappaF[g] = 0.0; + } else { + kappaF[g] /= denom; + } } return kappaF; } @@ -1181,38 +1202,60 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxInDiffusionLimit(con template AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeDustTemperature(double const T_gas, double const T_d_init, double const rho, quokka::valarray const &Erad, - amrex::GpuArray const &rad_boundaries) -> double + amrex::GpuArray const &rad_boundaries, amrex::GpuArray const &rad_boundary_ratios) -> double { - quokka::valarray fourPiBoverC{}; - quokka::valarray chiPVec{}; - quokka::valarray chiEVec{}; + quokka::valarray kappaPVec{}; + quokka::valarray kappaEVec{}; const double num_density = rho / mean_molecular_mass_; const double Lambda_compare = Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas) * T_gas; - // solve for dust temperature: Newton iteration on the dust temperature T_d + amrex::GpuArray alpha_quant_minus_one{}; + for (int g = 0; g < nGroups_; ++g) { + alpha_quant_minus_one[g] = -1.0; + } + + // solve for dust temperature T_d using Newton iteration double T_d = T_d_init; const double lambda_rel_tol = 1.0e-8; const int max_ite_td = 100; int ite_td = 0; - double delta_T_d = NAN; for (; ite_td < max_ite_td; ++ite_td) { - chiPVec[0] = ComputePlanckOpacity(rho, T_d) * rho; - chiEVec[0] = ComputeEnergyMeanOpacity(rho, T_d) * rho; - - fourPiBoverC = ComputeThermalRadiation(T_d, rad_boundaries); - + const auto fourPiBoverC = ComputeThermalRadiation(T_d, rad_boundaries); + if constexpr (opacity_model_ == OpacityModel::single_group) { + kappaPVec[0] = ComputePlanckOpacity(rho, T_d); + kappaEVec[0] = ComputeEnergyMeanOpacity(rho, T_d); + } else { + const auto kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(rad_boundaries, rho, T_d); + if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { + for (int g = 0; g < nGroups_; ++g) { + kappaPVec[g] = kappa_expo_and_lower_value[1][g]; + kappaEVec[g] = kappa_expo_and_lower_value[1][g]; + } + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { + kappaPVec = + ComputeGroupMeanOpacity(kappa_expo_and_lower_value, rad_boundary_ratios, alpha_quant_minus_one); + kappaEVec = kappaPVec; + } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { + const auto alpha_B = ComputeRadQuantityExponents(fourPiBoverC, rad_boundaries); + const auto alpha_E = ComputeRadQuantityExponents(Erad, rad_boundaries); + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, rad_boundary_ratios, alpha_B); + kappaEVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, rad_boundary_ratios, alpha_E); + } + } + AMREX_ASSERT(!kappaPVec.hasnan()); + AMREX_ASSERT(!kappaEVec.hasnan()); - const double LHS = c_light_ * (chiEVec[0] * Erad[0] - chiPVec[0] * fourPiBoverC[0]) + Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas) * (T_gas - T_d); + const double LHS = c_light_ * rho * sum(kappaEVec * Erad - kappaPVec * fourPiBoverC) + Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas) * (T_gas - T_d); if (std::abs(LHS) < lambda_rel_tol * std::abs(Lambda_compare)) { break; } const auto d_fourpib_over_c_d_t = ComputeThermalRadiationTempDerivative(T_d, rad_boundaries); - const double dLHS_dTd = - c_light_ * chiPVec[0] * d_fourpib_over_c_d_t[0] - Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas); - delta_T_d = LHS / dLHS_dTd; + const double dLHS_dTd = - c_light_ * rho * sum(kappaPVec * d_fourpib_over_c_d_t) - Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas); + const double delta_T_d = LHS / dLHS_dTd; T_d -= delta_T_d; if (ite_td > 0) { @@ -1425,7 +1468,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne T_d = T_gas; } else { if (n == 0) { - T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, radBoundaries_g_copy); // TODO + T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, radBoundaries_g_copy, radBoundaryRatios_copy); } else { const auto Lambda_gd = sum(Rvec) / (dt * chat / c); T_d = T_gas - Lambda_gd / (Lambda_gd_0 * num_den * num_den * std::sqrt(T_gas)); @@ -1440,16 +1483,6 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne kappaEVec[0] = ComputeEnergyMeanOpacity(rho, T_d); } else { kappa_expo_and_lower_value = DefineOpacityExponentsAndLowerValues(radBoundaries_g_copy, rho, T_d); - for (int g = 0; g < nGroups_; ++g) { - auto const nu_L = radBoundaries_g_copy[g]; - auto const nu_R = radBoundaries_g_copy[g + 1]; - auto const B_L = PlanckFunction(nu_L, T_d); // 4 pi B(nu) / c - auto const B_R = PlanckFunction(nu_R, T_d); // 4 pi B(nu) / c - auto const kappa_L = kappa_expo_and_lower_value[1][g]; - auto const kappa_R = kappa_L * std::pow(nu_R / nu_L, kappa_expo_and_lower_value[0][g]); - delta_nu_kappa_B_at_edge[g] = nu_R * kappa_R * B_R - nu_L * kappa_L * B_L; - delta_nu_B_at_edge[g] = nu_R * B_R - nu_L * B_L; - } if constexpr (opacity_model_ == OpacityModel::piecewise_constant_opacity) { for (int g = 0; g < nGroups_; ++g) { kappaPVec[g] = kappa_expo_and_lower_value[1][g]; From cd9f4ee6bd40ca8488cf35cd2bbe07748cacb423 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 16 Aug 2024 12:17:04 +1000 Subject: [PATCH 13/25] dust_MG working!!! --- src/problems/RadDustMG/test_rad_dust_MG.cpp | 5 +- src/radiation/radiation_system.hpp | 59 +++++++++++++-------- 2 files changed, 40 insertions(+), 24 deletions(-) diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp index 966909a19..8f07a00ab 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.cpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -141,7 +141,10 @@ template <> void QuokkaSimulation::computeAfterTimestep() const amrex::Real x3GasMom = values.at(RadSystem::x3GasMomentum_index)[0]; const amrex::Real rho = values.at(RadSystem::gasDensity_index)[0]; const amrex::Real Egas_i = RadSystem::ComputeEintFromEgas(rho, x1GasMom, x2GasMom, x3GasMom, Etot_i); - const amrex::Real Erad_i = values.at(RadSystem::radEnergy_index)[0]; + double Erad_i = 0.0; + for (int g = 0; g < Physics_Traits::nGroups; ++g) { + Erad_i += values.at(RadSystem::radEnergy_index + Physics_NumVars::numRadVars * g)[0]; + } // userData_.Trad_vec_.push_back(std::pow(Erad_i / a_rad, 1. / 4.)); userData_.Trad_vec_.push_back(Erad_i / a_rad); userData_.Tgas_vec_.push_back(quokka::EOS::ComputeTgasFromEint(rho, Egas_i)); diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index fdab8b305..0876bf10b 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -1448,11 +1448,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // dF_{D,i} / dD_i = - (1 / (chat * dt * rho * kappa_{E,i}) + 1) * tau0_i = - ((1 / tau_i)(kappa_Pi / kappa_Ei) + 1) * tau0_i double F_G = NAN; - double dFG_dEgas = NAN; double deltaEgas = NAN; - quokka::valarray dFG_dD{}; - quokka::valarray dFR_dEgas{}; - quokka::valarray dFR_i_dD_i{}; quokka::valarray deltaD{}; quokka::valarray F_D{}; @@ -1636,42 +1632,59 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne const auto d_fourpiboverc_d_t = ComputeThermalRadiationTempDerivative(T_d, radBoundaries_g_copy); AMREX_ASSERT(!d_fourpiboverc_d_t.hasnan()); - const auto d_Td_d_T = 3. / 2. - T_d / (2. * T_gas); - // compute Jacobian elements // I assume (kappaPVec / kappaEVec) is constant here. This is usually a reasonable assumption. Note that this assumption // only affects the convergence rate of the Newton-Raphson iteration and does not affect the converged solution at all. - dFG_dEgas = 1.0; - const double coeff_n = (chat / c) * dt * Lambda_gd_0 * num_den * num_den; + + const double cscale = c / chat; + auto dEg_dT = kappaPoverE * d_fourpiboverc_d_t; + + const double y0 = -F_G; + auto yg = -1. * F_D; + + quokka::valarray dF0_dXg{}; + quokka::valarray dFg_dX0{}; + quokka::valarray dFg_dXg{}; + + // M_00 + const double dF0_dX0 = 1.0; + // M_0g + dF0_dXg.fillin(cscale); + // M_g0 + if constexpr (!use_dust) { + dFg_dX0 = 1.0 / c_v * dEg_dT; + } else { + const double d_Td_d_T = 3. / 2. - T_d / (2. * T_gas); + const double coeff_n = dt * Lambda_gd_0 * num_den * num_den / cscale; + 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; + dFg_dX0 = 1.0 / c_v * dEg_dT - 1.0 / cscale * rg * dF0_dX0; + yg = yg - 1.0 / cscale * rg * y0; + } + // M_gg, same for dust and dust-free cases for (int g = 0; g < nGroups_; ++g) { if (tau[g] <= 0.0) { - dFR_i_dD_i[g] = -std::numeric_limits::infinity(); + dFg_dXg[g] = -std::numeric_limits::infinity(); } else { - double dBg_dRg = NAN; - if (use_dust) { - dBg_dRg = d_fourpiboverc_d_t[g] * (-1.0 / (coeff_n * std::sqrt(T_gas))); - } else { - dBg_dRg = 0.0; - } - dFR_i_dD_i[g] = kappaPoverE[g] * (dBg_dRg - 1.0 / tau[g]) - 1.0; + dFg_dXg[g] = -1.0 * kappaPoverE[g] / tau[g] - 1.0; } } + if constexpr (use_D_as_base) { - dFG_dD = (c / chat) * tau0; - dFR_i_dD_i = dFR_i_dD_i * tau0; - } else { - dFG_dD.fillin(c / chat); + dF0_dXg = dF0_dXg * tau0; + dFg_dXg = dFg_dXg * tau0; } - dFR_dEgas = 1.0 / c_v * kappaPoverE * d_fourpiboverc_d_t * d_Td_d_T; // update variables - RadSystem::SolveLinearEqs(dFG_dEgas, dFG_dD, dFR_dEgas, dFR_i_dD_i, -F_G, -1. * F_D, deltaEgas, deltaD); + RadSystem::SolveLinearEqs(dF0_dX0, dF0_dXg, dFg_dX0, dFg_dXg, y0, yg, deltaEgas, deltaD); AMREX_ASSERT(!std::isnan(deltaEgas)); AMREX_ASSERT(!deltaD.hasnan()); Egas_guess += deltaEgas; if constexpr (use_D_as_base) { - D += deltaD; + D += deltaD; // TODO: remove D and replace with Rvec += deltaD * tau0 + // Rvec = tau0 * D; } else { Rvec += deltaD; } From 125124f2856535325f278bc2d01ad5e566aa95d7 Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 16 Aug 2024 12:19:42 +1000 Subject: [PATCH 14/25] removed D --- src/radiation/radiation_system.hpp | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 0876bf10b..d008f0f0f 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -1344,7 +1344,6 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne quokka::valarray kappaPoverE{}; quokka::valarray tau0{}; // optical depth across c * dt at old state quokka::valarray tau{}; // optical depth across c * dt at new state - quokka::valarray D{}; // D = S / tau0 quokka::valarray work{}; quokka::valarray work_prev{}; amrex::GpuArray dMomentum{}; @@ -1584,13 +1583,9 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne tau0[g] = 1.0; } } - D = Rvec / tau0; } } else { // in the second and later loops, calculate tau and E (given R) tau = dt * rho * kappaPVec * chat * lorentz_factor; - if constexpr (use_D_as_base) { - Rvec = tau0 * D; - } for (int g = 0; g < nGroups_; ++g) { // If tau = 0.0, Erad_guess shouldn't change if (tau[g] > 0.0) { @@ -1683,8 +1678,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne Egas_guess += deltaEgas; if constexpr (use_D_as_base) { - D += deltaD; // TODO: remove D and replace with Rvec += deltaD * tau0 - // Rvec = tau0 * D; + Rvec += tau0 * deltaD; } else { Rvec += deltaD; } From e95c3a69a6fce6b772612d9cf8bb1f9079db325a Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 16 Aug 2024 12:24:03 +1000 Subject: [PATCH 15/25] cleanup --- src/problems/RadDust/test_rad_dust.cpp | 17 +++-------------- src/problems/RadDustMG/test_rad_dust_MG.cpp | 8 +++----- 2 files changed, 6 insertions(+), 19 deletions(-) diff --git a/src/problems/RadDust/test_rad_dust.cpp b/src/problems/RadDust/test_rad_dust.cpp index c3dad98c5..0f09fd373 100644 --- a/src/problems/RadDust/test_rad_dust.cpp +++ b/src/problems/RadDust/test_rad_dust.cpp @@ -1,5 +1,5 @@ /// \file test_rad_dust.cpp -/// \brief Defines a test problem for radiation advection in a uniform medium with grey radiation. +/// \brief Defines a single-group test problem for gas-dust-radiation coupling in uniform medium. /// #include "test_rad_dust.hpp" @@ -14,8 +14,6 @@ struct DustProblem { }; // dummy type to allow compile-type polymorphism via template specialization -// CGS - constexpr int beta_order_ = 1; // order of beta in the radiation four-force constexpr double c = 1.0e8; constexpr double chat = c; @@ -104,19 +102,10 @@ template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokk const double Egas = quokka::EOS::ComputeEintFromTgas(rho0, T0); - double erad = NAN; - double frad = NAN; - erad = Erad0; - frad = 0.0; - - // test - erad = erad_floor; - frad = 0.0; - // loop over the grid and set the initial condition amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - state_cc(i, j, k, RadSystem::radEnergy_index) = erad; - state_cc(i, j, k, RadSystem::x1RadFlux_index) = frad; + state_cc(i, j, k, RadSystem::radEnergy_index) = erad_floor; + state_cc(i, j, k, RadSystem::x1RadFlux_index) = 0; state_cc(i, j, k, RadSystem::x2RadFlux_index) = 0; state_cc(i, j, k, RadSystem::x3RadFlux_index) = 0; state_cc(i, j, k, RadSystem::gasEnergy_index) = Egas + 0.5 * rho0 * v0 * v0; diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp index 8f07a00ab..117afb14c 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.cpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -1,5 +1,5 @@ /// \file test_rad_dust_MG.cpp -/// \brief Defines a test problem for radiation advection in a uniform medium with grey radiation. +/// \brief Defines a multigroup test problem for gas-dust-radiation coupling in uniform medium. /// #include "test_rad_dust_MG.hpp" @@ -14,8 +14,6 @@ struct DustProblem { }; // dummy type to allow compile-type polymorphism via template specialization -// CGS - constexpr int beta_order_ = 1; // order of beta in the radiation four-force constexpr double c = 1.0e8; constexpr double chat = c; @@ -55,7 +53,7 @@ template <> struct Physics_Traits { static constexpr bool is_radiation_enabled = true; // face-centred static constexpr bool is_mhd_enabled = false; - static constexpr int nGroups = 2; + static constexpr int nGroups = 4; }; template <> struct RadSystem_Traits { @@ -65,7 +63,7 @@ template <> struct RadSystem_Traits { static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; static constexpr double energy_unit = 1.; - static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{1.0e-3, 1.0, 1.0e3}; + static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{1.0e-3, 0.1, 1.0, 10.0, 1.0e3}; // static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum; }; From 49a16bc25c041e3faf9e620b3495f5a7bf1ecbfa Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 16 Aug 2024 12:32:15 +1000 Subject: [PATCH 16/25] cleanup radiation_system.hpp --- src/radiation/radiation_system.hpp | 17 +++-------------- 1 file changed, 3 insertions(+), 14 deletions(-) diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index d008f0f0f..5ee692991 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -298,18 +298,11 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckEnergyFractions(am return radEnergyFractions; } else { amrex::Real const energy_unit_over_kT = RadSystem_Traits::energy_unit / (boltzmann_constant_ * temperature); - // amrex::Real x = boundaries[0] * energy_unit_over_kT; - // amrex::Real previous = NAN; - // if (x >= 100.) { - // previous = 1.0; - // } else { - // previous = integrate_planck_from_0_to_x(x); - // } amrex::Real y = NAN; amrex::Real previous = 0.0; for (int g = 0; g < nGroups_ - 1; ++g) { const amrex::Real x = boundaries[g + 1] * energy_unit_over_kT; - if (x >= 100.) { + if (x >= 100.) { // 100. is the upper limit of x in the table y = 1.0; } else { y = integrate_planck_from_0_to_x(x); @@ -317,15 +310,11 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckEnergyFractions(am radEnergyFractions[g] = y - previous; previous = y; } - // last group + // last group, enforcing the total fraction to be 1.0 y = 1.0; radEnergyFractions[nGroups_ - 1] = y - previous; + AMREX_ASSERT(std::abs(sum(radEnergyFractions) - 1.0) < 1.0e-10); - const auto tote = sum(radEnergyFractions); - AMREX_ALWAYS_ASSERT(std::abs(tote - 1.0) < 1.0e-10); - // AMREX_ALWAYS_ASSERT(tote <= 1.0); - // AMREX_ALWAYS_ASSERT(tote > 0.9999); - // radEnergyFractions /= tote; return radEnergyFractions; } } From 7599b4d88ff1469084cfa711d379d1958f98786c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 16 Aug 2024 02:39:27 +0000 Subject: [PATCH 17/25] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/problems/RadDust/test_rad_dust.cpp | 7 +++--- src/problems/RadDustMG/test_rad_dust_MG.cpp | 9 ++++---- src/radiation/radiation_system.hpp | 24 +++++++++++---------- 3 files changed, 20 insertions(+), 20 deletions(-) diff --git a/src/problems/RadDust/test_rad_dust.cpp b/src/problems/RadDust/test_rad_dust.cpp index 0f09fd373..03df3d7d7 100644 --- a/src/problems/RadDust/test_rad_dust.cpp +++ b/src/problems/RadDust/test_rad_dust.cpp @@ -1,5 +1,5 @@ /// \file test_rad_dust.cpp -/// \brief Defines a single-group test problem for gas-dust-radiation coupling in uniform medium. +/// \brief Defines a single-group test problem for gas-dust-radiation coupling in uniform medium. /// #include "test_rad_dust.hpp" @@ -85,9 +85,8 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiation(amrex } template <> -AMREX_GPU_HOST_DEVICE auto -RadSystem::ComputeThermalRadiationTempDerivative(amrex::Real temperature, - amrex::GpuArray const &boundaries) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiationTempDerivative( + amrex::Real temperature, amrex::GpuArray const &boundaries) -> quokka::valarray { auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); const double d_power_dt = radiation_constant_; diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp index 117afb14c..306c84941 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.cpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -1,5 +1,5 @@ /// \file test_rad_dust_MG.cpp -/// \brief Defines a multigroup test problem for gas-dust-radiation coupling in uniform medium. +/// \brief Defines a multigroup test problem for gas-dust-radiation coupling in uniform medium. /// #include "test_rad_dust_MG.hpp" @@ -71,7 +71,7 @@ template <> struct RadSystem_Traits { template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, const double rho, - const double /*Tgas*/) -> amrex::GpuArray, 2> + const double /*Tgas*/) -> amrex::GpuArray, 2> { amrex::GpuArray, 2> exponents_and_values{}; for (int i = 0; i < nGroups_ + 1; ++i) { @@ -92,9 +92,8 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiation(amrex } template <> -AMREX_GPU_HOST_DEVICE auto -RadSystem::ComputeThermalRadiationTempDerivative(amrex::Real temperature, - amrex::GpuArray const &boundaries) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiationTempDerivative( + amrex::Real temperature, amrex::GpuArray const &boundaries) -> quokka::valarray { auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); const double d_power_dt = radiation_constant_; diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 5ee692991..bb76ff62f 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -266,9 +266,9 @@ template class RadSystem : public HyperbolicSystem const &boundaries) -> quokka::valarray; - AMREX_GPU_HOST_DEVICE static auto - ComputeDustTemperature(double T_gas, double T_d_init, double rho, quokka::valarray const &Erad, - amrex::GpuArray const &rad_boundaries, amrex::GpuArray const &rad_boundary_ratios) -> double; + AMREX_GPU_HOST_DEVICE static auto ComputeDustTemperature(double T_gas, double T_d_init, double rho, quokka::valarray const &Erad, + amrex::GpuArray const &rad_boundaries, + amrex::GpuArray const &rad_boundary_ratios) -> double; template AMREX_GPU_DEVICE static auto @@ -1189,9 +1189,10 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxInDiffusionLimit(con } template -AMREX_GPU_HOST_DEVICE auto -RadSystem::ComputeDustTemperature(double const T_gas, double const T_d_init, double const rho, quokka::valarray const &Erad, - amrex::GpuArray const &rad_boundaries, amrex::GpuArray const &rad_boundary_ratios) -> double +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeDustTemperature(double const T_gas, double const T_d_init, double const rho, + quokka::valarray const &Erad, + amrex::GpuArray const &rad_boundaries, + amrex::GpuArray const &rad_boundary_ratios) -> double { quokka::valarray kappaPVec{}; quokka::valarray kappaEVec{}; @@ -1223,8 +1224,7 @@ RadSystem::ComputeDustTemperature(double const T_gas, double const T_ kappaEVec[g] = kappa_expo_and_lower_value[1][g]; } } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_fixed_slope_spectrum) { - kappaPVec = - ComputeGroupMeanOpacity(kappa_expo_and_lower_value, rad_boundary_ratios, alpha_quant_minus_one); + kappaPVec = ComputeGroupMeanOpacity(kappa_expo_and_lower_value, rad_boundary_ratios, alpha_quant_minus_one); kappaEVec = kappaPVec; } else if constexpr (opacity_model_ == OpacityModel::PPL_opacity_full_spectrum) { const auto alpha_B = ComputeRadQuantityExponents(fourPiBoverC, rad_boundaries); @@ -1236,14 +1236,15 @@ RadSystem::ComputeDustTemperature(double const T_gas, double const T_ AMREX_ASSERT(!kappaPVec.hasnan()); AMREX_ASSERT(!kappaEVec.hasnan()); - const double LHS = c_light_ * rho * sum(kappaEVec * Erad - kappaPVec * fourPiBoverC) + Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas) * (T_gas - T_d); + const double LHS = c_light_ * rho * sum(kappaEVec * Erad - kappaPVec * fourPiBoverC) + + Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas) * (T_gas - T_d); if (std::abs(LHS) < lambda_rel_tol * std::abs(Lambda_compare)) { break; } const auto d_fourpib_over_c_d_t = ComputeThermalRadiationTempDerivative(T_d, rad_boundaries); - const double dLHS_dTd = - c_light_ * rho * sum(kappaPVec * d_fourpib_over_c_d_t) - Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas); + const double dLHS_dTd = -c_light_ * rho * sum(kappaPVec * d_fourpib_over_c_d_t) - Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas); const double delta_T_d = LHS / dLHS_dTd; T_d -= delta_T_d; @@ -1452,7 +1453,8 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne T_d = T_gas; } else { if (n == 0) { - T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, radBoundaries_g_copy, radBoundaryRatios_copy); + T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, radBoundaries_g_copy, + radBoundaryRatios_copy); } else { const auto Lambda_gd = sum(Rvec) / (dt * chat / c); T_d = T_gas - Lambda_gd / (Lambda_gd_0 * num_den * num_den * std::sqrt(T_gas)); From 2bae32940c2f385f931568616e51be3b182c27ea Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 16 Aug 2024 13:25:07 +1000 Subject: [PATCH 18/25] add "enable_dust" constexpr parameter --- src/problems/RadBeam/test_radiation_beam.cpp | 1 + src/problems/RadDust/test_rad_dust.cpp | 1 + src/problems/RadDustMG/test_rad_dust_MG.cpp | 1 + src/problems/RadForce/test_radiation_force.cpp | 1 + src/problems/RadMarshak/test_radiation_marshak.cpp | 1 + .../test_radiation_marshak_asymptotic.cpp | 1 + .../RadMarshakCGS/test_radiation_marshak_cgs.cpp | 1 + .../RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp | 1 + .../test_radiation_matter_coupling.cpp | 1 + .../test_radiation_matter_coupling_rsla.cpp | 1 + src/problems/RadPulse/test_radiation_pulse.cpp | 1 + src/problems/RadShadow/test_radiation_shadow.cpp | 1 + src/problems/RadStreaming/test_radiation_streaming.cpp | 1 + .../RadStreamingY/test_radiation_streaming_y.cpp | 1 + src/problems/RadSuOlson/test_radiation_SuOlson.cpp | 1 + src/problems/RadTophat/test_radiation_tophat.cpp | 1 + src/problems/RadTube/test_radiation_tube.cpp | 1 + src/problems/RadhydroBB/test_radhydro_bb.cpp | 1 + src/problems/RadhydroPulse/test_radhydro_pulse.cpp | 2 ++ .../RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp | 2 ++ .../RadhydroPulseGrey/test_radhydro_pulse_grey.cpp | 2 ++ .../test_radhydro_pulse_MG_const_kappa.cpp | 2 ++ .../RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp | 2 ++ src/problems/RadhydroShell/test_radhydro_shell.cpp | 1 + src/problems/RadhydroShock/test_radhydro_shock.cpp | 1 + .../RadhydroShockCGS/test_radhydro_shock_cgs.cpp | 1 + .../test_radhydro_shock_multigroup.cpp | 1 + .../test_radhydro_uniform_advecting.cpp | 1 + src/radiation/radiation_system.hpp | 10 ++++++---- 29 files changed, 39 insertions(+), 4 deletions(-) diff --git a/src/problems/RadBeam/test_radiation_beam.cpp b/src/problems/RadBeam/test_radiation_beam.cpp index df9cc02a8..53947dc20 100644 --- a/src/problems/RadBeam/test_radiation_beam.cpp +++ b/src/problems/RadBeam/test_radiation_beam.cpp @@ -38,6 +38,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadDust/test_rad_dust.cpp b/src/problems/RadDust/test_rad_dust.cpp index 0f09fd373..7b95e1f2c 100644 --- a/src/problems/RadDust/test_rad_dust.cpp +++ b/src/problems/RadDust/test_rad_dust.cpp @@ -51,6 +51,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; + static constexpr bool enable_dust = true; }; template <> struct Physics_Traits { diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp index 117afb14c..081dc57e8 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.cpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -66,6 +66,7 @@ template <> struct RadSystem_Traits { static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{1.0e-3, 0.1, 1.0, 10.0, 1.0e3}; // static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum; + static constexpr bool enable_dust = true; }; template <> diff --git a/src/problems/RadForce/test_radiation_force.cpp b/src/problems/RadForce/test_radiation_force.cpp index d5934ebb4..fada3df80 100644 --- a/src/problems/RadForce/test_radiation_force.cpp +++ b/src/problems/RadForce/test_radiation_force.cpp @@ -71,6 +71,7 @@ template <> struct RadSystem_Traits { static constexpr double Erad_floor = 0.; static constexpr double energy_unit = C::ev2erg; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real { return 0.; } diff --git a/src/problems/RadMarshak/test_radiation_marshak.cpp b/src/problems/RadMarshak/test_radiation_marshak.cpp index 007af7b89..910f7e38a 100644 --- a/src/problems/RadMarshak/test_radiation_marshak.cpp +++ b/src/problems/RadMarshak/test_radiation_marshak.cpp @@ -42,6 +42,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 0; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp b/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp index e32e29abd..80cb7c386 100644 --- a/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp +++ b/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp @@ -39,6 +39,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = Erad_floor_; static constexpr int beta_order = 0; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp b/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp index aed8d04bb..fa9e08206 100644 --- a/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp +++ b/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp @@ -44,6 +44,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp b/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp index d9d8027b4..a4c3737ac 100644 --- a/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp +++ b/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp @@ -124,6 +124,7 @@ template <> struct RadSystem_Traits { static constexpr double energy_unit = C::hplanck; // set boundary unit to Hz static constexpr amrex::GpuArray radBoundaries = group_edges_; static constexpr OpacityModel opacity_model = opacity_model_; + static constexpr bool enable_dust = false; }; template <> diff --git a/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp b/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp index 5010bb1ed..a27a9555c 100644 --- a/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp +++ b/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp @@ -42,6 +42,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp b/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp index 475c2a284..873de5f7d 100644 --- a/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp +++ b/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp @@ -44,6 +44,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadPulse/test_radiation_pulse.cpp b/src/problems/RadPulse/test_radiation_pulse.cpp index cde6a25cc..33a5b73ae 100644 --- a/src/problems/RadPulse/test_radiation_pulse.cpp +++ b/src/problems/RadPulse/test_radiation_pulse.cpp @@ -40,6 +40,7 @@ template <> struct RadSystem_Traits { 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 = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadShadow/test_radiation_shadow.cpp b/src/problems/RadShadow/test_radiation_shadow.cpp index f4e583ac2..91adf2cf0 100644 --- a/src/problems/RadShadow/test_radiation_shadow.cpp +++ b/src/problems/RadShadow/test_radiation_shadow.cpp @@ -41,6 +41,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadStreaming/test_radiation_streaming.cpp b/src/problems/RadStreaming/test_radiation_streaming.cpp index 1746e9e5b..f15211ca8 100644 --- a/src/problems/RadStreaming/test_radiation_streaming.cpp +++ b/src/problems/RadStreaming/test_radiation_streaming.cpp @@ -46,6 +46,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = 1.0; static constexpr double Erad_floor = initial_Erad; static constexpr int beta_order = 0; + static constexpr bool enable_dust = false; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp index 2c6894291..e401b034d 100644 --- a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp +++ b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp @@ -49,6 +49,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = 1.0; static constexpr double Erad_floor = initial_Erad; static constexpr int beta_order = 0; + static constexpr bool enable_dust = false; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadSuOlson/test_radiation_SuOlson.cpp b/src/problems/RadSuOlson/test_radiation_SuOlson.cpp index b6204e5db..e850bbb6e 100644 --- a/src/problems/RadSuOlson/test_radiation_SuOlson.cpp +++ b/src/problems/RadSuOlson/test_radiation_SuOlson.cpp @@ -50,6 +50,7 @@ template <> struct RadSystem_Traits { static constexpr double gamma = 5. / 3.; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 0; + static constexpr bool enable_dust = false; }; template <> struct quokka::EOS_Traits { diff --git a/src/problems/RadTophat/test_radiation_tophat.cpp b/src/problems/RadTophat/test_radiation_tophat.cpp index eff81e77d..2bb39f64b 100644 --- a/src/problems/RadTophat/test_radiation_tophat.cpp +++ b/src/problems/RadTophat/test_radiation_tophat.cpp @@ -46,6 +46,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 0; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadTube/test_radiation_tube.cpp b/src/problems/RadTube/test_radiation_tube.cpp index 80c9770a4..ee5a27977 100644 --- a/src/problems/RadTube/test_radiation_tube.cpp +++ b/src/problems/RadTube/test_radiation_tube.cpp @@ -70,6 +70,7 @@ template <> struct RadSystem_Traits { // 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; + static constexpr bool enable_dust = false; }; template <> diff --git a/src/problems/RadhydroBB/test_radhydro_bb.cpp b/src/problems/RadhydroBB/test_radhydro_bb.cpp index a969fe3ef..0d0463743 100644 --- a/src/problems/RadhydroBB/test_radhydro_bb.cpp +++ b/src/problems/RadhydroBB/test_radhydro_bb.cpp @@ -132,6 +132,7 @@ template <> struct RadSystem_Traits { static constexpr double energy_unit = nu_unit; static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; + static constexpr bool enable_dust = false; }; template <> diff --git a/src/problems/RadhydroPulse/test_radhydro_pulse.cpp b/src/problems/RadhydroPulse/test_radhydro_pulse.cpp index 181ed8821..9be49e793 100644 --- a/src/problems/RadhydroPulse/test_radhydro_pulse.cpp +++ b/src/problems/RadhydroPulse/test_radhydro_pulse.cpp @@ -49,6 +49,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; + static constexpr bool enable_dust = false; }; template <> struct RadSystem_Traits { static constexpr double c_light = c; @@ -56,6 +57,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp b/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp index 2415c3e84..a1ca6b38b 100644 --- a/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp +++ b/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp @@ -50,6 +50,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; + static constexpr bool enable_dust = false; }; template <> struct RadSystem_Traits { static constexpr double c_light = c; @@ -57,6 +58,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp b/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp index 8c6536e0f..5d93cfdd7 100644 --- a/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp +++ b/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp @@ -52,6 +52,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> struct RadSystem_Traits { static constexpr double c_light = c; @@ -59,6 +60,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 2; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp b/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp index b920b6fdf..d856a1f88 100644 --- a/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp +++ b/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp @@ -111,6 +111,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real { return kappa0; } @@ -181,6 +182,7 @@ template <> struct RadSystem_Traits { // static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum; // static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_full_spectrum; + static constexpr bool enable_dust = false; }; template <> diff --git a/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp b/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp index 7be6a2a14..e8fefd340 100644 --- a/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp +++ b/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp @@ -132,6 +132,7 @@ template <> struct RadSystem_Traits { static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; static constexpr int beta_order = 1; static constexpr OpacityModel opacity_model = opacity_model_; + static constexpr bool enable_dust = false; }; template <> struct RadSystem_Traits { static constexpr double c_light = c; @@ -141,6 +142,7 @@ template <> struct RadSystem_Traits { static constexpr bool compute_v_over_c_terms = true; static constexpr int beta_order = 1; static constexpr OpacityModel opacity_model = OpacityModel::single_group; + static constexpr bool enable_dust = false; }; AMREX_GPU_HOST_DEVICE diff --git a/src/problems/RadhydroShell/test_radhydro_shell.cpp b/src/problems/RadhydroShell/test_radhydro_shell.cpp index 13bfc1e86..649b79dac 100644 --- a/src/problems/RadhydroShell/test_radhydro_shell.cpp +++ b/src/problems/RadhydroShell/test_radhydro_shell.cpp @@ -56,6 +56,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> struct HydroSystem_Traits { diff --git a/src/problems/RadhydroShock/test_radhydro_shock.cpp b/src/problems/RadhydroShock/test_radhydro_shock.cpp index dd04e87a9..a5d4eb2d6 100644 --- a/src/problems/RadhydroShock/test_radhydro_shock.cpp +++ b/src/problems/RadhydroShock/test_radhydro_shock.cpp @@ -59,6 +59,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> struct quokka::EOS_Traits { diff --git a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp index 940054e7b..4c8148cba 100644 --- a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp +++ b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp @@ -60,6 +60,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; + static constexpr bool enable_dust = false; }; template <> struct quokka::EOS_Traits { diff --git a/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp b/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp index d45ac1246..d0654a0de 100644 --- a/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp +++ b/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp @@ -66,6 +66,7 @@ template <> struct RadSystem_Traits { // static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum; // static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_full_spectrum; + static constexpr bool enable_dust = false; }; template <> struct quokka::EOS_Traits { diff --git a/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp b/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp index c1085a384..02e5004e4 100644 --- a/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp +++ b/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp @@ -62,6 +62,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.0; static constexpr int beta_order = beta_order_; + static constexpr bool enable_dust = false; }; template <> struct Physics_Traits { diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 5ee692991..ed4e073dd 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -62,7 +62,6 @@ static constexpr double c_light_cgs_ = C::c_light; // cgs static constexpr double radiation_constant_cgs_ = C::a_rad; // cgs static constexpr double inf = std::numeric_limits::max(); -static constexpr bool use_dust = true; static const double Lambda_gd_0 = 1e6; // static const double Lambda_gd_0 = 2.5e-34; // erg cm^3 s^−1 K^−3/2 @@ -150,6 +149,8 @@ template class RadSystem : public HyperbolicSystem::beta_order; + static constexpr bool enable_dust_ = RadSystem_Traits::enable_dust; + static constexpr int nGroups_ = Physics_Traits::nGroups; static constexpr amrex::GpuArray radBoundaries_ = []() constexpr { if constexpr (nGroups_ > 1) { @@ -159,6 +160,7 @@ template class RadSystem : public HyperbolicSystem::Erad_floor / nGroups_; static constexpr OpacityModel opacity_model_ = []() constexpr { @@ -1448,7 +1450,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne AMREX_ASSERT(T_gas >= 0.); // dust temperature - if constexpr (!use_dust) { + if constexpr (!enable_dust_) { T_d = T_gas; } else { if (n == 0) { @@ -1635,7 +1637,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // M_0g dF0_dXg.fillin(cscale); // M_g0 - if constexpr (!use_dust) { + if constexpr (!enable_dust_) { dFg_dX0 = 1.0 / c_v * dEg_dT; } else { const double d_Td_d_T = 3. / 2. - T_d / (2. * T_gas); @@ -1679,7 +1681,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne } // END NEWTON-RAPHSON LOOP AMREX_ALWAYS_ASSERT_WITH_MESSAGE(n < maxIter, "Newton-Raphson iteration failed to converge!"); - std::cout << "Newton-Raphson converged after " << n << " it." << std::endl; + // std::cout << "Newton-Raphson converged after " << n << " it." << std::endl; AMREX_ALWAYS_ASSERT(Egas_guess > 0.0); AMREX_ALWAYS_ASSERT(min(EradVec_guess) >= 0.0); From b1a8109f1bdb768b0b94f449299d12d9359f2bcc Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 16 Aug 2024 13:27:44 +1000 Subject: [PATCH 19/25] github-advanced-security --- src/problems/RadDust/test_rad_dust.cpp | 2 +- src/problems/RadDustMG/test_rad_dust_MG.cpp | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/problems/RadDust/test_rad_dust.cpp b/src/problems/RadDust/test_rad_dust.cpp index 0b247ac95..b32861351 100644 --- a/src/problems/RadDust/test_rad_dust.cpp +++ b/src/problems/RadDust/test_rad_dust.cpp @@ -119,7 +119,7 @@ template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokk template <> void QuokkaSimulation::computeAfterTimestep() { - auto [position, values] = fextract(state_new_cc_[0], Geom(0), 0, 0.5); + auto [position, values] = fextract(state_new_cc_[0], Geom(0), 0, 0.5); // NOLINT if (amrex::ParallelDescriptor::IOProcessor()) { userData_.t_vec_.push_back(tNew_[0]); diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp index bbdf01ab6..3420311f2 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.cpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -1,5 +1,5 @@ /// \file test_rad_dust_MG.cpp -/// \brief Defines a multigroup test problem for gas-dust-radiation coupling in uniform medium. +/// \brief Defines a multigroup test problem for gas-dust-radiation coupling in uniform medium. /// #include "test_rad_dust_MG.hpp" @@ -72,7 +72,7 @@ template <> struct RadSystem_Traits { template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, const double rho, - const double /*Tgas*/) -> amrex::GpuArray, 2> + const double /*Tgas*/) -> amrex::GpuArray, 2> { amrex::GpuArray, 2> exponents_and_values{}; for (int i = 0; i < nGroups_ + 1; ++i) { @@ -93,8 +93,9 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiation(amrex } template <> -AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiationTempDerivative( - amrex::Real temperature, amrex::GpuArray const &boundaries) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto +RadSystem::ComputeThermalRadiationTempDerivative(amrex::Real temperature, + amrex::GpuArray const &boundaries) -> quokka::valarray { auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); const double d_power_dt = radiation_constant_; @@ -128,7 +129,7 @@ template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokk template <> void QuokkaSimulation::computeAfterTimestep() { - auto [position, values] = fextract(state_new_cc_[0], Geom(0), 0, 0.5); + auto [position, values] = fextract(state_new_cc_[0], Geom(0), 0, 0.5); // NOLINT if (amrex::ParallelDescriptor::IOProcessor()) { userData_.t_vec_.push_back(tNew_[0]); From d47f20ec375564e9e5f0c0ac7ab5c669e252bfed Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 16 Aug 2024 03:28:04 +0000 Subject: [PATCH 20/25] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/problems/RadDustMG/test_rad_dust_MG.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp index 3420311f2..c6c7946f1 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.cpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -1,5 +1,5 @@ /// \file test_rad_dust_MG.cpp -/// \brief Defines a multigroup test problem for gas-dust-radiation coupling in uniform medium. +/// \brief Defines a multigroup test problem for gas-dust-radiation coupling in uniform medium. /// #include "test_rad_dust_MG.hpp" @@ -72,7 +72,7 @@ template <> struct RadSystem_Traits { template <> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto RadSystem::DefineOpacityExponentsAndLowerValues(amrex::GpuArray /*rad_boundaries*/, const double rho, - const double /*Tgas*/) -> amrex::GpuArray, 2> + const double /*Tgas*/) -> amrex::GpuArray, 2> { amrex::GpuArray, 2> exponents_and_values{}; for (int i = 0; i < nGroups_ + 1; ++i) { @@ -93,9 +93,8 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiation(amrex } template <> -AMREX_GPU_HOST_DEVICE auto -RadSystem::ComputeThermalRadiationTempDerivative(amrex::Real temperature, - amrex::GpuArray const &boundaries) -> quokka::valarray +AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeThermalRadiationTempDerivative( + amrex::Real temperature, amrex::GpuArray const &boundaries) -> quokka::valarray { auto radEnergyFractions = ComputePlanckEnergyFractions(boundaries, temperature); const double d_power_dt = radiation_constant_; From 214aeeb9af11f7ebf496ea8542bb1d021bd0e86e Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 16 Aug 2024 14:05:25 +1000 Subject: [PATCH 21/25] add radiation.dust_gas_interaction_coeff to .in --- src/QuokkaSimulation.hpp | 4 +++- src/problems/RadDust/test_rad_dust.cpp | 2 -- src/problems/RadDustMG/test_rad_dust_MG.cpp | 2 -- src/radiation/radiation_system.hpp | 23 +++++++++------------ tests/RadDust.in | 6 ++++++ 5 files changed, 19 insertions(+), 18 deletions(-) diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index 9e890aad2..1da6b21de 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -124,6 +124,7 @@ template class QuokkaSimulation : public AMRSimulation void QuokkaSimulation::readParmParse() amrex::ParmParse rpp("radiation"); rpp.query("reconstruction_order", radiationReconstructionOrder_); rpp.query("cfl", radiationCflNumber_); + rpp.query("dust_gas_interaction_coeff", dustGasInteractionCoeff_); } } @@ -1803,7 +1805,7 @@ void QuokkaSimulation::operatorSplitSourceTerms(amrex::Array4::SetRadEnergySource(radEnergySource.array(), indexRange, dx, prob_lo, prob_hi, time + dt); // cell-centered source terms - RadSystem::AddSourceTerms(stateNew, radEnergySource.const_array(), indexRange, dt, stage); + RadSystem::AddSourceTerms(stateNew, radEnergySource.const_array(), indexRange, dt, stage, dustGasInteractionCoeff_); } template diff --git a/src/problems/RadDust/test_rad_dust.cpp b/src/problems/RadDust/test_rad_dust.cpp index b32861351..24c248bbc 100644 --- a/src/problems/RadDust/test_rad_dust.cpp +++ b/src/problems/RadDust/test_rad_dust.cpp @@ -142,7 +142,6 @@ auto problem_main() -> int // Problem parameters const int max_timesteps = 1e6; const double CFL_number_gas = 0.8; - const double CFL_number_rad = 8.0; // Boundary conditions constexpr int nvars = RadSystem::nvar_; @@ -159,7 +158,6 @@ auto problem_main() -> int sim.radiationReconstructionOrder_ = 3; // PPM sim.stopTime_ = max_time; - sim.radiationCflNumber_ = CFL_number_rad; sim.cflNumber_ = CFL_number_gas; sim.maxTimesteps_ = max_timesteps; sim.plotfileInterval_ = -1; diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp index 3420311f2..1291e9998 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.cpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -155,7 +155,6 @@ auto problem_main() -> int // Problem parameters const int max_timesteps = 1e6; const double CFL_number_gas = 0.8; - const double CFL_number_rad = 8.0; // Boundary conditions constexpr int nvars = RadSystem::nvar_; @@ -172,7 +171,6 @@ auto problem_main() -> int sim.radiationReconstructionOrder_ = 3; // PPM sim.stopTime_ = max_time; - sim.radiationCflNumber_ = CFL_number_rad; sim.cflNumber_ = CFL_number_gas; sim.maxTimesteps_ = max_timesteps; sim.plotfileInterval_ = -1; diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 2dd1e636a..68c9115dc 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -62,9 +62,6 @@ static constexpr double c_light_cgs_ = C::c_light; // cgs static constexpr double radiation_constant_cgs_ = C::a_rad; // cgs static constexpr double inf = std::numeric_limits::max(); -static const double Lambda_gd_0 = 1e6; -// static const double Lambda_gd_0 = 2.5e-34; // erg cm^3 s^−1 K^−3/2 - // enum for opacity_model enum class OpacityModel { single_group = 0, // user-defined opacity for each group, given as a function of density and temperature. @@ -207,7 +204,7 @@ template class RadSystem : public HyperbolicSystem const &prob_lo, amrex::GpuArray const &prob_hi, amrex::Real time); - static void AddSourceTerms(array_t &consVar, arrayconst_t &radEnergySource, amrex::Box const &indexRange, amrex::Real dt, int stage); + static void AddSourceTerms(array_t &consVar, arrayconst_t &radEnergySource, amrex::Box const &indexRange, amrex::Real dt, int stage, double dustGasCoeff); static void balanceMatterRadiation(arrayconst_t &consPrev, array_t &consNew, amrex::Box const &indexRange); @@ -269,7 +266,7 @@ template class RadSystem : public HyperbolicSystem const &boundaries) -> quokka::valarray; AMREX_GPU_HOST_DEVICE static auto ComputeDustTemperature(double T_gas, double T_d_init, double rho, quokka::valarray const &Erad, - amrex::GpuArray const &rad_boundaries, + double dustGasCoeff, amrex::GpuArray const &rad_boundaries, amrex::GpuArray const &rad_boundary_ratios) -> double; template @@ -1192,7 +1189,7 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxInDiffusionLimit(con template AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeDustTemperature(double const T_gas, double const T_d_init, double const rho, - quokka::valarray const &Erad, + quokka::valarray const &Erad, double dustGasCoeff, amrex::GpuArray const &rad_boundaries, amrex::GpuArray const &rad_boundary_ratios) -> double { @@ -1200,7 +1197,7 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeDustTemperature(double c quokka::valarray kappaEVec{}; const double num_density = rho / mean_molecular_mass_; - const double Lambda_compare = Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas) * T_gas; + const double Lambda_compare = dustGasCoeff * num_density * num_density * std::sqrt(T_gas) * T_gas; amrex::GpuArray alpha_quant_minus_one{}; for (int g = 0; g < nGroups_; ++g) { @@ -1239,14 +1236,14 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeDustTemperature(double c AMREX_ASSERT(!kappaEVec.hasnan()); const double LHS = c_light_ * rho * sum(kappaEVec * Erad - kappaPVec * fourPiBoverC) + - Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas) * (T_gas - T_d); + dustGasCoeff * num_density * num_density * std::sqrt(T_gas) * (T_gas - T_d); if (std::abs(LHS) < lambda_rel_tol * std::abs(Lambda_compare)) { break; } const auto d_fourpib_over_c_d_t = ComputeThermalRadiationTempDerivative(T_d, rad_boundaries); - const double dLHS_dTd = -c_light_ * rho * sum(kappaPVec * d_fourpib_over_c_d_t) - Lambda_gd_0 * num_density * num_density * std::sqrt(T_gas); + const double dLHS_dTd = -c_light_ * rho * sum(kappaPVec * d_fourpib_over_c_d_t) - dustGasCoeff * num_density * num_density * std::sqrt(T_gas); const double delta_T_d = LHS / dLHS_dTd; T_d -= delta_T_d; @@ -1263,7 +1260,7 @@ AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeDustTemperature(double c template void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEnergySource, amrex::Box const &indexRange, amrex::Real dt_radiation, - const int stage) + const int stage, double dustGasCoeff) { arrayconst_t &consPrev = consVar; // make read-only array_t &consNew = consVar; @@ -1455,11 +1452,11 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne T_d = T_gas; } else { if (n == 0) { - T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, radBoundaries_g_copy, + T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, dustGasCoeff, radBoundaries_g_copy, radBoundaryRatios_copy); } else { const auto Lambda_gd = sum(Rvec) / (dt * chat / c); - T_d = T_gas - Lambda_gd / (Lambda_gd_0 * num_den * num_den * std::sqrt(T_gas)); + T_d = T_gas - Lambda_gd / (dustGasCoeff * num_den * num_den * std::sqrt(T_gas)); } } AMREX_ASSERT(T_d >= 0.); @@ -1643,7 +1640,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne dFg_dX0 = 1.0 / c_v * dEg_dT; } else { const double d_Td_d_T = 3. / 2. - T_d / (2. * T_gas); - const double coeff_n = dt * Lambda_gd_0 * num_den * num_den / cscale; + const double coeff_n = dt * dustGasCoeff * num_den * num_den / cscale; 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; diff --git a/tests/RadDust.in b/tests/RadDust.in index 82d70ab3f..f46a110d2 100644 --- a/tests/RadDust.in +++ b/tests/RadDust.in @@ -5,6 +5,12 @@ geometry.prob_lo = 0.0 0.0 0.0 geometry.prob_hi = 1.0 1.0 1.0 geometry.is_periodic = 1 1 1 +# ***************************************************************** +# RADIATION +# ***************************************************************** +radiation.cfl = 8.0 +radiation.dust_gas_interaction_coeff = 1.0e6 + # ***************************************************************** # VERBOSITY # ***************************************************************** From 76a0f8d4eadc85808343e8c7c2ce873f0205c95a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 16 Aug 2024 04:05:52 +0000 Subject: [PATCH 22/25] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/QuokkaSimulation.hpp | 2 +- src/radiation/radiation_system.hpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/QuokkaSimulation.hpp b/src/QuokkaSimulation.hpp index 1da6b21de..a1278dc44 100644 --- a/src/QuokkaSimulation.hpp +++ b/src/QuokkaSimulation.hpp @@ -123,7 +123,7 @@ template class QuokkaSimulation : public AMRSimulation::nstartHyperbolic_; amrex::Real radiationCflNumber_ = 0.3; - int maxSubsteps_ = 10; // maximum number of radiation subcycles per hydro step + int maxSubsteps_ = 10; // maximum number of radiation subcycles per hydro step amrex::Real dustGasInteractionCoeff_ = 2.5e-34; // erg cm^3 s^−1 K^−3/2 bool computeReferenceSolution_ = false; diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 68c9115dc..4d4361d1e 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -204,7 +204,8 @@ template class RadSystem : public HyperbolicSystem const &prob_lo, amrex::GpuArray const &prob_hi, amrex::Real time); - static void AddSourceTerms(array_t &consVar, arrayconst_t &radEnergySource, amrex::Box const &indexRange, amrex::Real dt, int stage, double dustGasCoeff); + static void AddSourceTerms(array_t &consVar, arrayconst_t &radEnergySource, amrex::Box const &indexRange, amrex::Real dt, int stage, + double dustGasCoeff); static void balanceMatterRadiation(arrayconst_t &consPrev, array_t &consNew, amrex::Box const &indexRange); From 5fe36106b6886d4b50badf85604f2747a66970eb Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 16 Aug 2024 15:07:35 +1000 Subject: [PATCH 23/25] rename enable_dust --- src/problems/RadBeam/test_radiation_beam.cpp | 2 +- src/problems/RadDust/test_rad_dust.cpp | 2 +- src/problems/RadDustMG/test_rad_dust_MG.cpp | 2 +- src/problems/RadForce/test_radiation_force.cpp | 2 +- src/problems/RadMarshak/test_radiation_marshak.cpp | 2 +- .../test_radiation_marshak_asymptotic.cpp | 2 +- src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp | 2 +- .../RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp | 2 +- .../RadMatterCoupling/test_radiation_matter_coupling.cpp | 2 +- .../test_radiation_matter_coupling_rsla.cpp | 2 +- src/problems/RadPulse/test_radiation_pulse.cpp | 2 +- src/problems/RadShadow/test_radiation_shadow.cpp | 2 +- src/problems/RadStreaming/test_radiation_streaming.cpp | 2 +- src/problems/RadStreamingY/test_radiation_streaming_y.cpp | 2 +- src/problems/RadSuOlson/test_radiation_SuOlson.cpp | 2 +- src/problems/RadTophat/test_radiation_tophat.cpp | 2 +- src/problems/RadTube/test_radiation_tube.cpp | 2 +- src/problems/RadhydroBB/test_radhydro_bb.cpp | 2 +- src/problems/RadhydroPulse/test_radhydro_pulse.cpp | 4 ++-- src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp | 4 ++-- src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp | 4 ++-- .../test_radhydro_pulse_MG_const_kappa.cpp | 4 ++-- .../RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp | 4 ++-- src/problems/RadhydroShell/test_radhydro_shell.cpp | 2 +- src/problems/RadhydroShock/test_radhydro_shock.cpp | 2 +- src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp | 2 +- .../test_radhydro_shock_multigroup.cpp | 2 +- .../test_radhydro_uniform_advecting.cpp | 2 +- src/radiation/radiation_system.hpp | 6 +++--- 29 files changed, 36 insertions(+), 36 deletions(-) diff --git a/src/problems/RadBeam/test_radiation_beam.cpp b/src/problems/RadBeam/test_radiation_beam.cpp index 53947dc20..e760f96b9 100644 --- a/src/problems/RadBeam/test_radiation_beam.cpp +++ b/src/problems/RadBeam/test_radiation_beam.cpp @@ -38,7 +38,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadDust/test_rad_dust.cpp b/src/problems/RadDust/test_rad_dust.cpp index 24c248bbc..a89996a96 100644 --- a/src/problems/RadDust/test_rad_dust.cpp +++ b/src/problems/RadDust/test_rad_dust.cpp @@ -51,7 +51,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; - static constexpr bool enable_dust = true; + static constexpr bool enable_dust_gas_thermal_coupling_model = true; }; template <> struct Physics_Traits { diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp index 0f6f5aa2b..41d30ae7d 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.cpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -66,7 +66,7 @@ template <> struct RadSystem_Traits { static constexpr amrex::GpuArray::nGroups + 1> radBoundaries{1.0e-3, 0.1, 1.0, 10.0, 1.0e3}; // static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum; - static constexpr bool enable_dust = true; + static constexpr bool enable_dust_gas_thermal_coupling_model = true; }; template <> diff --git a/src/problems/RadForce/test_radiation_force.cpp b/src/problems/RadForce/test_radiation_force.cpp index fada3df80..180b093ff 100644 --- a/src/problems/RadForce/test_radiation_force.cpp +++ b/src/problems/RadForce/test_radiation_force.cpp @@ -71,7 +71,7 @@ template <> struct RadSystem_Traits { static constexpr double Erad_floor = 0.; static constexpr double energy_unit = C::ev2erg; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real { return 0.; } diff --git a/src/problems/RadMarshak/test_radiation_marshak.cpp b/src/problems/RadMarshak/test_radiation_marshak.cpp index 910f7e38a..e7902309a 100644 --- a/src/problems/RadMarshak/test_radiation_marshak.cpp +++ b/src/problems/RadMarshak/test_radiation_marshak.cpp @@ -42,7 +42,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 0; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp b/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp index 80cb7c386..44ca3b0f0 100644 --- a/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp +++ b/src/problems/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp @@ -39,7 +39,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = Erad_floor_; static constexpr int beta_order = 0; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp b/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp index fa9e08206..fa4a1f997 100644 --- a/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp +++ b/src/problems/RadMarshakCGS/test_radiation_marshak_cgs.cpp @@ -44,7 +44,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp b/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp index a4c3737ac..6c9b33f11 100644 --- a/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp +++ b/src/problems/RadMarshakVaytet/test_radiation_marshak_Vaytet.cpp @@ -124,7 +124,7 @@ template <> struct RadSystem_Traits { static constexpr double energy_unit = C::hplanck; // set boundary unit to Hz static constexpr amrex::GpuArray radBoundaries = group_edges_; static constexpr OpacityModel opacity_model = opacity_model_; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> diff --git a/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp b/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp index a27a9555c..43b171c77 100644 --- a/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp +++ b/src/problems/RadMatterCoupling/test_radiation_matter_coupling.cpp @@ -42,7 +42,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp b/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp index 873de5f7d..465fb7edf 100644 --- a/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp +++ b/src/problems/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp @@ -44,7 +44,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadPulse/test_radiation_pulse.cpp b/src/problems/RadPulse/test_radiation_pulse.cpp index 33a5b73ae..0df12a9e6 100644 --- a/src/problems/RadPulse/test_radiation_pulse.cpp +++ b/src/problems/RadPulse/test_radiation_pulse.cpp @@ -40,7 +40,7 @@ template <> struct RadSystem_Traits { 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 = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadShadow/test_radiation_shadow.cpp b/src/problems/RadShadow/test_radiation_shadow.cpp index 91adf2cf0..bd7a5f8dd 100644 --- a/src/problems/RadShadow/test_radiation_shadow.cpp +++ b/src/problems/RadShadow/test_radiation_shadow.cpp @@ -41,7 +41,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadStreaming/test_radiation_streaming.cpp b/src/problems/RadStreaming/test_radiation_streaming.cpp index f15211ca8..4bf372adb 100644 --- a/src/problems/RadStreaming/test_radiation_streaming.cpp +++ b/src/problems/RadStreaming/test_radiation_streaming.cpp @@ -46,7 +46,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = 1.0; static constexpr double Erad_floor = initial_Erad; static constexpr int beta_order = 0; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp index e401b034d..87eadc783 100644 --- a/src/problems/RadStreamingY/test_radiation_streaming_y.cpp +++ b/src/problems/RadStreamingY/test_radiation_streaming_y.cpp @@ -49,7 +49,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = 1.0; static constexpr double Erad_floor = initial_Erad; static constexpr int beta_order = 0; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real diff --git a/src/problems/RadSuOlson/test_radiation_SuOlson.cpp b/src/problems/RadSuOlson/test_radiation_SuOlson.cpp index e850bbb6e..4e6570cf9 100644 --- a/src/problems/RadSuOlson/test_radiation_SuOlson.cpp +++ b/src/problems/RadSuOlson/test_radiation_SuOlson.cpp @@ -50,7 +50,7 @@ template <> struct RadSystem_Traits { static constexpr double gamma = 5. / 3.; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 0; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct quokka::EOS_Traits { diff --git a/src/problems/RadTophat/test_radiation_tophat.cpp b/src/problems/RadTophat/test_radiation_tophat.cpp index 2bb39f64b..03c6a4595 100644 --- a/src/problems/RadTophat/test_radiation_tophat.cpp +++ b/src/problems/RadTophat/test_radiation_tophat.cpp @@ -46,7 +46,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = radiation_constant_cgs_; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 0; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadTube/test_radiation_tube.cpp b/src/problems/RadTube/test_radiation_tube.cpp index ee5a27977..a17fb0a01 100644 --- a/src/problems/RadTube/test_radiation_tube.cpp +++ b/src/problems/RadTube/test_radiation_tube.cpp @@ -70,7 +70,7 @@ template <> struct RadSystem_Traits { // 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; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> diff --git a/src/problems/RadhydroBB/test_radhydro_bb.cpp b/src/problems/RadhydroBB/test_radhydro_bb.cpp index 0d0463743..a0a631d5e 100644 --- a/src/problems/RadhydroBB/test_radhydro_bb.cpp +++ b/src/problems/RadhydroBB/test_radhydro_bb.cpp @@ -132,7 +132,7 @@ template <> struct RadSystem_Traits { static constexpr double energy_unit = nu_unit; static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> diff --git a/src/problems/RadhydroPulse/test_radhydro_pulse.cpp b/src/problems/RadhydroPulse/test_radhydro_pulse.cpp index 9be49e793..c101927ac 100644 --- a/src/problems/RadhydroPulse/test_radhydro_pulse.cpp +++ b/src/problems/RadhydroPulse/test_radhydro_pulse.cpp @@ -49,7 +49,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct RadSystem_Traits { static constexpr double c_light = c; @@ -57,7 +57,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp b/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp index a1ca6b38b..5f6bd12d7 100644 --- a/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp +++ b/src/problems/RadhydroPulseDyn/test_radhydro_pulse_dyn.cpp @@ -50,7 +50,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct RadSystem_Traits { static constexpr double c_light = c; @@ -58,7 +58,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = beta_order_; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp b/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp index 5d93cfdd7..83e4fdf71 100644 --- a/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp +++ b/src/problems/RadhydroPulseGrey/test_radhydro_pulse_grey.cpp @@ -52,7 +52,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct RadSystem_Traits { static constexpr double c_light = c; @@ -60,7 +60,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 2; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp b/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp index d856a1f88..d7f4fd79e 100644 --- a/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp +++ b/src/problems/RadhydroPulseMGconst/test_radhydro_pulse_MG_const_kappa.cpp @@ -111,7 +111,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = erad_floor; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputePlanckOpacity(const double /*rho*/, const double /*Tgas*/) -> amrex::Real { return kappa0; } @@ -182,7 +182,7 @@ template <> struct RadSystem_Traits { // static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum; // static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_full_spectrum; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> diff --git a/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp b/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp index e8fefd340..c0f028dc8 100644 --- a/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp +++ b/src/problems/RadhydroPulseMGint/test_radhydro_pulse_MG_int.cpp @@ -132,7 +132,7 @@ template <> struct RadSystem_Traits { static constexpr amrex::GpuArray radBoundaries = rad_boundaries_; static constexpr int beta_order = 1; static constexpr OpacityModel opacity_model = opacity_model_; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct RadSystem_Traits { static constexpr double c_light = c; @@ -142,7 +142,7 @@ template <> struct RadSystem_Traits { static constexpr bool compute_v_over_c_terms = true; static constexpr int beta_order = 1; static constexpr OpacityModel opacity_model = OpacityModel::single_group; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; AMREX_GPU_HOST_DEVICE diff --git a/src/problems/RadhydroShell/test_radhydro_shell.cpp b/src/problems/RadhydroShell/test_radhydro_shell.cpp index 649b79dac..9acbd5352 100644 --- a/src/problems/RadhydroShell/test_radhydro_shell.cpp +++ b/src/problems/RadhydroShell/test_radhydro_shell.cpp @@ -56,7 +56,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct HydroSystem_Traits { diff --git a/src/problems/RadhydroShock/test_radhydro_shock.cpp b/src/problems/RadhydroShock/test_radhydro_shock.cpp index a5d4eb2d6..75cf62af4 100644 --- a/src/problems/RadhydroShock/test_radhydro_shock.cpp +++ b/src/problems/RadhydroShock/test_radhydro_shock.cpp @@ -59,7 +59,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct quokka::EOS_Traits { diff --git a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp index 4c8148cba..4e1ae9999 100644 --- a/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp +++ b/src/problems/RadhydroShockCGS/test_radhydro_shock_cgs.cpp @@ -60,7 +60,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.; static constexpr int beta_order = 1; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct quokka::EOS_Traits { diff --git a/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp b/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp index d0654a0de..cb2c03655 100644 --- a/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp +++ b/src/problems/RadhydroShockMultigroup/test_radhydro_shock_multigroup.cpp @@ -66,7 +66,7 @@ template <> struct RadSystem_Traits { // static constexpr OpacityModel opacity_model = OpacityModel::piecewise_constant_opacity; static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_fixed_slope_spectrum; // static constexpr OpacityModel opacity_model = OpacityModel::PPL_opacity_full_spectrum; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct quokka::EOS_Traits { diff --git a/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp b/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp index 02e5004e4..79f55b4db 100644 --- a/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp +++ b/src/problems/RadhydroUniformAdvecting/test_radhydro_uniform_advecting.cpp @@ -62,7 +62,7 @@ template <> struct RadSystem_Traits { static constexpr double radiation_constant = a_rad; static constexpr double Erad_floor = 0.0; static constexpr int beta_order = beta_order_; - static constexpr bool enable_dust = false; + static constexpr bool enable_dust_gas_thermal_coupling_model = false; }; template <> struct Physics_Traits { diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 68c9115dc..4119bc44c 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -146,7 +146,7 @@ template class RadSystem : public HyperbolicSystem::beta_order; - static constexpr bool enable_dust_ = RadSystem_Traits::enable_dust; + static constexpr bool enable_dust_gas_thermal_coupling_model_ = RadSystem_Traits::enable_dust_gas_thermal_coupling_model; static constexpr int nGroups_ = Physics_Traits::nGroups; static constexpr amrex::GpuArray radBoundaries_ = []() constexpr { @@ -1448,7 +1448,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne AMREX_ASSERT(T_gas >= 0.); // dust temperature - if constexpr (!enable_dust_) { + if constexpr (!enable_dust_gas_thermal_coupling_model_) { T_d = T_gas; } else { if (n == 0) { @@ -1636,7 +1636,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne // M_0g dF0_dXg.fillin(cscale); // M_g0 - if constexpr (!enable_dust_) { + if constexpr (!enable_dust_gas_thermal_coupling_model_) { dFg_dX0 = 1.0 / c_v * dEg_dT; } else { const double d_Td_d_T = 3. / 2. - T_d / (2. * T_gas); From 050761cd0c735b3726751ca2e44b87d0cd5bc12e Mon Sep 17 00:00:00 2001 From: ChongChong He Date: Fri, 16 Aug 2024 15:15:33 +1000 Subject: [PATCH 24/25] fix first-capture error --- src/radiation/radiation_system.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index 0b0a19b2d..bbafb5f1c 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -1266,6 +1266,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne arrayconst_t &consPrev = consVar; // make read-only array_t &consNew = consVar; auto dt = dt_radiation; + const double dustGasCoeff_local = dustGasCoeff; if (stage == 2) { dt = (1.0 - IMEX_a32) * dt_radiation; } @@ -1453,11 +1454,11 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne T_d = T_gas; } else { if (n == 0) { - T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, dustGasCoeff, radBoundaries_g_copy, + T_d = ComputeDustTemperature(T_gas, T_gas, rho, EradVec_guess, dustGasCoeff_local, radBoundaries_g_copy, radBoundaryRatios_copy); } else { const auto Lambda_gd = sum(Rvec) / (dt * chat / c); - T_d = T_gas - Lambda_gd / (dustGasCoeff * num_den * num_den * std::sqrt(T_gas)); + T_d = T_gas - Lambda_gd / (dustGasCoeff_local * num_den * num_den * std::sqrt(T_gas)); } } AMREX_ASSERT(T_d >= 0.); @@ -1641,7 +1642,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne dFg_dX0 = 1.0 / c_v * dEg_dT; } else { const double d_Td_d_T = 3. / 2. - T_d / (2. * T_gas); - const double coeff_n = dt * dustGasCoeff * num_den * num_den / cscale; + const double coeff_n = dt * dustGasCoeff_local * num_den * num_den / cscale; 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; From 452c141d74f62285b7c9d85c5b41fe4c9b12d2e3 Mon Sep 17 00:00:00 2001 From: chongchonghe Date: Fri, 16 Aug 2024 15:46:48 +1000 Subject: [PATCH 25/25] fix first-capture error --- src/problems/RadDust/test_rad_dust.cpp | 1 - src/problems/RadDustMG/test_rad_dust_MG.cpp | 1 - src/radiation/radiation_system.hpp | 18 +++++++----------- 3 files changed, 7 insertions(+), 13 deletions(-) diff --git a/src/problems/RadDust/test_rad_dust.cpp b/src/problems/RadDust/test_rad_dust.cpp index a89996a96..b6e23e976 100644 --- a/src/problems/RadDust/test_rad_dust.cpp +++ b/src/problems/RadDust/test_rad_dust.cpp @@ -21,7 +21,6 @@ constexpr double v0 = 0.0; constexpr double chi0 = 10000.0; constexpr double T0 = 1.0; -constexpr double T_equi = 0.7680325; constexpr double rho0 = 1.0; constexpr double a_rad = 1.0; constexpr double mu = 1.0; diff --git a/src/problems/RadDustMG/test_rad_dust_MG.cpp b/src/problems/RadDustMG/test_rad_dust_MG.cpp index 41d30ae7d..5fa6045f2 100644 --- a/src/problems/RadDustMG/test_rad_dust_MG.cpp +++ b/src/problems/RadDustMG/test_rad_dust_MG.cpp @@ -21,7 +21,6 @@ constexpr double v0 = 0.0; constexpr double chi0 = 10000.0; constexpr double T0 = 1.0; -constexpr double T_equi = 0.7680325; constexpr double rho0 = 1.0; constexpr double a_rad = 1.0; constexpr double mu = 1.0; diff --git a/src/radiation/radiation_system.hpp b/src/radiation/radiation_system.hpp index bbafb5f1c..cabcf063f 100644 --- a/src/radiation/radiation_system.hpp +++ b/src/radiation/radiation_system.hpp @@ -1266,20 +1266,11 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne arrayconst_t &consPrev = consVar; // make read-only array_t &consNew = consVar; auto dt = dt_radiation; - const double dustGasCoeff_local = dustGasCoeff; if (stage == 2) { dt = (1.0 - IMEX_a32) * dt_radiation; } amrex::GpuArray radBoundaries_g = radBoundaries_; - amrex::GpuArray radBoundaryRatios{}; - if constexpr (nGroups_ > 1) { - if constexpr (static_cast(opacity_model_) > 0) { - for (int g = 0; g < nGroups_; ++g) { - radBoundaryRatios[g] = radBoundaries_g[g + 1] / radBoundaries_g[g]; - } - } - } // Add source terms @@ -1290,6 +1281,7 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { const double c = c_light_; const double chat = c_hat_; + const double dustGasCoeff_local = dustGasCoeff; // load fluid properties const double rho = consPrev(i, j, k, gasDensity_index); @@ -1355,8 +1347,12 @@ void RadSystem::AddSourceTerms(array_t &consVar, arrayconst_t &radEne amrex::GpuArray radBoundaryRatios_copy{}; for (int g = 0; g < nGroups_ + 1; ++g) { radBoundaries_g_copy[g] = radBoundaries_g[g]; - if (g < nGroups_) { - radBoundaryRatios_copy[g] = radBoundaryRatios[g]; + } + if constexpr (nGroups_ > 1) { + if constexpr (static_cast(opacity_model_) > 0) { + for (int g = 0; g < nGroups_; ++g) { + radBoundaryRatios_copy[g] = radBoundaries_g_copy[g + 1] / radBoundaries_g_copy[g]; + } } }