From d1ad36edbeef9488a34da00c0888d1c0a109bb49 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Mon, 13 Jul 2026 07:36:49 -0600 Subject: [PATCH 01/14] Update full coupling --- src/radiation/moments/matter_coupling.hpp | 2098 ++++++++++++++++++--- 1 file changed, 1803 insertions(+), 295 deletions(-) diff --git a/src/radiation/moments/matter_coupling.hpp b/src/radiation/moments/matter_coupling.hpp index c471f3a4..bece8327 100644 --- a/src/radiation/moments/matter_coupling.hpp +++ b/src/radiation/moments/matter_coupling.hpp @@ -1,15 +1,19 @@ //======================================================================================== -// (C) (or copyright) 2023-2025. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2023-2026. Triad National Security, LLC. All rights +// reserved. // -// This program was produced under U.S. Government contract 89233218CNA000001 for Los -// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC -// for the U.S. Department of Energy/National Nuclear Security Administration. All rights -// in the program are reserved by Triad National Security, LLC, and the U.S. Department -// of Energy/National Nuclear Security Administration. The Government is granted for -// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide -// license in this material to reproduce, prepare derivative works, distribute copies to -// the public, perform publicly and display publicly, and to permit others to do so. +// This program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National Nuclear +// Security Administration. All rights in the program are reserved by Triad +// National Security, LLC, and the U.S. Department of Energy/National Nuclear +// Security Administration. The Government is granted for itself and others +// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license +// in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do +// so. //======================================================================================== +// This file was modified in part with the assistance of generative AI. #ifndef RADIATION_MOMENTS_MATTER_COUPLING_HPP_ #define RADIATION_MOMENTS_MATTER_COUPLING_HPP_ @@ -28,177 +32,1060 @@ using ArtemisUtils::VI; namespace Moments { -//---------------------------------------------------------------------------------------- -//! \fn TaskStatus Moments::MatterCouplingSimpleImpl -//! \brief Implementation for simple radiation-matter coupling source -template -TaskStatus MatterCouplingSimpleImpl(MeshData *u0, const Real dt) { - PARTHENON_INSTRUMENT - using parthenon::MakePackDescriptor; - using parthenon::variable_names::any; - auto pm = u0->GetParentPointer(); - auto &resolved_pkgs = pm->resolved_packages; +KOKKOS_INLINE_FUNCTION Real FiniteOrZero(const Real x) { + return std::isfinite(x) ? x : 0.0; +} - // Extract gas package and params - auto &gas_pkg = pm->packages.Get("gas"); - auto eos_d = gas_pkg->template Param("eos_d"); - auto opac_d = gas_pkg->template Param("opacity_d"); - auto scat_d = gas_pkg->template Param("scattering_d"); - auto dflr = gas_pkg->template Param("dfloor"); - auto de_switch = gas_pkg->template Param("de_switch"); +KOKKOS_INLINE_FUNCTION Real CouplingFluxMagnitude(const std::array &F) { + return std::sqrt(SQR(F[0]) + SQR(F[1]) + SQR(F[2])); +} - // Extract radiation and moments package and params - auto &moments_pkg = pm->packages.Get("moments"); - const auto chat = moments_pkg->template Param("chat"); - const auto c = moments_pkg->template Param("c"); - const auto arad = moments_pkg->template Param("arad"); - const auto tfloor = moments_pkg->template Param("tfloor"); - const auto Bfloor = arad * SQR(SQR(tfloor)); - const auto efloor = Bfloor; - const auto outer_max = moments_pkg->template Param("outer_iteration_max"); - const auto inner_max = moments_pkg->template Param("inner_iteration_max"); - const auto outer_tol = moments_pkg->template Param("outer_iteration_tol"); - const auto inner_tol = moments_pkg->template Param("inner_iteration_tol"); - const auto &cpars = - pm->packages.Get("artemis")->template Param("coord_params"); +KOKKOS_INLINE_FUNCTION bool ComputeCouplingEnergyCoefficients( + const Real sigp, const Real sigs, const Real g, const Real g2, + const Real beta2, const Real bdbdp, const Real bdf, Real &ca, Real &cb, + Real &cd) { + // Algebraically equivalent to + // ca = g * (sigp + sigs - g2 * sigs * (1 + bdbdp)); + // cb = g * sigp; + // cd = -g * bdf * (sigp + sigs - 2 * g2 * sigs); + // but avoids subtracting nearly equal O(sigs) terms when beta << 1. + ca = g * (sigp - g2 * sigs * (beta2 + bdbdp)); + cb = g * sigp; + cd = g * bdf * (g2 * sigs * (1.0 + beta2) - sigp); + return std::isfinite(ca) && std::isfinite(cb) && std::isfinite(cd); +} - const auto fatal_if_unconverged = - moments_pkg->template Param("fatal_if_unconverged"); +KOKKOS_INLINE_FUNCTION bool ComputeCouplingEquilibriumEnergy( + const Real E0, const Real B, const Real ca, const Real cb, + const Real cd, Real &Eeq) { + // The implicit radiation-energy equation is + // (E - E0) + ca E - cb B + cd = 0. + // Evaluate its solution using coefficient ratios. This remains accurate in + // the optically thick limit, where ca and cb may be O(1e10) and directly + // evaluating ca*E - cb*B loses all useful digits near LTE. + const Real denom = 1.0 + ca; + if (!std::isfinite(denom) || std::abs(denom) <= Fuzz()) + return false; + const Real inv_denom = 1.0 / denom; + Eeq = E0 * inv_denom + (cb * inv_denom) * B - cd * inv_denom; + return std::isfinite(Eeq); +} - // Extract rotating frame quantities - Real om0 = 0.0; - Real qshear = 0.0; - Real gm_bg = 0.0; - if (pm->packages.Get("artemis")->template Param("do_orbital_advection") || - pm->packages.Get("artemis")->template Param("do_rotating_frame")) { - auto &rframe_pkg = pm->packages.Get("rotating_frame"); - qshear = rframe_pkg->template Param("qshear"); - om0 = rframe_pkg->template Param("omega"); - gm_bg = rframe_pkg->template Param("gm"); +KOKKOS_INLINE_FUNCTION std::array +ProjectCouplingFluxToEnergy(const std::array &Fin, const Real E) { + const Real fmag = CouplingFluxMagnitude(Fin); + if (!std::isfinite(fmag) || !std::isfinite(E) || E <= 0.0) + return {0.0, 0.0, 0.0}; + if (fmag <= E || fmag <= Fuzz()) + return Fin; + const Real scale = E / fmag; + return {scale * Fin[0], scale * Fin[1], scale * Fin[2]}; +} + +KOKKOS_INLINE_FUNCTION bool SolveCouplingDense3x3( + Real A[3][3], Real rhs[3], std::array &x) { + // Small partial-pivoting Gaussian elimination used by the local momentum + // Newton solve. The system is only 3x3, so forming and factorizing it is + // cheaper and more robust than attempting an analytic inverse. + for (int col = 0; col < 3; ++col) { + int pivot = col; + Real pivot_abs = std::abs(A[col][col]); + for (int row = col + 1; row < 3; ++row) { + const Real candidate = std::abs(A[row][col]); + if (candidate > pivot_abs) { + pivot = row; + pivot_abs = candidate; + } + } + if (!std::isfinite(pivot_abs) || pivot_abs <= Fuzz()) + return false; + if (pivot != col) { + for (int j = col; j < 3; ++j) { + const Real tmp = A[col][j]; + A[col][j] = A[pivot][j]; + A[pivot][j] = tmp; + } + const Real tmp = rhs[col]; + rhs[col] = rhs[pivot]; + rhs[pivot] = tmp; + } + + const Real inv_pivot = 1.0 / A[col][col]; + for (int row = col + 1; row < 3; ++row) { + const Real factor = A[row][col] * inv_pivot; + A[row][col] = 0.0; + for (int j = col + 1; j < 3; ++j) + A[row][j] -= factor * A[col][j]; + rhs[row] -= factor * rhs[col]; + } } - const bool do_raytrace = - pm->packages.Get("artemis")->template Param("do_raytrace"); - // Packing and indexing - static auto desc = parthenon::MakePackDescriptor< - rad::cons::energy, rad::cons::flux, gas::cons::density, gas::cons::momentum, - gas::cons::internal_energy, gas::cons::total_energy, gas::src::energy>( - resolved_pkgs.get()); - const auto v0 = desc.GetPack(u0); - static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); - auto vg = desc_g.GetPack(u0); - const auto ib = u0->GetBoundsI(IndexDomain::interior); - const auto jb = u0->GetBoundsJ(IndexDomain::interior); - const auto kb = u0->GetBoundsK(IndexDomain::interior); + for (int row = 2; row >= 0; --row) { + Real value = rhs[row]; + for (int j = row + 1; j < 3; ++j) + value -= A[row][j] * x[j]; + if (!std::isfinite(A[row][row]) || + std::abs(A[row][row]) <= Fuzz()) + return false; + x[row] = value / A[row][row]; + if (!std::isfinite(x[row])) + return false; + } + return true; +} - parthenon::par_for( - DEFAULT_LOOP_PATTERN, "MatterCoupling", DevExecSpace(), 0, u0->NumBlocks() - 1, - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, - KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - geometry::Coords coords(cpars, v0.GetCoordinates(b), k, j, i); - const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); - // y = U^(0) + dt S(y) +template +KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualBeta( + const std::array &beta, + const std::array &beta0, + const std::array &Fr0, const Real E0, const Real efloor, + const Real B, const Real dEg, const Real Q, const Real dens, + const Real eref, const Real c, const Real chat, const Real sigp, + const Real sigs, const Real energy_floor_scale, + std::array &F, Real &E, Real &dEk, + std::array &residual, Real &residual_norm) { + // Enforce gas+radiation pseudo-momentum conservation algebraically: + // rho c (beta-beta0) + (eref/chat) (F-Fr0) = 0, + // where F is normalized by c*eref. This relation is the source of the + // instability in a direct Picard update when eref/(rho*c*chat) is large. + const Real mu = eref / (dens * c * chat); + if (!std::isfinite(mu) || mu <= 0.0) + return false; + + const Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); + if (!std::isfinite(beta2) || beta2 >= 1.0) + return false; + + std::array delta_F{0.0, 0.0, 0.0}; + Real delta_F2 = 0.0; + Real beta0_dot_delta_F = 0.0; + for (int d = 0; d < 3; ++d) { + F[d] = Fr0[d] + (beta0[d] - beta[d]) / mu; + delta_F[d] = F[d] - Fr0[d]; + delta_F2 += SQR(delta_F[d]); + beta0_dot_delta_F += beta0[d] * delta_F[d]; + } - // U^(0) values - const Real dens = v0(b, gas::cons::density(), k, j, i); - Real Q = 0.0; - if (do_raytrace) Q = dt * v0(b, gas::src::energy(), k, j, i); - Real e0 = v0(b, gas::cons::internal_energy(), k, j, i); - const auto vb = RotatingFrame::BackgroundVelocity( - qshear, om0, gm_bg, coords.GetCellCenter(vg, b, k, j, i)); - std::array v{ - vb[0] + v0(b, gas::cons::momentum(0), k, j, i) / (hx[0] * dens), - vb[1] + v0(b, gas::cons::momentum(1), k, j, i) / (hx[1] * dens), - vb[2] + v0(b, gas::cons::momentum(2), k, j, i) / (hx[2] * dens)}; - const Real Er0 = v0(b, rad::cons::energy(), k, j, i); - std::array Fr0{v0(b, rad::cons::flux(0), k, j, i) / hx[0], - v0(b, rad::cons::flux(1), k, j, i) / hx[1], - v0(b, rad::cons::flux(2), k, j, i) / hx[2]}; - - // Note(AMD): There is some floating point difference between the internal energy - // used to compute the temperature and the internal energy obtained from that - // temperature: T = eos_d.TemperatureFromDensityInternalEnergy(dens, eg/dens); eg - // /= dens * eos_d.InternalEnergyFromDensityTemperature(dens,T) - // - // Because of this, zero opacity problems will not result in zero change as - // expected. Thus, we recalculate the internal and total energies from the - // temperature. This does not affect energy conservation because at the end of the - // step we update the energy with an increment. + // From beta-beta0 = -mu*(F-Fr0), evaluate the kinetic-energy + // change without subtracting two nearly equal beta^2 values: + // dEk = (c/chat)[-beta0.dF + 0.5*mu*|dF|^2]. + dEk = c / chat * + (-beta0_dot_delta_F + 0.5 * mu * delta_F2); + E = E0 + chat / c * (Q - dEg - dEk); + const Real floor_slop = + 128.0 * std::numeric_limits::epsilon() * + std::max(1.0, std::max(std::abs(E0), std::abs(E))); + if (!std::isfinite(E) || E < efloor - floor_slop || E <= 0.0) + return false; + E = std::max(E, efloor); + + const Real fmag = CouplingFluxMagnitude(F); + const Real realizability_slop = + 256.0 * std::numeric_limits::epsilon() * + std::max(1.0, E); + if (!std::isfinite(fmag) || fmag > E + realizability_slop) + return false; + + const Real g2 = 1.0 / (1.0 - beta2); + const Real g = std::sqrt(g2); + const auto fedd = EddingtonTensor( + {F[0] / E, F[1] / E, F[2] / E}); + const std::array bdp{ + beta[0] * fedd[TensIdx::X11] + beta[1] * fedd[TensIdx::X12] + + beta[2] * fedd[TensIdx::X13], + beta[0] * fedd[TensIdx::X12] + beta[1] * fedd[TensIdx::X22] + + beta[2] * fedd[TensIdx::X23], + beta[0] * fedd[TensIdx::X13] + beta[1] * fedd[TensIdx::X23] + + beta[2] * fedd[TensIdx::X33]}; + const Real bdbdp = + beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; + + const Real sigf = sigp + sigs; + const Real a = g * sigf; + const Real bcoef = 2.0 * g2 * g * sigs; + const Real d1 = + g * (sigp * B + g2 * sigs * (1.0 + bdbdp) * E); + const Real d2 = g * sigf * E; + const std::array rhs{ + Fr0[0] + d1 * beta[0] + d2 * bdp[0], + Fr0[1] + d1 * beta[1] + d2 * bdp[1], + Fr0[2] + d1 * beta[2] + d2 * bdp[2]}; + const auto Fsolve = SolveRadFlux(1.0 + a, bcoef, beta, rhs); + if (!std::isfinite(Fsolve[0]) || !std::isfinite(Fsolve[1]) || + !std::isfinite(Fsolve[2])) + return false; + const auto Ftarget = ProjectCouplingFluxToEnergy(Fsolve, E); + + const Real scale = std::max( + energy_floor_scale, + std::max(E, std::max(CouplingFluxMagnitude(Fr0), + std::max(fmag, CouplingFluxMagnitude(Ftarget))))); + residual_norm = 0.0; + for (int d = 0; d < 3; ++d) { + residual[d] = F[d] - Ftarget[d]; + residual_norm = + std::max(residual_norm, std::abs(residual[d]) / scale); + } + return std::isfinite(residual_norm); +} - Real T = eos_d.TemperatureFromDensityInternalEnergy(dens, e0 / dens); - e0 = dens * eos_d.InternalEnergyFromDensityTemperature(dens, T); - Real e = e0; - Real B = arad * SQR(SQR(T)); +template +KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( + const std::array &beta0, + const std::array &Fr0, const Real E0, const Real efloor, + const Real B, const Real dEg, const Real Q, const Real dens, + const Real eref, const Real c, const Real chat, const Real sigp, + const Real sigs, const Real energy_floor_scale, const Real tolerance, + std::array &beta, std::array &F, Real &E, + Real &dEk, Real &momentum_error, int &iterations) { + const Real mu = eref / (dens * c * chat); + if (!std::isfinite(mu) || mu <= 0.0 || !std::isfinite(sigp) || + !std::isfinite(sigs) || sigp < 0.0 || sigs < 0.0) + return false; + + const Real solve_tol = + std::max(tolerance, 64.0 * std::numeric_limits::epsilon()); + constexpr int max_iterations = 32; + constexpr int max_line_search = 20; + const Real fd_factor = + std::pow(std::numeric_limits::epsilon(), 1.0 / 3.0); + + std::array residual{0.0, 0.0, 0.0}; + std::array best_beta = beta; + std::array best_F = F; + Real best_E = E; + Real best_dEk = dEk; + Real best_error = std::numeric_limits::max(); + + for (iterations = 1; iterations <= max_iterations; ++iterations) { + Real current_error = std::numeric_limits::max(); + if (!EvaluateCouplingMomentumEnergyResidualBeta( + beta, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, + chat, sigp, sigs, energy_floor_scale, F, E, dEk, residual, + current_error)) + break; + + if (current_error < best_error) { + best_error = current_error; + best_beta = beta; + best_F = F; + best_E = E; + best_dEk = dEk; + } + if (current_error <= solve_tol) { + momentum_error = current_error; + return true; + } + + Real J[3][3]{{0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0}}; + bool jacobian_valid = true; + for (int col = 0; col < 3; ++col) { + const Real h = std::max( + 32.0 * std::numeric_limits::epsilon(), + fd_factor * std::max(1.0e-3, std::abs(beta[col]))); + auto beta_plus = beta; + auto beta_minus = beta; + beta_plus[col] += h; + beta_minus[col] -= h; + + std::array Fplus, Fminus; + std::array Rplus, Rminus; + Real Eplus = E; + Real Eminus = E; + Real dEkplus = dEk; + Real dEkminus = dEk; + Real err_plus = 0.0; + Real err_minus = 0.0; + const bool plus_valid = + EvaluateCouplingMomentumEnergyResidualBeta( + beta_plus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, + eref, c, chat, sigp, sigs, energy_floor_scale, Fplus, + Eplus, dEkplus, Rplus, err_plus); + const bool minus_valid = + EvaluateCouplingMomentumEnergyResidualBeta( + beta_minus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, + eref, c, chat, sigp, sigs, energy_floor_scale, Fminus, + Eminus, dEkminus, Rminus, err_minus); + + if (plus_valid && minus_valid) { + for (int row = 0; row < 3; ++row) + J[row][col] = (Rplus[row] - Rminus[row]) / (2.0 * h); + } else if (plus_valid) { + for (int row = 0; row < 3; ++row) + J[row][col] = (Rplus[row] - residual[row]) / h; + } else if (minus_valid) { + for (int row = 0; row < 3; ++row) + J[row][col] = (residual[row] - Rminus[row]) / h; + } else { + jacobian_valid = false; + break; + } + } + + std::array step{0.0, 0.0, 0.0}; + if (jacobian_valid) { + Real rhs[3]{-residual[0], -residual[1], -residual[2]}; + jacobian_valid = SolveCouplingDense3x3(J, rhs, step); + } + + // If the numerical Jacobian is singular, use the analytically motivated + // damped Picard direction. In the diffusion limit the unstable Picard + // gain is approximately -(4/3) mu E; the denominator below removes that + // stiffness while retaining the correct fixed point. + if (!jacobian_valid) { + const Real relax = 1.0 / (1.0 + (4.0 / 3.0) * mu * + std::max(E, energy_floor_scale)); + for (int d = 0; d < 3; ++d) + step[d] = relax * mu * residual[d]; + } + + const Real step_mag = + std::sqrt(SQR(step[0]) + SQR(step[1]) + SQR(step[2])); + if (!std::isfinite(step_mag) || step_mag <= Fuzz()) + break; + if (step_mag > 0.25) { + const Real scale = 0.25 / step_mag; + for (int d = 0; d < 3; ++d) + step[d] *= scale; + } + + bool accepted = false; + Real alpha = 1.0; + for (int ls = 0; ls < max_line_search; ++ls) { + std::array beta_trial{ + beta[0] + alpha * step[0], beta[1] + alpha * step[1], + beta[2] + alpha * step[2]}; + std::array Ftrial, Rtrial; + Real Etrial = E; + Real dEktrial = dEk; + Real trial_error = std::numeric_limits::max(); + const bool trial_valid = + EvaluateCouplingMomentumEnergyResidualBeta( + beta_trial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, + eref, c, chat, sigp, sigs, energy_floor_scale, Ftrial, + Etrial, dEktrial, Rtrial, trial_error); + if (trial_valid && + trial_error < current_error * (1.0 - 1.0e-4 * alpha)) { + beta = beta_trial; + F = Ftrial; + E = Etrial; + dEk = dEktrial; + accepted = true; + break; + } + alpha *= 0.5; + } + + if (!accepted) { + // Retry with the stiffness-aware Picard direction even when the Newton + // Jacobian existed but produced a poor globalization step. + const Real relax = 1.0 / (1.0 + (4.0 / 3.0) * mu * + std::max(E, energy_floor_scale)); + for (int d = 0; d < 3; ++d) + step[d] = relax * mu * residual[d]; + alpha = 1.0; + for (int ls = 0; ls < max_line_search; ++ls) { + std::array beta_trial{ + beta[0] + alpha * step[0], beta[1] + alpha * step[1], + beta[2] + alpha * step[2]}; + std::array Ftrial, Rtrial; + Real Etrial = E; + Real dEktrial = dEk; + Real trial_error = std::numeric_limits::max(); + const bool trial_valid = + EvaluateCouplingMomentumEnergyResidualBeta( + beta_trial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, + eref, c, chat, sigp, sigs, energy_floor_scale, Ftrial, + Etrial, dEktrial, Rtrial, trial_error); + if (trial_valid && trial_error < current_error) { + beta = beta_trial; + F = Ftrial; + E = Etrial; + dEk = dEktrial; + accepted = true; + break; + } + alpha *= 0.5; + } + } + if (!accepted) + break; + } - Real E = Er0; - Real etot = e + c / chat * E; + beta = best_beta; + F = best_F; + E = best_E; + dEk = best_dEk; + momentum_error = best_error; + return best_error <= solve_tol; +} - // S(y) = sigma*chat*(E-B) - // - int inner_iter = 0; - Real inner_err = 0.; - for (inner_iter = 0; inner_iter < inner_max; inner_iter++) { - T = std::pow(B / arad, 0.25); - e = eos_d.InternalEnergyFromDensityTemperature(dens, T) * dens; - const Real Cv = dens * eos_d.SpecificHeatFromDensityTemperature(dens, T); - const Real a = chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); - const Real fleck = FleckFactor(arad, T, Cv); - - const Real Ri = a * (E - B); - const Real Fi = (e - e0) - c / chat * Ri - Q; - const Real Fr = (E - Er0) + Ri; - const Real idet = 1. / (1. + a + c / chat * fleck * a); - Real dE = ((1. + c / chat * fleck * a) * (-Fr) + a * (-fleck * Fi)) * idet; - Real dB = ((c / chat * fleck * a) * (-Fr) + (1. + a) * (-fleck * Fi)) * idet; - - Real Enew = E + dE; - E = (Enew < efloor) ? efloor : Enew; - Real Bnew = B + dB; - B = (Bnew < Bfloor) ? Bfloor : Bnew; - - inner_err = std::max((std::abs(Fi) / etot), (c / chat * std::abs(Fr) / etot)); - if (inner_err <= inner_tol) { - break; - } + +// In the high-matter-inertia limit mu = eref/(rho*c*chat) << 1, the +// physically resolvable momentum unknown is the radiation-flux change, not +// beta. A finite flux correction corresponds to a beta correction of order +// mu*dF, which can be smaller than one floating-point ulp of beta. Solving in +// beta then cannot represent the root even though F remains well resolved. +template +KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualFlux( + const std::array &F, + const std::array &beta0, + const std::array &Fr0, const Real E0, const Real efloor, + const Real B, const Real dEg, const Real Q, const Real dens, + const Real eref, const Real c, const Real chat, const Real sigp, + const Real sigs, const Real energy_floor_scale, + std::array &beta, Real &E, Real &dEk, + std::array &residual, Real &residual_norm) { + const Real mu = eref / (dens * c * chat); + if (!std::isfinite(mu) || mu <= 0.0) + return false; + + std::array delta_F{0.0, 0.0, 0.0}; + Real delta_F2 = 0.0; + Real beta0_dot_delta_F = 0.0; + for (int d = 0; d < 3; ++d) { + delta_F[d] = F[d] - Fr0[d]; + beta[d] = beta0[d] - mu * delta_F[d]; + delta_F2 += SQR(delta_F[d]); + beta0_dot_delta_F += beta0[d] * delta_F[d]; + } + + const Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); + if (!std::isfinite(beta2) || beta2 >= 1.0) + return false; + + // This form remains accurate even when beta-beta0 is below the resolution + // of a stored velocity component. + dEk = c / chat * + (-beta0_dot_delta_F + 0.5 * mu * delta_F2); + E = E0 + chat / c * (Q - dEg - dEk); + const Real floor_slop = + 128.0 * std::numeric_limits::epsilon() * + std::max(1.0, std::max(std::abs(E0), std::abs(E))); + if (!std::isfinite(E) || E < efloor - floor_slop || E <= 0.0) + return false; + E = std::max(E, efloor); + + const Real fmag = CouplingFluxMagnitude(F); + const Real realizability_slop = + 256.0 * std::numeric_limits::epsilon() * + std::max(1.0, E); + if (!std::isfinite(fmag) || fmag > E + realizability_slop) + return false; + + const Real g2 = 1.0 / (1.0 - beta2); + const Real g = std::sqrt(g2); + const auto fedd = EddingtonTensor( + {F[0] / E, F[1] / E, F[2] / E}); + const std::array bdp{ + beta[0] * fedd[TensIdx::X11] + beta[1] * fedd[TensIdx::X12] + + beta[2] * fedd[TensIdx::X13], + beta[0] * fedd[TensIdx::X12] + beta[1] * fedd[TensIdx::X22] + + beta[2] * fedd[TensIdx::X23], + beta[0] * fedd[TensIdx::X13] + beta[1] * fedd[TensIdx::X23] + + beta[2] * fedd[TensIdx::X33]}; + const Real bdbdp = + beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; + + const Real sigf = sigp + sigs; + const Real a = g * sigf; + const Real bcoef = 2.0 * g2 * g * sigs; + const Real d1 = + g * (sigp * B + g2 * sigs * (1.0 + bdbdp) * E); + const Real d2 = g * sigf * E; + const std::array rhs{ + Fr0[0] + d1 * beta[0] + d2 * bdp[0], + Fr0[1] + d1 * beta[1] + d2 * bdp[1], + Fr0[2] + d1 * beta[2] + d2 * bdp[2]}; + const auto Fsolve = SolveRadFlux(1.0 + a, bcoef, beta, rhs); + if (!std::isfinite(Fsolve[0]) || !std::isfinite(Fsolve[1]) || + !std::isfinite(Fsolve[2])) + return false; + const auto Ftarget = ProjectCouplingFluxToEnergy(Fsolve, E); + + const Real scale = std::max( + energy_floor_scale, + std::max(E, std::max(CouplingFluxMagnitude(Fr0), + std::max(fmag, CouplingFluxMagnitude(Ftarget))))); + residual_norm = 0.0; + for (int d = 0; d < 3; ++d) { + residual[d] = F[d] - Ftarget[d]; + residual_norm = + std::max(residual_norm, std::abs(residual[d]) / scale); + } + return std::isfinite(residual_norm); +} + +template +KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( + const std::array &beta0, + const std::array &Fr0, const Real E0, const Real efloor, + const Real B, const Real dEg, const Real Q, const Real dens, + const Real eref, const Real c, const Real chat, const Real sigp, + const Real sigs, const Real energy_floor_scale, const Real tolerance, + std::array &beta, std::array &F, Real &E, + Real &dEk, Real &momentum_error, int &iterations) { + const Real mu = eref / (dens * c * chat); + if (!std::isfinite(mu) || mu <= 0.0 || !std::isfinite(sigp) || + !std::isfinite(sigs) || sigp < 0.0 || sigs < 0.0) + return false; + + const Real solve_tol = + std::max(tolerance, 64.0 * std::numeric_limits::epsilon()); + constexpr int max_iterations = 32; + constexpr int max_line_search = 20; + const Real fd_factor = + std::pow(std::numeric_limits::epsilon(), 1.0 / 3.0); + + std::array residual{0.0, 0.0, 0.0}; + std::array best_beta = beta; + std::array best_F = F; + Real best_E = E; + Real best_dEk = dEk; + Real best_error = std::numeric_limits::max(); + + for (iterations = 1; iterations <= max_iterations; ++iterations) { + Real current_error = std::numeric_limits::max(); + if (!EvaluateCouplingMomentumEnergyResidualFlux( + F, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, + chat, sigp, sigs, energy_floor_scale, beta, E, dEk, + residual, current_error)) + break; + + if (current_error < best_error) { + best_error = current_error; + best_beta = beta; + best_F = F; + best_E = E; + best_dEk = dEk; + } + if (current_error <= solve_tol) { + momentum_error = current_error; + return true; + } + + Real J[3][3]{{0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0}}; + bool jacobian_valid = true; + for (int col = 0; col < 3; ++col) { + const Real variable_scale = std::max( + energy_floor_scale, + std::max(std::abs(F[col]), + std::max(std::abs(Fr0[col]), 1.0e-6 * E))); + const Real h = std::max( + 64.0 * std::numeric_limits::epsilon() * + std::max(1.0, std::abs(F[col])), + fd_factor * variable_scale); + auto Fplus = F; + auto Fminus = F; + Fplus[col] += h; + Fminus[col] -= h; + + std::array beta_plus, beta_minus; + std::array Rplus, Rminus; + Real Eplus = E; + Real Eminus = E; + Real dEkplus = dEk; + Real dEkminus = dEk; + Real err_plus = 0.0; + Real err_minus = 0.0; + const bool plus_valid = + EvaluateCouplingMomentumEnergyResidualFlux( + Fplus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, + eref, c, chat, sigp, sigs, energy_floor_scale, + beta_plus, Eplus, dEkplus, Rplus, err_plus); + const bool minus_valid = + EvaluateCouplingMomentumEnergyResidualFlux( + Fminus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, + eref, c, chat, sigp, sigs, energy_floor_scale, + beta_minus, Eminus, dEkminus, Rminus, err_minus); + + if (plus_valid && minus_valid) { + for (int row = 0; row < 3; ++row) + J[row][col] = (Rplus[row] - Rminus[row]) / (2.0 * h); + } else if (plus_valid) { + for (int row = 0; row < 3; ++row) + J[row][col] = (Rplus[row] - residual[row]) / h; + } else if (minus_valid) { + for (int row = 0; row < 3; ++row) + J[row][col] = (residual[row] - Rminus[row]) / h; + } else { + jacobian_valid = false; + break; + } + } + + std::array step{0.0, 0.0, 0.0}; + if (jacobian_valid) { + Real rhs[3]{-residual[0], -residual[1], -residual[2]}; + jacobian_valid = SolveCouplingDense3x3(J, rhs, step); + } + if (!jacobian_valid) { + for (int d = 0; d < 3; ++d) + step[d] = -residual[d]; + } + + const Real step_mag = + std::sqrt(SQR(step[0]) + SQR(step[1]) + SQR(step[2])); + const Real state_scale = std::max( + energy_floor_scale, + std::max(E, std::max(CouplingFluxMagnitude(F), + CouplingFluxMagnitude(Fr0)))); + if (!std::isfinite(step_mag) || + step_mag <= 64.0 * std::numeric_limits::epsilon() * + std::max(1.0, state_scale)) + break; + if (step_mag > 0.25 * state_scale) { + const Real scale = 0.25 * state_scale / step_mag; + for (int d = 0; d < 3; ++d) + step[d] *= scale; + } + + bool accepted = false; + Real alpha = 1.0; + for (int ls = 0; ls < max_line_search; ++ls) { + std::array Ftrial{ + F[0] + alpha * step[0], F[1] + alpha * step[1], + F[2] + alpha * step[2]}; + std::array beta_trial, Rtrial; + Real Etrial = E; + Real dEktrial = dEk; + Real trial_error = std::numeric_limits::max(); + const bool trial_valid = + EvaluateCouplingMomentumEnergyResidualFlux( + Ftrial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, + eref, c, chat, sigp, sigs, energy_floor_scale, + beta_trial, Etrial, dEktrial, Rtrial, trial_error); + if (trial_valid && + trial_error < current_error * (1.0 - 1.0e-4 * alpha)) { + F = Ftrial; + beta = beta_trial; + E = Etrial; + dEk = dEktrial; + accepted = true; + break; + } + alpha *= 0.5; + } + + if (!accepted) { + // A fixed-point step in F is well conditioned in this branch. + for (int d = 0; d < 3; ++d) + step[d] = -residual[d]; + alpha = 1.0; + for (int ls = 0; ls < max_line_search; ++ls) { + std::array Ftrial{ + F[0] + alpha * step[0], F[1] + alpha * step[1], + F[2] + alpha * step[2]}; + std::array beta_trial, Rtrial; + Real Etrial = E; + Real dEktrial = dEk; + Real trial_error = std::numeric_limits::max(); + const bool trial_valid = + EvaluateCouplingMomentumEnergyResidualFlux( + Ftrial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, + eref, c, chat, sigp, sigs, energy_floor_scale, + beta_trial, Etrial, dEktrial, Rtrial, trial_error); + if (trial_valid && trial_error < current_error) { + F = Ftrial; + beta = beta_trial; + E = Etrial; + dEk = dEktrial; + accepted = true; + break; } - if ((inner_iter == inner_max) && (fatal_if_unconverged)) { - printf("(%d,%d,%d,%d) %lg > %lg after %d iterations\n", b, k, j, i, inner_err, - inner_tol, inner_max); - PARTHENON_FAIL("Radiation matter coupling did not converge!"); + alpha *= 0.5; + } + } + if (!accepted) + break; + } + + beta = best_beta; + F = best_F; + E = best_E; + dEk = best_dEk; + momentum_error = best_error; + return best_error <= solve_tol; +} + +template +KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergy( + const std::array &beta0, + const std::array &Fr0, const Real E0, const Real efloor, + const Real B, const Real dEg, const Real Q, const Real dens, + const Real eref, const Real c, const Real chat, const Real sigp, + const Real sigs, const Real energy_floor_scale, const Real tolerance, + std::array &beta, std::array &F, Real &E, + Real &dEk, Real &momentum_error, int &iterations) { + const Real mu = eref / (dens * c * chat); + // Choose the nonlinear variable from the relative pseudo-inertia. For + // mu <= 1 the gas velocity change can be much less resolvable than the + // radiation-flux change, so solve in F. For mu > 1 solve in beta to avoid + // magnifying flux perturbations into large velocity perturbations. + if (mu <= 1.0) { + return SolveCouplingMomentumEnergyFlux( + beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, + sigp, sigs, energy_floor_scale, tolerance, beta, F, E, dEk, + momentum_error, iterations); + } + return SolveCouplingMomentumEnergyBeta( + beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, + sigp, sigs, energy_floor_scale, tolerance, beta, F, E, dEk, + momentum_error, iterations); +} + +template +KOKKOS_INLINE_FUNCTION bool EvaluateCouplingInnerScalarResidual( + const Real eg, const Real eg0, const Real dEk, const Real Q, + const Real E0, const Real efloor, const Real Bfloor, const Real dens, + const Real eref, const Real arad, const Real tfloor, const Real dt, + const Real chat, const Real c, const Real g, const Real g2, + const Real beta2, const Real bdbdp, const Real bdf, const EOSType &eos, + const OpacityType &opacity, const ScatteringType &scattering, Real &E, + Real &B, Real &residual) { + // Exact reduced-speed-of-light total-energy conservation at fixed velocity: + // dEg + dEk + (c/chat) dEr = Q. + E = E0 + chat / c * (Q - dEk - (eg - eg0)); + const Real floor_slop = + 64.0 * std::numeric_limits::epsilon() * + std::max(1.0, std::max(std::abs(E0), std::abs(E))); + if (!std::isfinite(E) || E < efloor - floor_slop) + return false; + E = std::max(E, efloor); + + const Real eint = eg * eref / dens; + const Real T = std::max( + tfloor, eos.TemperatureFromDensityInternalEnergy(dens, eint)); + B = std::max(Bfloor, arad * SQR(SQR(T)) / eref); + const Real sigp = + chat * dt * opacity.PlanckMeanAbsorptionCoefficient(dens, T); + const Real sigs = + chat * dt * scattering.RosselandMeanTotalScatteringCoefficient(dens, T); + if (!std::isfinite(T) || !std::isfinite(B) || !std::isfinite(sigp) || + !std::isfinite(sigs) || sigp < 0.0 || sigs < 0.0) + return false; + + Real ca = 0.0; + Real cb = 0.0; + Real cd = 0.0; + Real Eeq = E; + if (!ComputeCouplingEnergyCoefficients(sigp, sigs, g, g2, beta2, + bdbdp, bdf, ca, cb, cd) || + !ComputeCouplingEquilibriumEnergy(E0, B, ca, cb, cd, Eeq)) + return false; + + // Compare two independently constructed radiation energies: + // E: exact total-energy conservation, + // Eeq: the implicit radiation source equation. + // Do not evaluate (1+ca)E - (E0+cbB-cd) here. In an optically thick + // near-LTE cell that expression subtracts O(sigp) numbers and its absolute + // roundoff can be orders of magnitude larger than the requested tolerance. + residual = E - Eeq; + return std::isfinite(residual); +} + +KOKKOS_INLINE_FUNCTION Real CouplingInnerScalarError( + const Real eg, const Real eg0, const Real dEk, const Real Q, + const Real E, const Real E0, const Real chat, const Real c, + const Real energy_floor_scale, const Real residual) { + const Real scale = std::max( + energy_floor_scale, + std::max( + std::max(std::abs(eg0), std::abs(eg)), + std::max(std::abs(Q), + std::max(c / chat * (std::abs(E0) + std::abs(E)), + std::abs(dEk))))); + return c / chat * std::abs(residual) / scale; +} + +KOKKOS_INLINE_FUNCTION bool CouplingResidualChangesSign(const Real a, + const Real b) { + return (a <= 0.0 && b >= 0.0) || (a >= 0.0 && b <= 0.0); +} + +template +KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( + const Real eg0, const Real dEk, const Real Q, const Real E0, + const Real efloor, const Real Bfloor, const Real dens, const Real eref, + const Real arad, const Real tfloor, const Real dt, const Real chat, + const Real c, const Real g, const Real g2, const Real beta2, + const Real bdbdp, const Real bdf, const Real energy_floor_scale, + const Real tolerance, + const EOSType &eos, const OpacityType &opacity, + const ScatteringType &scattering, Real &E, Real &B, Real &inner_err, + int &iterations) { + const Real eg_floor = + dens * eos.InternalEnergyFromDensityTemperature(dens, tfloor) / eref; + const Real eg_ceiling = + eg0 + Q - dEk + c / chat * (E0 - efloor); + if (!std::isfinite(eg_floor) || !std::isfinite(eg_ceiling) || + eg_ceiling < eg_floor) + return false; + + const Real Tguess = std::max( + tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); + Real eg_guess = + dens * eos.InternalEnergyFromDensityTemperature(dens, Tguess) / eref; + eg_guess = std::max(eg_floor, std::min(eg_ceiling, eg_guess)); + + Real Eguess = E; + Real Bguess = B; + Real Rguess = 0.0; + bool guess_valid = EvaluateCouplingInnerScalarResidual( + eg_guess, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, + tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, + scattering, + Eguess, Bguess, Rguess); + + Real best_eg = eg_guess; + Real best_E = Eguess; + Real best_B = Bguess; + Real best_error = + guess_valid + ? CouplingInnerScalarError(eg_guess, eg0, dEk, Q, Eguess, E0, + chat, c, energy_floor_scale, Rguess) + : std::numeric_limits::max(); + if (guess_valid && best_error <= tolerance) { + E = Eguess; + B = Bguess; + inner_err = best_error; + iterations = 1; + return true; + } + + Real Elo = E; + Real Blo = B; + Real Rlo = 0.0; + const bool lo_valid = EvaluateCouplingInnerScalarResidual( + eg_floor, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, + tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, + scattering, + Elo, Blo, Rlo); + if (lo_valid) { + const Real err = CouplingInnerScalarError(eg_floor, eg0, dEk, Q, Elo, E0, chat, c, + energy_floor_scale, Rlo); + if (err < best_error) { + best_eg = eg_floor; + best_E = Elo; + best_B = Blo; + best_error = err; + } + } + + Real Ehi = E; + Real Bhi = B; + Real Rhi = 0.0; + const bool hi_valid = EvaluateCouplingInnerScalarResidual( + eg_ceiling, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, + tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, + scattering, + Ehi, Bhi, Rhi); + if (hi_valid) { + const Real err = CouplingInnerScalarError(eg_ceiling, eg0, dEk, Q, Ehi, E0, chat, c, + energy_floor_scale, Rhi); + if (err < best_error) { + best_eg = eg_ceiling; + best_E = Ehi; + best_B = Bhi; + best_error = err; + } + } + + const bool global_bracketed = + lo_valid && hi_valid && CouplingResidualChangesSign(Rlo, Rhi); + bool bracketed = false; + Real bracket_lo = eg_floor; + Real bracket_hi = eg_ceiling; + Real bracket_Rlo = Rlo; + Real bracket_Rhi = Rhi; + + // Search outward from the current state before using the full admissible + // interval. Geometric spacing is important in low-density cells: a large + // change in B = a_r T^4 can correspond to an extremely small change in gas + // energy, while the upper energy bound may be many orders of magnitude away. + if (guess_valid) { + constexpr int scan_points = 64; + Real left_eg = eg_guess; + Real left_R = Rguess; + Real right_eg = eg_guess; + Real right_R = Rguess; + Real weight = 5.42101086242752217004e-20; // 2^-64 + for (int n = 1; n <= scan_points && !bracketed; ++n) { + weight = (n == scan_points) ? 1.0 : 2.0 * weight; + + const Real next_left = eg_guess - weight * (eg_guess - eg_floor); + Real Eleft = E; + Real Bleft = B; + Real Rleft = 0.0; + const bool left_valid = EvaluateCouplingInnerScalarResidual( + next_left, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, + tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, + scattering, Eleft, Bleft, Rleft); + if (left_valid) { + const Real err = CouplingInnerScalarError( + next_left, eg0, dEk, Q, Eleft, E0, chat, c, + energy_floor_scale, Rleft); + if (err < best_error) { + best_eg = next_left; + best_E = Eleft; + best_B = Bleft; + best_error = err; + } + if (CouplingResidualChangesSign(Rleft, left_R)) { + bracket_lo = next_left; + bracket_hi = left_eg; + bracket_Rlo = Rleft; + bracket_Rhi = left_R; + bracketed = true; + break; + } + left_eg = next_left; + left_R = Rleft; + } + + const Real next_right = eg_guess + weight * (eg_ceiling - eg_guess); + Real Eright = E; + Real Bright = B; + Real Rright = 0.0; + const bool right_valid = EvaluateCouplingInnerScalarResidual( + next_right, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, + tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, + scattering, Eright, Bright, Rright); + if (right_valid) { + const Real err = CouplingInnerScalarError( + next_right, eg0, dEk, Q, Eright, E0, chat, c, + energy_floor_scale, Rright); + if (err < best_error) { + best_eg = next_right; + best_E = Eright; + best_B = Bright; + best_error = err; + } + if (CouplingResidualChangesSign(right_R, Rright)) { + bracket_lo = right_eg; + bracket_hi = next_right; + bracket_Rlo = right_R; + bracket_Rhi = Rright; + bracketed = true; + break; } - T = std::pow(B / arad, 0.25); - e = eos_d.InternalEnergyFromDensityTemperature(dens, T) * dens; - const Real dEg = e - e0; - Real a = chat * dt * - (opac_d.RosselandMeanAbsorptionCoefficient(dens, T) + - scat_d.RosselandMeanTotalScatteringCoefficient(dens, T)); - std::array dF{-a / (1. + a) * Fr0[0], -a / (1. + a) * Fr0[1], - -a / (1. + a) * Fr0[2]}; - const Real icc = -1. / (c * chat * dens); - std::array dv{-icc * dF[0], -icc * dF[1], -icc * dF[2]}; - std::array vn{v[0] + dv[0], v[1] + dv[1], v[2] + dv[2]}; - - const Real dEk = 0.5 * dens * - ((SQR(vn[0]) - SQR(v[0])) + (SQR(vn[1]) - SQR(v[1])) + - (SQR(vn[2]) - SQR(v[2]))); + right_eg = next_right; + right_R = Rright; + } + } + } - // Update state vector (both gas and radiation) - v0(b, rad::cons::energy(), k, j, i) += chat / c * (dEk - dEg); - v0(b, gas::cons::total_energy(), k, j, i) += dEg + dEk; - v0(b, gas::cons::internal_energy(), k, j, i) += dEg; - v0(b, rad::cons::flux(0), k, j, i) += dF[0] * hx[0]; - v0(b, gas::cons::momentum(0), k, j, i) += dv[0] * dens * hx[0]; - v0(b, rad::cons::flux(1), k, j, i) += dF[1] * hx[1]; - v0(b, gas::cons::momentum(1), k, j, i) += dv[1] * dens * hx[1]; - v0(b, rad::cons::flux(2), k, j, i) += dF[2] * hx[2]; - v0(b, gas::cons::momentum(2), k, j, i) += dv[2] * dens * hx[2]; - }); - return TaskStatus::complete; + if (!bracketed && global_bracketed) { + bracket_lo = eg_floor; + bracket_hi = eg_ceiling; + bracket_Rlo = Rlo; + bracket_Rhi = Rhi; + bracketed = true; + } + + iterations = 0; + if (bracketed) { + for (iterations = 1; iterations <= 96; ++iterations) { + // Regula falsi when it lies safely inside the bracket; otherwise use the + // bisection midpoint. This retains bracketing while accelerating smooth + // cases. + Real eg_trial = 0.5 * (bracket_lo + bracket_hi); + const Real denom = bracket_Rhi - bracket_Rlo; + if (std::isfinite(denom) && std::abs(denom) > Fuzz()) { + const Real eg_secant = + (bracket_lo * bracket_Rhi - bracket_hi * bracket_Rlo) / denom; + const Real width = bracket_hi - bracket_lo; + if (eg_secant > bracket_lo + 0.1 * width && + eg_secant < bracket_hi - 0.1 * width) + eg_trial = eg_secant; + } + + Real Etrial = E; + Real Btrial = B; + Real Rtrial = 0.0; + const bool trial_valid = EvaluateCouplingInnerScalarResidual( + eg_trial, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, + tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, + scattering, Etrial, Btrial, Rtrial); + if (!trial_valid) + break; + + const Real err = CouplingInnerScalarError( + eg_trial, eg0, dEk, Q, Etrial, E0, chat, c, energy_floor_scale, + Rtrial); + if (err < best_error) { + best_eg = eg_trial; + best_E = Etrial; + best_B = Btrial; + best_error = err; + } + if (err <= tolerance) + break; + + if (CouplingResidualChangesSign(bracket_Rlo, Rtrial)) { + bracket_hi = eg_trial; + bracket_Rhi = Rtrial; + } else { + bracket_lo = eg_trial; + bracket_Rlo = Rtrial; + } + + const Real width_scale = std::max( + energy_floor_scale, + std::abs(bracket_lo) + std::abs(bracket_hi)); + if ((bracket_hi - bracket_lo) / width_scale <= + 8.0 * std::numeric_limits::epsilon()) + break; + } + } + + E = best_E; + B = best_B; + inner_err = best_error; + const Real attainable_tolerance = + std::max(tolerance, 32.0 * std::numeric_limits::epsilon()); + return std::isfinite(best_eg) && std::isfinite(E) && std::isfinite(B) && + std::isfinite(inner_err) && inner_err <= attainable_tolerance; +} + +template +KOKKOS_INLINE_FUNCTION bool EvaluateCouplingEnergyOnlyResidual( + const Real eg, const Real eg0, const Real Q, const Real E0, + const Real Emin, const Real dens, const Real eref, const Real arad, + const Real tfloor, const Real dt, const Real chat, const Real c, + const EOSType &eos, const OpacityType &opacity, Real &E, Real &B, + Real &residual) { + E = E0 + chat / c * (Q - (eg - eg0)); + if (!std::isfinite(E) || E < Emin) + return false; + + const Real eint = eg * eref / dens; + const Real T = std::max( + tfloor, eos.TemperatureFromDensityInternalEnergy(dens, eint)); + B = arad * SQR(SQR(T)) / eref; + const Real sigp = + chat * dt * opacity.PlanckMeanAbsorptionCoefficient(dens, T); + if (!std::isfinite(T) || !std::isfinite(B) || !std::isfinite(sigp) || + sigp < 0.0) + return false; + + // This fallback intentionally retains only thermal absorption/emission. + // Scattering energy exchange is work associated with the momentum source, + // which is frozen in this path. Solve the implicit radiation equation in + // ratio form; sigp*(E-B) is not a resolvable residual when sigp >> 1 and + // E and B agree to machine precision. + const Real inv_denom = 1.0 / (1.0 + sigp); + const Real Eeq = E0 * inv_denom + (sigp * inv_denom) * B; + residual = E - Eeq; + return std::isfinite(Eeq) && std::isfinite(residual); } //---------------------------------------------------------------------------------------- -//! \fn TaskStatus Moments::MatterCouplingSimpleImpl +//! \fn TaskStatus Moments::MatterCouplingFullSingleImpl //! \brief Implementation for "full" radiation-matter coupling source template TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { @@ -214,18 +1101,25 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { auto opac_d = gas_pkg->template Param("opacity_d"); auto scat_d = gas_pkg->template Param("scattering_d"); auto dflr = gas_pkg->template Param("dfloor"); - auto de_switch = gas_pkg->template Param("de_switch"); // Extract radiation package and params auto &moments_pkg = pm->packages.Get("moments"); const auto chat = moments_pkg->template Param("chat"); const auto c = moments_pkg->template Param("c"); const auto arad = moments_pkg->template Param("arad"); + const auto rad_efloor = moments_pkg->template Param("efloor"); const auto tfloor = moments_pkg->template Param("tfloor"); - const auto outer_max = moments_pkg->template Param("outer_iteration_max"); - const auto inner_max = moments_pkg->template Param("inner_iteration_max"); - const auto outer_tol = moments_pkg->template Param("outer_iteration_tol"); - const auto inner_tol = moments_pkg->template Param("inner_iteration_tol"); + const auto Bfloor_phys = arad * SQR(SQR(tfloor)); + const auto outer_max = + moments_pkg->template Param("outer_iteration_max"); + const auto inner_max = + moments_pkg->template Param("inner_iteration_max"); + const auto outer_tol = + moments_pkg->template Param("outer_iteration_tol"); + const auto inner_tol = + moments_pkg->template Param("inner_iteration_tol"); + const Real nonlinear_roundoff_tol = + 64.0 * std::numeric_limits::epsilon(); const auto fatal_if_unconverged = moments_pkg->template Param("fatal_if_unconverged"); @@ -233,7 +1127,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { Real om0 = 0.0; Real qshear = 0.0; Real gm_bg = 0.0; - if (pm->packages.Get("artemis")->template Param("do_orbital_advection") || + if (pm->packages.Get("artemis")->template Param( + "do_orbital_advection") || pm->packages.Get("artemis")->template Param("do_rotating_frame")) { auto &rframe_pkg = pm->packages.Get("rotating_frame"); qshear = rframe_pkg->template Param("qshear"); @@ -241,51 +1136,71 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { gm_bg = rframe_pkg->template Param("gm"); } const auto &cpars = - pm->packages.Get("artemis")->template Param("coord_params"); + pm->packages.Get("artemis")->template Param( + "coord_params"); const bool do_raytrace = pm->packages.Get("artemis")->template Param("do_raytrace"); // Packing and indexing - static auto desc = parthenon::MakePackDescriptor< - rad::cons::energy, rad::cons::flux, gas::cons::density, gas::cons::momentum, - gas::cons::internal_energy, gas::cons::total_energy, gas::src::energy>( - resolved_pkgs.get()); + static auto desc = + parthenon::MakePackDescriptor( + resolved_pkgs.get()); const auto v0 = desc.GetPack(u0); - static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); + static auto desc_g = + MakePackDescriptor(resolved_pkgs.get()); auto vg = desc_g.GetPack(u0); const auto ib = u0->GetBoundsI(IndexDomain::interior); const auto jb = u0->GetBoundsJ(IndexDomain::interior); const auto kb = u0->GetBoundsK(IndexDomain::interior); parthenon::par_for( - DEFAULT_LOOP_PATTERN, "MatterCoupling", DevExecSpace(), 0, u0->NumBlocks() - 1, - kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + DEFAULT_LOOP_PATTERN, "MatterCoupling", DevExecSpace(), 0, + u0->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { geometry::Coords coords(cpars, v0.GetCoordinates(b), k, j, i); const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); // y = U^(0) + dt S(y) // U^(0) values - const Real &dens = v0(b, gas::cons::density(), k, j, i); + const Real dens_raw = v0(b, gas::cons::density(), k, j, i); + const Real dens = + (std::isfinite(dens_raw) && dens_raw > dflr) ? dens_raw : dflr; Real Q = 0.0; - if (do_raytrace) Q = dt * v0(b, gas::src::energy(), k, j, i); - - // Note(AMD): There is some floating point difference between the internal energy - // used to compute the temperature and the internal energy obtained from that - // temperature: T = eos_d.TemperatureFromDensityInternalEnergy(dens, eg/dens); eg + if (do_raytrace) + Q = FiniteOrZero(dt * v0(b, gas::src::energy(), k, j, i)); + + // In the numerical atmosphere, retain thermal emission/absorption but + // suppress the momentum update. This avoids accelerating floor-density + // gas while still allowing absorbed raytraced energy to reradiate into + // the moment field through the conservative energy-only fallback. + constexpr Real atmosphere_floor_factor = 100.0; + const bool numerical_atmosphere = + !std::isfinite(dens_raw) || + dens_raw <= atmosphere_floor_factor * dflr; + + // Note(AMD): There is some floating point difference between the + // internal energy used to compute the temperature and the internal + // energy obtained from that temperature: T = + // eos_d.TemperatureFromDensityInternalEnergy(dens, eg/dens); eg // /= dens * eos_d.InternalEnergyFromDensityTemperature(dens,T) // - // Because of this, zero opacity problems will not result in zero change as - // expected. Thus, we recalculate the internal and total energies from the - // temperature. This does not affect energy conservation because at the end of the - // step we update the energy with an increment. - Real T = eos_d.TemperatureFromDensityInternalEnergy( - dens, v0(b, gas::cons::internal_energy(), k, j, i) / dens); + // Because of this, zero opacity problems will not result in zero change + // as expected. Thus, we recalculate the internal and total energies + // from the temperature. This does not affect energy conservation + // because at the end of the step we update the energy with an + // increment. + const Real eint0 = + std::max(0.0, v0(b, gas::cons::internal_energy(), k, j, i) / dens); + Real T = std::max( + tfloor, eos_d.TemperatureFromDensityInternalEnergy(dens, eint0)); Real eg0 = dens * eos_d.InternalEnergyFromDensityTemperature(dens, T); - Real B = arad * SQR(SQR(T)); + Real B = std::max(Bfloor_phys, arad * SQR(SQR(T))); const auto vb = RotatingFrame::BackgroundVelocity( qshear, om0, gm_bg, coords.GetCellCenter(vg, b, k, j, i)); @@ -294,172 +1209,765 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { vb[1] * dens + v0(b, gas::cons::momentum(1), k, j, i) / hx[1], vb[2] * dens + v0(b, gas::cons::momentum(2), k, j, i) / hx[2]}; - Real E0 = v0(b, rad::cons::energy(), k, j, i); - // choose the ref scale - - Real eref = std::sqrt(E0 * B); - if (eref == 0.0) eref = 0.5 * (E0 + B); + // The radiation-energy floor and material-temperature floor are + // independent constraints. In particular, an optically thin cell may + // have E_r << a_r T_floor^4. Promoting E_r to Bfloor_phys here creates + // an artificial LTE radiation bath and can make the constrained + // matter-coupling equations inconsistent at the first timestep. + const Real Er_state = + FiniteOrZero(v0(b, rad::cons::energy(), k, j, i)); + Real E0 = std::max(rad_efloor, Er_state); + // Use the arithmetic mean as the reference scale. Unlike the + // geometric mean, this keeps both normalized energies bounded when the + // gas and radiation temperatures are initially very different. + const Real eref_max = std::max(E0, B); + const Real eref_min = std::min(E0, B); + Real eref = + 0.5 * eref_max * (1.0 + eref_min / std::max(eref_max, Fuzz())); + if (!std::isfinite(eref) || eref <= 0.0) + eref = eref_max; + if (!std::isfinite(eref) || eref <= 0.0) { + v0(b, gas::cons::internal_energy(), k, j, i) += Q; + v0(b, gas::cons::total_energy(), k, j, i) += Q; + v0(b, rad::cons::energy(), k, j, i) = rad_efloor; + v0(b, rad::cons::flux(0), k, j, i) = 0.0; + v0(b, rad::cons::flux(1), k, j, i) = 0.0; + v0(b, rad::cons::flux(2), k, j, i) = 0.0; + return; + } const Real fref = c * eref; - const Real efloor = SQR(SQR(tfloor)) * arad / eref; - const Real Bfloor = efloor; + // Keep the radiation floor independent of the material temperature + // floor. Bfloor constrains T; efloor constrains E_r. + const Real efloor = rad_efloor / eref; + const Real Bfloor = Bfloor_phys / eref; Q /= eref; E0 /= eref; eg0 /= eref; B /= eref; - const std::array Fr0{v0(b, rad::cons::flux(0), k, j, i) / hx[0] / fref, - v0(b, rad::cons::flux(1), k, j, i) / hx[1] / fref, - v0(b, rad::cons::flux(2), k, j, i) / hx[2] / fref}; + const auto fred0 = NormalizeFlux( + FiniteOrZero(v0(b, rad::cons::flux(0), k, j, i) / + (hx[0] * fref * E0)), + FiniteOrZero(v0(b, rad::cons::flux(1), k, j, i) / + (hx[1] * fref * E0)), + FiniteOrZero(v0(b, rad::cons::flux(2), k, j, i) / + (hx[2] * fref * E0))); + const std::array Fr0{E0 * fred0[0], E0 * fred0[1], + E0 * fred0[2]}; std::array v{p0[0] / dens, p0[1] / dens, p0[2] / dens}; - const Real ke0 = 0.5 * dens * (SQR(v[0]) + SQR(v[1]) + SQR(v[2])) / eref; - const Real et0 = ke0 + eg0; + const std::array beta0{v[0] / c, v[1] / c, v[2] / c}; + const Real beta20 = SQR(beta0[0]) + SQR(beta0[1]) + SQR(beta0[2]); + if (!std::isfinite(beta20) || beta20 >= 1.0) { + v0(b, gas::cons::internal_energy(), k, j, i) += Q * eref; + v0(b, gas::cons::total_energy(), k, j, i) += Q * eref; + v0(b, rad::cons::energy(), k, j, i) = E0 * eref; + v0(b, rad::cons::flux(0), k, j, i) = Fr0[0] * hx[0] * fref; + v0(b, rad::cons::flux(1), k, j, i) = Fr0[1] * hx[1] * fref; + v0(b, rad::cons::flux(2), k, j, i) = Fr0[2] * hx[2] * fref; + return; + } + const Real ke0 = + 0.5 * dens * (SQR(v[0]) + SQR(v[1]) + SQR(v[2])) / eref; Real E = E0; auto F = Fr0; // Start outer iteration - int outer_iter = 0; int inner_iter = 0; - Real outer_err = 0.0; - Real inner_err = 0.0; + Real outer_err = std::numeric_limits::max(); + Real inner_err = std::numeric_limits::max(); - std::array dF{0., 0., 0.}; - std::array dv{0., 0., 0.}; + std::array dv{0.0, 0.0, 0.0}; Real dEk = 0.0; Real dEg = 0.0; Real dEr = 0.0; - const Real icc = 1. / (c * chat * dens); - Real escale = et0 + c / chat * E; - - for (outer_iter = 1; outer_iter <= outer_max; outer_iter++) { + const Real energy_floor_scale = + std::max(efloor + Bfloor, Fuzz()); + Real escale = std::max( + energy_floor_scale, + std::max(std::abs(eg0), + std::max(std::abs(Q), c / chat * std::abs(E0)))); + bool solve_valid = !numerical_atmosphere && std::isfinite(escale) && + std::isfinite(ke0) && escale > 0.0; + bool outer_converged = false; + bool have_complete_iterate = false; + Real best_outer_err = std::numeric_limits::max(); + std::array best_F = Fr0; + std::array best_dv{0.0, 0.0, 0.0}; + Real best_dEg = 0.0; + Real best_dEk = 0.0; + Real best_dEr = 0.0; + bool scalar_inner_used = false; + int momentum_predictor_count = 0; + int momentum_iteration_count = 0; + Real momentum_error = std::numeric_limits::max(); + + for (outer_iter = 1; outer_iter <= outer_max; ++outer_iter) { + if (!solve_valid) + break; - // Set some v and F quantities - Real ke = 0.5 * dens * (SQR(v[0]) + SQR(v[1]) + SQR(v[2])) / eref; + const Real ke = + 0.5 * dens * (SQR(v[0]) + SQR(v[1]) + SQR(v[2])) / eref; std::array beta{v[0] / c, v[1] / c, v[2] / c}; const Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); - const Real g2 = 1. / (1. - beta2); + if (!std::isfinite(beta2) || beta2 >= 1.0 || E <= 0.0) { + solve_valid = false; + break; + } + const Real g2 = 1.0 / (1.0 - beta2); const Real g = std::sqrt(g2); - auto fedd = - EddingtonTensor({F[0] / (c * E), F[1] / (c * E), F[2] / (c * E)}); - - std::array bdp{ + // F is normalized by c*eref, so F/E is the reduced flux. + const auto fedd = + EddingtonTensor({F[0] / E, F[1] / E, F[2] / E}); + const std::array bdp{ beta[0] * fedd[TensIdx::X11] + beta[1] * fedd[TensIdx::X12] + beta[2] * fedd[TensIdx::X13], beta[0] * fedd[TensIdx::X12] + beta[1] * fedd[TensIdx::X22] + beta[2] * fedd[TensIdx::X23], beta[0] * fedd[TensIdx::X13] + beta[1] * fedd[TensIdx::X23] + beta[2] * fedd[TensIdx::X33]}; - const Real bdbdp = beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; - const Real bdf = beta[0] * F[0] + beta[1] * F[1] + beta[2] * F[2]; // 1/c - - // start inner iteration for (B,E) - for (inner_iter = 1; inner_iter <= inner_max; inner_iter++) { - T = std::pow(eref * B / arad, 0.25); - Real eint = dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / eref; - Real et = ke + eint; - const Real Cv = dens * eos_d.SpecificHeatFromDensityTemperature(dens, T); + const Real bdbdp = + beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; + const Real bdf = + beta[0] * F[0] + beta[1] * F[1] + beta[2] * F[2]; + + // Damped quasi-Newton solve for (B,E). Opacity derivatives are not + // available, so accept only residual-decreasing trial steps. + Real previous_inner_err = std::numeric_limits::max(); + int stalled_iterations = 0; + bool inner_converged = false; + for (inner_iter = 1; inner_iter <= inner_max; ++inner_iter) { + T = std::max( + tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); + const Real eint = + dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / + eref; + const Real Cv = + dens * eos_d.SpecificHeatFromDensityTemperature(dens, T); const Real fleck = FleckFactor(arad, T, Cv); - const Real sigp = chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); + const Real sigp = + chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); const Real sigs = - chat * dt * scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); - const Real sigf = sigp + sigs; - - const Real ca = g * (sigf - g2 * sigs * (1. + bdbdp)); - const Real cb = g * sigp; - const Real cd = -g * bdf * (sigf - 2. * g2 * sigs); - + chat * dt * + scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); + Real ca = 0.0; + Real cb = 0.0; + Real cd = 0.0; + const bool coeff_valid = ComputeCouplingEnergyCoefficients( + sigp, sigs, g, g2, beta2, bdbdp, bdf, ca, cb, cd); + Real Eeq_inner = E; + const bool equilibrium_valid = + coeff_valid && ComputeCouplingEquilibriumEnergy( + E0, B, ca, cb, cd, Eeq_inner); + + // Retain the original residuals only to form the quasi-Newton + // search direction. Use the conservation/equilibrium form below + // for convergence and line-search acceptance. const Real G0 = ca * E - cb * B + cd; - const Real Fi = (et - et0) - c / chat * G0 - Q; + const Real Fi = + (ke - ke0) + (eint - eg0) - c / chat * G0 - Q; const Real Fr = (E - E0) + G0; - - // not converged yet - const Real dfac = 1. + c / chat * fleck * cb; - Real dE = dfac / (dfac + ca) * (-Fr) + fleck / (dfac + ca) * (-Fi * cb); - Real dB = c / chat * fleck / (dfac + ca) * (-ca * Fr) + - (1. + ca) * fleck / (dfac + ca) * (-Fi); - Real Enew = E + dE; - E = (Enew < efloor) ? efloor : Enew; - Real Bnew = B + dB; - B = (Bnew < Bfloor) ? Bfloor : Bnew; - + const Real conservation_inner = + (ke - ke0) + (eint - eg0) + c / chat * (E - E0) - Q; + const Real source_inner = E - Eeq_inner; + + const Real escale_inner = std::max( + energy_floor_scale, + std::max( + std::max(std::abs(eg0), std::abs(eint)), + std::max(std::abs(Q), + std::max(c / chat * + (std::abs(E0) + std::abs(E)), + std::abs(ke - ke0))))); inner_err = - std::max((std::abs(Fi) / escale), (c / chat * std::abs(Fr) / escale)); + std::max(std::abs(conservation_inner) / escale_inner, + c / chat * std::abs(source_inner) / escale_inner); + solve_valid = + solve_valid && equilibrium_valid && + std::isfinite(inner_err) && std::isfinite(fleck) && + std::isfinite(eint) && Cv > 0.0 && fleck >= 0.0 && + sigp >= 0.0 && sigs >= 0.0; + if (!solve_valid) + break; if (inner_err <= inner_tol) { - // converged, so don't compute new E and B; + inner_converged = true; + break; + } + + if (inner_err >= previous_inner_err * (1.0 - 1.0e-6)) + ++stalled_iterations; + else + stalled_iterations = 0; + previous_inner_err = inner_err; + if (stalled_iterations >= 8) break; + + const Real dfac = 1.0 + c / chat * fleck * cb; + const Real denom = dfac + ca; + if (!std::isfinite(denom) || std::abs(denom) <= Fuzz()) { + solve_valid = false; + break; + } + const Real dE = + dfac / denom * (-Fr) + fleck / denom * (-Fi * cb); + const Real dB = c / chat * fleck / denom * (-ca * Fr) + + (1.0 + ca) * fleck / denom * (-Fi); + + bool accepted = false; + Real best_trial_err = inner_err; + Real best_trial_E = E; + Real best_trial_B = B; + Real alpha = 1.0; + for (int ls = 0; ls < 12; ++ls) { + const Real Etrial = std::max(efloor, E + alpha * dE); + const Real Btrial = std::max(Bfloor, B + alpha * dB); + const Real Ttrial = std::max( + tfloor, + std::pow(std::max(eref * Btrial / arad, 0.0), 0.25)); + const Real eint_trial = + dens * + eos_d.InternalEnergyFromDensityTemperature(dens, Ttrial) / + eref; + const Real sigp_trial = + chat * dt * + opac_d.PlanckMeanAbsorptionCoefficient(dens, Ttrial); + const Real sigs_trial = + chat * dt * + scat_d.RosselandMeanTotalScatteringCoefficient(dens, Ttrial); + Real ca_trial = 0.0; + Real cb_trial = 0.0; + Real cd_trial = 0.0; + const bool coeff_trial_valid = + ComputeCouplingEnergyCoefficients( + sigp_trial, sigs_trial, g, g2, beta2, bdbdp, bdf, + ca_trial, cb_trial, cd_trial); + Real Eeq_trial = Etrial; + const bool equilibrium_trial_valid = + coeff_trial_valid && ComputeCouplingEquilibriumEnergy( + E0, Btrial, ca_trial, cb_trial, + cd_trial, Eeq_trial); + const Real conservation_trial = + (ke - ke0) + (eint_trial - eg0) + + c / chat * (Etrial - E0) - Q; + const Real source_trial = Etrial - Eeq_trial; + const Real escale_trial = std::max( + energy_floor_scale, + std::max( + std::max(std::abs(eg0), std::abs(eint_trial)), + std::max( + std::abs(Q), + std::max(c / chat * + (std::abs(E0) + std::abs(Etrial)), + std::abs(ke - ke0))))); + const Real trial_err = + std::max(std::abs(conservation_trial) / escale_trial, + c / chat * std::abs(source_trial) / escale_trial); + const bool trial_valid = + std::isfinite(Ttrial) && std::isfinite(eint_trial) && + std::isfinite(sigp_trial) && std::isfinite(sigs_trial) && + equilibrium_trial_valid && std::isfinite(trial_err) && + sigp_trial >= 0.0 && sigs_trial >= 0.0; + if (trial_valid && trial_err < best_trial_err) { + best_trial_err = trial_err; + best_trial_E = Etrial; + best_trial_B = Btrial; + } + if (trial_valid && + (trial_err <= inner_tol || + trial_err <= inner_err * (1.0 - 1.0e-4 * alpha))) { + E = Etrial; + B = Btrial; + accepted = true; + break; + } + alpha *= 0.5; } + if (!accepted) { + if (best_trial_err < inner_err) { + E = best_trial_E; + B = best_trial_B; + } else { + break; + } + } + } // inner iteration - } // inner_iter - if ((inner_iter > inner_max) && (fatal_if_unconverged)) { - printf("(%d,%d,%d,%d) %lg > %lg after %d iterations\n", b, k, j, i, - inner_err, inner_tol, inner_max); - PARTHENON_FAIL("Inner not converged"); + if (!inner_converged && solve_valid && + (inner_iter > inner_max || stalled_iterations >= 8) && + inner_err <= std::max(inner_tol, nonlinear_roundoff_tol)) { + inner_converged = true; } - // Have new E and T + // The opacity-frozen quasi-Newton Jacobian can stall when opacity is + // strongly temperature dependent or the solution lies on a floor. + // In that case eliminate E using exact total-energy conservation and + // solve the remaining full mixed-frame radiation equation as a + // bracketed scalar problem in gas internal energy. The scalar + // residual compares the energy from exact conservation with the + // energy from the implicit source equation, avoiding opacity- + // amplified cancellation near LTE. + if (!inner_converged && solve_valid) { + scalar_inner_used = true; + int scalar_iterations = 0; + const bool scalar_converged = SolveCouplingInnerScalar( + eg0, ke - ke0, Q, E0, efloor, Bfloor, dens, eref, arad, + tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, + energy_floor_scale, inner_tol, eos_d, opac_d, scat_d, E, B, + inner_err, scalar_iterations); + inner_iter += scalar_iterations; + inner_converged = scalar_converged; + } + if (!inner_converged && solve_valid && + momentum_predictor_count < outer_max) { + // Bootstrap the kinetic-work term with the same coupled momentum + // solve used below. Unlike the old full Picard predictor, this + // cannot jump to a superluminal intermediate state when the + // reduced-c radiation pseudo-inertia dominates the gas inertia. + const Real Tpred = std::max( + tfloor, + std::pow(std::max(eref * B / arad, 0.0), 0.25)); + const Real eg_pred = + dens * eos_d.InternalEnergyFromDensityTemperature(dens, + Tpred) / + eref; + const Real dEg_pred = eg_pred - eg0; + const Real sigp_pred = + chat * dt * + opac_d.RosselandMeanAbsorptionCoefficient(dens, Tpred); + const Real sigs_pred = + chat * dt * + scat_d.RosselandMeanTotalScatteringCoefficient(dens, Tpred); + std::array beta_pred{v[0] / c, v[1] / c, v[2] / c}; + std::array Fpred = F; + Real Epred = E; + Real dEkpred = dEk; + Real predictor_error = std::numeric_limits::max(); + int predictor_iterations = 0; + const bool predictor_converged = + SolveCouplingMomentumEnergy( + beta0, Fr0, E0, efloor, B, dEg_pred, Q, dens, eref, + c, chat, sigp_pred, sigs_pred, energy_floor_scale, + std::max(outer_tol, nonlinear_roundoff_tol), beta_pred, + Fpred, Epred, dEkpred, predictor_error, + predictor_iterations); + momentum_iteration_count += predictor_iterations; + if (predictor_converged) { + const Real predictor_scale = std::max( + energy_floor_scale, + std::max(E, std::max(CouplingFluxMagnitude(F), + CouplingFluxMagnitude(Fpred)))); + Real predictor_change = std::abs(Epred - E) / predictor_scale; + for (int d = 0; d < 3; ++d) { + predictor_change = std::max( + predictor_change, + std::abs(Fpred[d] - F[d]) / predictor_scale); + predictor_change = std::max( + predictor_change, + std::abs(beta_pred[d] - v[d] / c)); + } + if (predictor_change > + 8.0 * std::numeric_limits::epsilon()) { + E = Epred; + F = Fpred; + dEk = dEkpred; + const Real dv_factor = eref / (dens * chat); + for (int d = 0; d < 3; ++d) { + dv[d] = -dv_factor * (F[d] - Fr0[d]); + v[d] = p0[d] / dens + dv[d]; + } + ++momentum_predictor_count; + --outer_iter; + continue; + } + } + } + if (!inner_converged) + break; - T = std::pow(eref * B / arad, 0.25); - Real eg = dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / eref; + T = std::max( + tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); + const Real eg = + dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / + eref; dEg = eg - eg0; - const Real sigp = - chat * dt * opac_d.RosselandMeanAbsorptionCoefficient(dens, T); - const Real sigs = - chat * dt * scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); - const Real sigf = sigp + sigs; - - const Real a = g * sigf; - const Real b = 2. * g2 * g * sigs; - const Real d1 = g * (sigp * B + g2 * sigs * (1. + bdbdp) * E); // * c - const Real d2 = g * sigf * E; // * c - const std::array rhs{Fr0[0] + d1 * beta[0] + d2 * bdp[0], - Fr0[1] + d1 * beta[1] + d2 * bdp[1], - Fr0[2] + d1 * beta[2] + d2 * bdp[2]}; - - F = SolveRadFlux(1. + a, b, beta, rhs); - - for (int d = 0; d < 3; d++) { - dF[d] = F[d] - Fr0[d]; - dv[d] = -icc * dF[d] * fref; - v[d] = p0[d] / dens + dv[d]; + const Real sigp_flux = + chat * dt * + opac_d.RosselandMeanAbsorptionCoefficient(dens, T); + const Real sigs_flux = + chat * dt * + scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); + if (!std::isfinite(sigp_flux) || !std::isfinite(sigs_flux) || + sigp_flux < 0.0 || sigs_flux < 0.0) { + solve_valid = false; + break; } - const Real dEk_prev = dEk; - dEk = 0.5 * dens * - (dv[0] * (v[0] + p0[0] / dens) + dv[1] * (v[1] + p0[1] / dens) + - dv[2] * (v[2] + p0[2] / dens)) / - eref; - dEr = -chat / c * (dEg + dEk); - E = E0 + dEr; + // Solve the radiation-flux source equation together with gas + // pseudo-momentum conservation and reduced-c total-energy + // conservation. A direct Picard update is unstable when + // mu = eref/(rho*c*chat) + // is large; that is exactly the regime encountered just above a + // density floor. The local 3D Newton solve removes this stiffness. + std::array beta_momentum{v[0] / c, v[1] / c, + v[2] / c}; + int momentum_iterations = 0; + const bool momentum_converged = + SolveCouplingMomentumEnergy( + beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, + chat, sigp_flux, sigs_flux, energy_floor_scale, + std::max(outer_tol, nonlinear_roundoff_tol), + beta_momentum, F, E, dEk, momentum_error, + momentum_iterations); + momentum_iteration_count += momentum_iterations; + if (!momentum_converged) { + solve_valid = false; + break; + } + const Real dv_factor = eref / (dens * chat); + for (int d = 0; d < 3; ++d) { + dv[d] = -dv_factor * (F[d] - Fr0[d]); + v[d] = p0[d] / dens + dv[d]; + } + dEr = E - E0; + + // Recompute the full coupled residual from the final state of this + // outer iteration. This avoids declaring convergence merely because + // the kinetic-energy correction changed little in a Keplerian flow. + const std::array beta_out{v[0] / c, v[1] / c, v[2] / c}; + const Real beta2_out = SQR(beta_out[0]) + SQR(beta_out[1]) + + SQR(beta_out[2]); + if (!std::isfinite(beta2_out) || beta2_out >= 1.0 || E <= 0.0) { + solve_valid = false; + break; + } + const Real g2_out = 1.0 / (1.0 - beta2_out); + const Real g_out = std::sqrt(g2_out); + const auto fedd_out = EddingtonTensor( + {F[0] / E, F[1] / E, F[2] / E}); + const std::array bdp_out{ + beta_out[0] * fedd_out[TensIdx::X11] + + beta_out[1] * fedd_out[TensIdx::X12] + + beta_out[2] * fedd_out[TensIdx::X13], + beta_out[0] * fedd_out[TensIdx::X12] + + beta_out[1] * fedd_out[TensIdx::X22] + + beta_out[2] * fedd_out[TensIdx::X23], + beta_out[0] * fedd_out[TensIdx::X13] + + beta_out[1] * fedd_out[TensIdx::X23] + + beta_out[2] * fedd_out[TensIdx::X33]}; + const Real bdbdp_out = beta_out[0] * bdp_out[0] + + beta_out[1] * bdp_out[1] + + beta_out[2] * bdp_out[2]; + const Real bdf_out = beta_out[0] * F[0] + beta_out[1] * F[1] + + beta_out[2] * F[2]; + + T = std::max( + tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); + const Real eg_out = + dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / + eref; + dEg = eg_out - eg0; + const Real sigp_energy = + chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); + const Real sigs_out = + chat * dt * + scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); + Real ca_out = 0.0; + Real cb_out = 0.0; + Real cd_out = 0.0; + Real Eeq_out = E; + const bool energy_coeff_valid = ComputeCouplingEnergyCoefficients( + sigp_energy, sigs_out, g_out, g2_out, beta2_out, + bdbdp_out, bdf_out, ca_out, cb_out, cd_out); + const bool equilibrium_valid = + energy_coeff_valid && ComputeCouplingEquilibriumEnergy( + E0, B, ca_out, cb_out, cd_out, + Eeq_out); + const Real conservation_residual = + dEg + dEk + c / chat * dEr - Q; + const Real source_residual = E - Eeq_out; + + const Real sigp_flux_out = + chat * dt * + opac_d.RosselandMeanAbsorptionCoefficient(dens, T); + const Real sigf_flux_out = sigp_flux_out + sigs_out; + const Real a_out = g_out * sigf_flux_out; + const Real b_out = 2.0 * g2_out * g_out * sigs_out; + const Real d1_out = + g_out * + (sigp_flux_out * B + + g2_out * sigs_out * (1.0 + bdbdp_out) * E); + const Real d2_out = g_out * sigf_flux_out * E; + const std::array rhs_out{ + Fr0[0] + d1_out * beta_out[0] + d2_out * bdp_out[0], + Fr0[1] + d1_out * beta_out[1] + d2_out * bdp_out[1], + Fr0[2] + d1_out * beta_out[2] + d2_out * bdp_out[2]}; + const auto Fsolve_out = + SolveRadFlux(1.0 + a_out, b_out, beta_out, rhs_out); + const auto Ftarget_out = ProjectCouplingFluxToEnergy(Fsolve_out, E); + + escale = std::max( + energy_floor_scale, + std::max( + std::max(std::abs(eg0), std::abs(eg_out)), + std::max( + std::abs(Q), + std::max(c / chat * + (std::abs(E0) + std::abs(E)), + std::abs(dEk))))); + const Real energy_residual = + std::max(std::abs(conservation_residual) / escale, + c / chat * std::abs(source_residual) / escale); + const Real flux_scale = std::max( + energy_floor_scale, + std::max( + E, std::max(CouplingFluxMagnitude(Fr0), + std::max(CouplingFluxMagnitude(F), + CouplingFluxMagnitude(Ftarget_out))))); + Real flux_residual = 0.0; + for (int d = 0; d < 3; ++d) { + flux_residual = + std::max(flux_residual, + std::abs(F[d] - Ftarget_out[d]) / flux_scale); + } + const Real realizability_residual = + std::max(0.0, (CouplingFluxMagnitude(F) - E) / + std::max(E, energy_floor_scale)); + // These are residuals of the coupled equations evaluated at the + // current state. Iterate-to-iterate changes are deliberately not an + // acceptance criterion: a converged nonlinear root need not move by + // less than a user tolerance tighter than the local arithmetic can + // resolve. + outer_err = std::max( + energy_residual, + std::max(flux_residual, realizability_residual)); + solve_valid = + solve_valid && equilibrium_valid && + std::isfinite(outer_err) && std::isfinite(dEr) && + std::isfinite(sigp_energy) && std::isfinite(sigs_out) && + std::isfinite(sigp_flux_out) && + sigp_energy >= 0.0 && sigs_out >= 0.0 && + sigp_flux_out >= 0.0; + if (!solve_valid) + break; - // This needs to be something else related to the change from the last iteration - outer_err = std::abs(dEk - dEk_prev) / escale; + have_complete_iterate = true; + if (outer_err < best_outer_err) { + best_outer_err = outer_err; + best_F = F; + best_dv = dv; + best_dEg = dEg; + best_dEk = dEk; + best_dEr = dEr; + } if (outer_err <= outer_tol) { + outer_converged = true; break; } + } // outer iteration + + // A hard floor turns the nonlinear equality problem into a + // constrained active-set problem. If the unconstrained root lies below + // either floor, a zero equality residual does not exist in the + // admissible state space. In that case use the conservative + // energy-only solve below with momentum frozen, rather than repeatedly + // applying a radiation drag whose kinetic-energy work cannot be paid by + // either the gas or radiation field. + const Real floor_active_tol = + 256.0 * std::numeric_limits::epsilon(); + const bool material_floor_active = + B <= Bfloor * (1.0 + floor_active_tol) + Fuzz(); + const bool radiation_floor_active = + E <= efloor * (1.0 + floor_active_tol) + Fuzz(); + const bool floor_constrained_failure = + solve_valid && !have_complete_iterate && + (material_floor_active || radiation_floor_active); + bool fallback_success = false; + Real fallback_error = std::numeric_limits::max(); + + if (have_complete_iterate) { + F = best_F; + dv = best_dv; + for (int d = 0; d < 3; ++d) + v[d] = p0[d] / dens + dv[d]; + dEg = best_dEg; + dEk = best_dEk; + dEr = best_dEr; + E = E0 + dEr; + } else { + // Conservative energy-only fallback. Hold radiation flux and gas + // momentum fixed, eliminate E using total-energy conservation, and + // bracket the remaining scalar radiation-energy residual in gas + // internal energy. This is also the intentional atmosphere path. + F = Fr0; + dv = {0.0, 0.0, 0.0}; + v = {p0[0] / dens, p0[1] / dens, p0[2] / dens}; + dEk = 0.0; + const Real Emin = std::max(efloor, CouplingFluxMagnitude(Fr0)); + const Real eg_floor = + dens * eos_d.InternalEnergyFromDensityTemperature(dens, + tfloor) / + eref; + const Real eg_max = + eg0 + Q + c / chat * (E0 - Emin); + + bool fallback_valid = std::isfinite(eg_floor) && + std::isfinite(eg_max) && eg_max >= eg_floor; + bool fallback_have = false; + bool bracketed = false; + Real bracket_lo = eg_floor; + Real bracket_hi = eg_max; + Real residual_lo = 0.0; + Real best_eg = eg0; + Real best_E_fallback = E0; + Real best_B_fallback = B; + Real best_residual = std::numeric_limits::max(); + Real previous_eg = eg_floor; + Real previous_residual = 0.0; + bool have_previous = false; + + if (fallback_valid) { + constexpr int fallback_scan_points = 64; + for (int n = 0; n <= fallback_scan_points; ++n) { + const Real frac = + static_cast(n) / fallback_scan_points; + const Real eg_trial = + eg_floor + frac * (eg_max - eg_floor); + Real E_trial = E0; + Real B_trial = B; + Real residual = 0.0; + const bool valid_trial = EvaluateCouplingEnergyOnlyResidual( + eg_trial, eg0, Q, E0, Emin, dens, eref, arad, tfloor, dt, + chat, c, eos_d, opac_d, E_trial, B_trial, residual); + if (!valid_trial) + continue; + + fallback_have = true; + if (std::abs(residual) < best_residual) { + best_residual = std::abs(residual); + best_eg = eg_trial; + best_E_fallback = E_trial; + best_B_fallback = B_trial; + } + if (have_previous && !bracketed && + ((previous_residual <= 0.0 && residual >= 0.0) || + (previous_residual >= 0.0 && residual <= 0.0))) { + bracketed = true; + bracket_lo = previous_eg; + bracket_hi = eg_trial; + residual_lo = previous_residual; + } + previous_eg = eg_trial; + previous_residual = residual; + have_previous = true; + } + } + + if (bracketed) { + for (int n = 0; n < 64; ++n) { + const Real eg_mid = 0.5 * (bracket_lo + bracket_hi); + Real E_mid = E0; + Real B_mid = B; + Real residual_mid = 0.0; + const bool valid_mid = EvaluateCouplingEnergyOnlyResidual( + eg_mid, eg0, Q, E0, Emin, dens, eref, arad, tfloor, dt, + chat, c, eos_d, opac_d, E_mid, B_mid, residual_mid); + if (!valid_mid) + break; + if (std::abs(residual_mid) < best_residual) { + best_residual = std::abs(residual_mid); + best_eg = eg_mid; + best_E_fallback = E_mid; + best_B_fallback = B_mid; + } + if ((residual_lo <= 0.0 && residual_mid >= 0.0) || + (residual_lo >= 0.0 && residual_mid <= 0.0)) { + bracket_hi = eg_mid; + } else { + bracket_lo = eg_mid; + residual_lo = residual_mid; + } + if ((bracket_hi - bracket_lo) / + std::max(energy_floor_scale, + std::abs(bracket_hi) + + std::abs(bracket_lo)) <= + nonlinear_roundoff_tol) + break; + } + } + + if (fallback_have) { + dEg = best_eg - eg0; + E = best_E_fallback; + B = best_B_fallback; + dEr = E - E0; + const Real fallback_scale = std::max( + energy_floor_scale, + std::max(std::abs(Q), + std::max(std::abs(eg0) + std::abs(best_eg), + c / chat * + (std::abs(E0) + std::abs(E))))); + fallback_error = + c / chat * best_residual / std::max(fallback_scale, Fuzz()); + fallback_success = std::isfinite(fallback_error); + } else { + // Last-resort finite update. This is intentionally limited to the + // external ray source and leaves the moment state unchanged. + dEg = Q; + dEr = 0.0; + E = E0; + F = Fr0; + fallback_success = false; + } + } - } // outer_iter - if ((outer_iter > outer_max) && (fatal_if_unconverged)) { - printf("(%d,%d,%d,%d) %lg > %lg after %d iterations\n", b, k, j, i, outer_err, - outer_tol, outer_max); + const bool accepted_floor_fallback = + floor_constrained_failure && fallback_success; + if (!outer_converged && fatal_if_unconverged && + !numerical_atmosphere && !accepted_floor_fallback) { + const Real beta_fail = + std::sqrt(SQR(v[0] / c) + SQR(v[1] / c) + SQR(v[2] / c)); + const Real fred_fail = + CouplingFluxMagnitude(F) / std::max(E, energy_floor_scale); + const Real mu_fail = eref / (dens * c * chat); + const Real delta_f_fail = CouplingFluxMagnitude( + {F[0] - Fr0[0], F[1] - Fr0[1], F[2] - Fr0[2]}); + printf("MatterCoupling full fail (%d,%d,%d,%d): outer=%.17e " + "(best=%.17e, tol=%.17e), inner=%.17e (tol=%.17e), " + "fallback=%.17e, E=%.17e, E0=%.17e, efloor=%.17e, " + "B=%.17e, Bfloor=%.17e, Q=%.17e, dEk=%.17e, " + "rho=%.17e, rho/dfloor=%.17e, beta0=%.17e, " + "beta=%.17e, fred=%.17e, dF=%.17e, chat/c=%.17e, " + "mu=%.17e, flux_variable=%d, outer_iter=%d, " + "inner_iter=%d, momentum_iter=%d, " + "momentum_err=%.17e, valid=%d, scalar=%d, " + "predictors=%d, complete=%d, floor_constrained=%d\n", + b, k, j, i, outer_err, best_outer_err, outer_tol, inner_err, + inner_tol, fallback_error, E, E0, efloor, B, Bfloor, Q, + dEk, dens, dens / dflr, std::sqrt(beta20), beta_fail, + fred_fail, delta_f_fail, chat / c, mu_fail, + static_cast(mu_fail <= 1.0), outer_iter, inner_iter, + momentum_iteration_count, momentum_error, + solve_valid, scalar_inner_used, momentum_predictor_count, + have_complete_iterate, floor_constrained_failure); PARTHENON_FAIL("Outer not converged"); } // Update state vector (both gas and radiation) v0(b, gas::cons::internal_energy(), k, j, i) += dEg * eref; v0(b, gas::cons::total_energy(), k, j, i) += (dEg + dEk) * eref; - v0(b, rad::cons::energy(), k, j, i) += dEr * eref; + v0(b, rad::cons::energy(), k, j, i) = (E0 + dEr) * eref; v0(b, gas::cons::momentum(0), k, j, i) += dv[0] * dens * hx[0]; v0(b, gas::cons::momentum(1), k, j, i) += dv[1] * dens * hx[1]; v0(b, gas::cons::momentum(2), k, j, i) += dv[2] * dens * hx[2]; - v0(b, rad::cons::flux(0), k, j, i) += dF[0] * hx[0] * fref; - v0(b, rad::cons::flux(1), k, j, i) += dF[1] * hx[1] * fref; - v0(b, rad::cons::flux(2), k, j, i) += dF[2] * hx[2] * fref; + v0(b, rad::cons::flux(0), k, j, i) = F[0] * hx[0] * fref; + v0(b, rad::cons::flux(1), k, j, i) = F[1] * hx[1] * fref; + v0(b, rad::cons::flux(2), k, j, i) = F[2] * hx[2] * fref; }); return TaskStatus::complete; From efb60a21328b6a48f1c287587bbf84f485727fbd Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Mon, 13 Jul 2026 07:37:12 -0600 Subject: [PATCH 02/14] Add sub cycling as an option --- src/radiation/moments/moments.cpp | 26 +++++++++++++++-- src/radiation/moments/moments.hpp | 46 +++++++++++++++++++++++++++++++ 2 files changed, 69 insertions(+), 3 deletions(-) diff --git a/src/radiation/moments/moments.cpp b/src/radiation/moments/moments.cpp index 4cb1ca85..e65cd8d2 100644 --- a/src/radiation/moments/moments.cpp +++ b/src/radiation/moments/moments.cpp @@ -18,6 +18,7 @@ #include "artemis.hpp" #include "geometry/geometry.hpp" #include "matter_coupling.hpp" +#include "matter_coupling_simple.hpp" #include "moments.hpp" #include "utils/artemis_utils.hpp" #include "utils/fluxes/fluid_fluxes.hpp" @@ -109,7 +110,7 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("tfloor", tfloor * units.GetTemperaturePhysicalToCode()); params.Add("use_opac", - pin->GetOrAddBoolean("radiation/moments", "init_with_opac", true)); + pin->GetOrAddBoolean("radiation/moment", "init_with_opac", true)); // Number of radiation species const int nspecies = pin->GetOrAddInteger("radiation/moment", "nspecies", 1); @@ -126,6 +127,25 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("inner_iteration_tol", pin->GetOrAddReal("radiation/moment", "inner_iteration_tol", 1e-10)); + const bool substep = pin->GetOrAddReal("radiation/moment", "substep", true); + if (!substep) { + if (coords == Coordinates::cartesian) { + moments->EstimateTimestepMesh = EstimateTimeStepMesh; + } else if (coords == Coordinates::spherical1D) { + moments->EstimateTimestepMesh = EstimateTimeStepMesh; + } else if (coords == Coordinates::spherical2D) { + moments->EstimateTimestepMesh = EstimateTimeStepMesh; + } else if (coords == Coordinates::spherical3D) { + moments->EstimateTimestepMesh = EstimateTimeStepMesh; + } else if (coords == Coordinates::cylindrical) { + moments->EstimateTimestepMesh = EstimateTimeStepMesh; + } else if (coords == Coordinates::axisymmetric) { + moments->EstimateTimestepMesh = EstimateTimeStepMesh; + } else { + PARTHENON_FAIL("Invalid artemis/coordinate system!"); + } + } + // Number of radiation "species" (i.e., groups) std::vector fluidids; for (int n = 0; n < nspecies; ++n) @@ -369,13 +389,13 @@ TaskStatus MatterCoupling(MeshData *u0, const Real dt) { if (full_coupling) { return MatterCouplingFullSingleImpl(u0, dt); } else { - return MatterCouplingSimpleImpl(u0, dt); + // return MatterCouplingSimpleImpl(u0, dt); } } else if (closure_type == Closure::p1) { if (full_coupling) { return MatterCouplingFullSingleImpl(u0, dt); } else { - return MatterCouplingSimpleImpl(u0, dt); + // return MatterCouplingSimpleImpl(u0, dt); } } return TaskStatus::complete; diff --git a/src/radiation/moments/moments.hpp b/src/radiation/moments/moments.hpp index 251b7d72..c69ab320 100644 --- a/src/radiation/moments/moments.hpp +++ b/src/radiation/moments/moments.hpp @@ -43,6 +43,52 @@ TaskCollection MomentsTasks(Mesh *pmesh, const SimTime &tm, template void InitMesh(parthenon::Mesh *pmesh); //---------------------------------------------------------------------------------------- +//! \fn Real Moments::EstimateTimeStepMesh +//! \brief Not enrolled in parthenon's determination for global dt +template +Real EstimateTimeStepMesh(MeshData *md) { + PARTHENON_INSTRUMENT + auto pm = md->GetParentPointer(); + auto &resolved_pkgs = pm->resolved_packages; + auto &moments_pkg = pm->packages.Get("moments"); + auto ¶ms = moments_pkg->AllParams(); + + Real dxmin = Big(); + + // Packing and Indexing + static auto desc_g = + MakePackDescriptor(resolved_pkgs.get()); + auto vg = desc_g.GetPack(md); + + IndexRange ib = md->GetBoundsI(IndexDomain::interior); + IndexRange jb = md->GetBoundsJ(IndexDomain::interior); + IndexRange kb = md->GetBoundsK(IndexDomain::interior); + const auto ndim = pm->ndim; + const auto &cpars = + pm->packages.Get("artemis")->template Param("coord_params"); + + // Compute minimum dx + Real min_dx = Big(); + parthenon::par_reduce( + parthenon::loop_pattern_mdrange_tag, "Moments::EstimateTimestepMesh", + DevExecSpace(), 0, md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &ldx_m) { + // Extract coordinates + geometry::Coords coords(cpars, vg.GetCoordinates(b), k, j, i); + const auto &dx = coords.GetCellWidths(vg, b, k, j, i); + for (int d = 0; d < ndim; d++) { + ldx_m = std::min(ldx_m, dx[d]); + } + }, + Kokkos::Min(min_dx)); + + dxmin = std::min(dxmin, min_dx); + + const auto chat = params.template Get("chat"); + const auto cfl = params.template Get("cfl"); + return cfl * dxmin / chat; +} +//---------------------------------------------------------------------------------------- //! \fn Real Moments::EstimateTimeStep //! \brief Not enrolled in parthenon's determination for global dt template From 6a8310174d341b3d913ac66807185ecdacc9c473 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Mon, 13 Jul 2026 07:38:17 -0600 Subject: [PATCH 03/14] substep is a bool not real --- src/radiation/moments/moments.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/radiation/moments/moments.cpp b/src/radiation/moments/moments.cpp index e65cd8d2..1e7d7a66 100644 --- a/src/radiation/moments/moments.cpp +++ b/src/radiation/moments/moments.cpp @@ -127,7 +127,7 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("inner_iteration_tol", pin->GetOrAddReal("radiation/moment", "inner_iteration_tol", 1e-10)); - const bool substep = pin->GetOrAddReal("radiation/moment", "substep", true); + const bool substep = pin->GetOrAddBoolean("radiation/moment", "substep", true); if (!substep) { if (coords == Coordinates::cartesian) { moments->EstimateTimestepMesh = EstimateTimeStepMesh; From 71cd205f70275bc53f1ed4e333ccf248aa7e5b77 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Mon, 13 Jul 2026 07:39:01 -0600 Subject: [PATCH 04/14] Be careful of floors and tolerances --- src/radiation/raytrace/raytrace.cpp | 24 +++++++++++++++++------- src/radiation/raytrace/raytrace.hpp | 27 ++++++++++++++++++++------- 2 files changed, 37 insertions(+), 14 deletions(-) diff --git a/src/radiation/raytrace/raytrace.cpp b/src/radiation/raytrace/raytrace.cpp index 3f888bba..0ee3502b 100644 --- a/src/radiation/raytrace/raytrace.cpp +++ b/src/radiation/raytrace/raytrace.cpp @@ -39,6 +39,7 @@ std::shared_ptr Initialize(ParameterInput *pin, ArtemisUtils::Opacity opacity; if (opacity_model_name == "constant") { const Real kappa_a = pin->GetOrAddReal(block_name, "kappa_a", 0.0); + PARTHENON_REQUIRE(kappa_a >= 0.0, "Raytrace absorption opacity must be non-negative"); opacity = NonCGSUnits(Gray(kappa_a), time, mass, length, temp); } else { PARTHENON_FAIL("Opacity model not recognized!"); @@ -48,6 +49,8 @@ std::shared_ptr Initialize(ParameterInput *pin, const Real stellar_temp = pin->GetReal("radiation/raytrace", "temperature_cgs"); Real stellar_radius = pin->GetReal("radiation/raytrace", "radius_solar"); + PARTHENON_REQUIRE(stellar_temp > 0.0, "Raytrace stellar temperature must be positive"); + PARTHENON_REQUIRE(stellar_radius > 0.0, "Raytrace stellar radius must be positive"); stellar_radius *= constants.GetRsolarPhysical(); const Real sb = 0.25 * constants.GetCPhysical() * constants.GetARPhysical(); Real luminosity_cgs = 4 * M_PI * SQR(stellar_radius) * sb * SQR(SQR(stellar_temp)); @@ -57,12 +60,18 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("stellar_radius", stellar_radius * units.GetLengthPhysicalToCode()); const Real radius_factor = pin->GetOrAddReal("radiation/raytrace", "radius_factor", 6.); + PARTHENON_REQUIRE(radius_factor >= 0.0, + "Raytrace stellar radius factor must be non-negative"); params.Add("radius_factor", radius_factor); - params.Add("max_iterations", - pin->GetOrAddInteger("radiation/raytrace", "max_iterations", 1000)); + const int max_iterations = + pin->GetOrAddInteger("radiation/raytrace", "max_iterations", 1000); + PARTHENON_REQUIRE(max_iterations > 0, "Raytrace max_iterations must be positive"); + params.Add("max_iterations", max_iterations); - params.Add("efloor", pin->GetOrAddReal("radiation/raytrace", "efloor", 1e-10)); + const Real efloor = pin->GetOrAddReal("radiation/raytrace", "efloor", 1e-10); + PARTHENON_REQUIRE(efloor >= 0.0, "Raytrace photon energy floor must be non-negative"); + params.Add("efloor", efloor); // Note that these pull the real x1 values, not the ones from the artemis package params.Add("x1min", pin->GetReal("parthenon/mesh", "x1min")); @@ -242,6 +251,8 @@ TaskStatus CheckCompletion(MeshData *md) { auto &rt_pkg = pm->packages.Get("raytrace"); auto x1max = rt_pkg->Param("x1max"); + const Real x1tol = + 32.0 * std::numeric_limits::epsilon() * std::max(1.0, std::abs(x1max)); int num_unfinished = 0; parthenon::par_reduce( @@ -252,9 +263,7 @@ TaskStatus CheckCompletion(MeshData *md) { if (swarm_d.IsActive(n)) { const Real &xp = ppack_r(b, swarm_position::x(), n); const bool alive = ppack_r(b, rad::star::flux(), n) > 0.0; - const bool outside = - (xp >= x1max) || - (std::abs(xp - x1max) < 10 * std::numeric_limits::epsilon()); + const bool outside = (xp >= x1max) || (std::abs(xp - x1max) <= x1tol); num_unfinished += (alive && !outside); } }, @@ -298,8 +307,9 @@ TaskStatus EvalOpac(MeshData *md) { // Evaluated at T* //%%%%%%%%%%%%%%%% const Real temp = eos_d.TemperatureFromDensityInternalEnergy(rho, sie); + const Real alpha = opacity_d.AbsorptionCoefficient(rho, temp, 1.0); vmesh(b, rad::star::absorption(), k, j, i) = - opacity_d.AbsorptionCoefficient(rho, temp, 1.0); + (std::isfinite(alpha) && alpha > 0.0) ? alpha : 0.0; }); return TaskStatus::complete; diff --git a/src/radiation/raytrace/raytrace.hpp b/src/radiation/raytrace/raytrace.hpp index 3a548ebc..021c722d 100644 --- a/src/radiation/raytrace/raytrace.hpp +++ b/src/radiation/raytrace/raytrace.hpp @@ -13,6 +13,9 @@ #ifndef RADIATION_RAYTRACE_RAYTRACE_HPP_ #define RADIATION_RAYTRACE_RAYTRACE_HPP_ +// C++ headers +#include + // Artemis includes #include "artemis.hpp" #include "geometry/geometry.hpp" @@ -92,6 +95,8 @@ TaskStatus PushParticlesImpl(MeshData *md, const geometry::CoordParams &cp const int ngh = parthenon::Globals::nghost; const Real rmin = (LOGR) ? std::exp(x1min) : x1min; + const Real x1tol = 32.0 * std::numeric_limits::epsilon() * + std::max(1.0, std::max(std::abs(x1min), std::abs(x1max))); parthenon::par_for( DEFAULT_LOOP_PATTERN, "TransportPhotons", DevExecSpace(), 0, nparticles_per_pack, @@ -109,19 +114,25 @@ TaskStatus PushParticlesImpl(MeshData *md, const geometry::CoordParams &cp const auto &pco = vmesh.GetCoordinates(b); const auto inds = GetIndices(pco, {xp, yp, zp}); i = ib.s + inds[0] - ngh; + // A ray arrives exactly on a block face. Roundoff in floor() can otherwise + // select the last ghost cell and deposit outside the active mesh. + if (i == ib.s - 1 && std::abs(xp - pco.template Xf(ib.s)) <= x1tol) { + i = ib.s; + } if (multi_d) j = jb.s + inds[1] - ngh; if (three_d) k = kb.s + inds[2] - ngh; - while ((i <= ib.e) && (ee > 0.0)) { + while ((i >= ib.s) && (i <= ib.e) && (ee > 0.0)) { geometry::Coords coords(cpars, pco, k, j, i); // Deposit energy for this cell and decrement the photon energy const auto dx = coords.bnds.x1[1] - coords.bnds.x1[0]; - Real dtau = vmesh(b, rad::star::absorption(), k, j, i); + Real dtau = std::max(0.0, vmesh(b, rad::star::absorption(), k, j, i)); // Corrections for additional extinction inside the inner boundary - if (xp <= x1min + 1e-10) { - Real dtau_i = dtau * (rmin - zero_rad); + if (xp <= x1min + x1tol) { + const Real inner_path = std::max(0.0, rmin - zero_rad); + Real dtau_i = dtau * inner_path; const Real efac = (dtau_i > 100.) ? 0.0 : std::exp(-dtau_i); ee *= efac; if (ee < efloor) ee = 0.0; @@ -130,9 +141,11 @@ TaskStatus PushParticlesImpl(MeshData *md, const geometry::CoordParams &cp break; } } - dtau *= dx; + dtau *= std::max(0.0, dx); const Real efac = (dtau > 100.) ? 0.0 : std::exp(-dtau); - const Real reduc = (dtau <= 1e-4) ? dtau - 0.5 * SQR(dtau) : (1. - efac); + const Real reduc = (dtau <= 1e-4) + ? dtau - 0.5 * SQR(dtau) + dtau * SQR(dtau) / 6.0 + : (1. - efac); Real dE = ee * reduc; ee *= efac; @@ -148,7 +161,7 @@ TaskStatus PushParticlesImpl(MeshData *md, const geometry::CoordParams &cp } else { xp = coords.bnds.x1[1]; } - if (std::abs(xp - x1max) <= 1e-10) xp = x1max; + if (std::abs(xp - x1max) <= x1tol) xp = x1max; if ((ee == 0.0) || (xp >= x1max)) { swarm_d.MarkParticleForRemoval(n); break; From 3dce9cca788842b1e21fcf0a5c2e7a4bbb3188f2 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Mon, 13 Jul 2026 09:14:44 -0600 Subject: [PATCH 05/14] Use better tolerances --- src/radiation/moments/matter_coupling.hpp | 1049 ++++++++------------- 1 file changed, 412 insertions(+), 637 deletions(-) diff --git a/src/radiation/moments/matter_coupling.hpp b/src/radiation/moments/matter_coupling.hpp index bece8327..a7e8660e 100644 --- a/src/radiation/moments/matter_coupling.hpp +++ b/src/radiation/moments/matter_coupling.hpp @@ -40,10 +40,10 @@ KOKKOS_INLINE_FUNCTION Real CouplingFluxMagnitude(const std::array &F) return std::sqrt(SQR(F[0]) + SQR(F[1]) + SQR(F[2])); } -KOKKOS_INLINE_FUNCTION bool ComputeCouplingEnergyCoefficients( - const Real sigp, const Real sigs, const Real g, const Real g2, - const Real beta2, const Real bdbdp, const Real bdf, Real &ca, Real &cb, - Real &cd) { +KOKKOS_INLINE_FUNCTION bool +ComputeCouplingEnergyCoefficients(const Real sigp, const Real sigs, const Real g, + const Real g2, const Real beta2, const Real bdbdp, + const Real bdf, Real &ca, Real &cb, Real &cd) { // Algebraically equivalent to // ca = g * (sigp + sigs - g2 * sigs * (1 + bdbdp)); // cb = g * sigp; @@ -55,17 +55,16 @@ KOKKOS_INLINE_FUNCTION bool ComputeCouplingEnergyCoefficients( return std::isfinite(ca) && std::isfinite(cb) && std::isfinite(cd); } -KOKKOS_INLINE_FUNCTION bool ComputeCouplingEquilibriumEnergy( - const Real E0, const Real B, const Real ca, const Real cb, - const Real cd, Real &Eeq) { +KOKKOS_INLINE_FUNCTION bool ComputeCouplingEquilibriumEnergy(const Real E0, const Real B, + const Real ca, const Real cb, + const Real cd, Real &Eeq) { // The implicit radiation-energy equation is // (E - E0) + ca E - cb B + cd = 0. // Evaluate its solution using coefficient ratios. This remains accurate in // the optically thick limit, where ca and cb may be O(1e10) and directly // evaluating ca*E - cb*B loses all useful digits near LTE. const Real denom = 1.0 + ca; - if (!std::isfinite(denom) || std::abs(denom) <= Fuzz()) - return false; + if (!std::isfinite(denom) || std::abs(denom) <= Fuzz()) return false; const Real inv_denom = 1.0 / denom; Eeq = E0 * inv_denom + (cb * inv_denom) * B - cd * inv_denom; return std::isfinite(Eeq); @@ -74,16 +73,14 @@ KOKKOS_INLINE_FUNCTION bool ComputeCouplingEquilibriumEnergy( KOKKOS_INLINE_FUNCTION std::array ProjectCouplingFluxToEnergy(const std::array &Fin, const Real E) { const Real fmag = CouplingFluxMagnitude(Fin); - if (!std::isfinite(fmag) || !std::isfinite(E) || E <= 0.0) - return {0.0, 0.0, 0.0}; - if (fmag <= E || fmag <= Fuzz()) - return Fin; + if (!std::isfinite(fmag) || !std::isfinite(E) || E <= 0.0) return {0.0, 0.0, 0.0}; + if (fmag <= E || fmag <= Fuzz()) return Fin; const Real scale = E / fmag; return {scale * Fin[0], scale * Fin[1], scale * Fin[2]}; } -KOKKOS_INLINE_FUNCTION bool SolveCouplingDense3x3( - Real A[3][3], Real rhs[3], std::array &x) { +KOKKOS_INLINE_FUNCTION bool SolveCouplingDense3x3(Real A[3][3], Real rhs[3], + std::array &x) { // Small partial-pivoting Gaussian elimination used by the local momentum // Newton solve. The system is only 3x3, so forming and factorizing it is // cheaper and more robust than attempting an analytic inverse. @@ -97,8 +94,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingDense3x3( pivot_abs = candidate; } } - if (!std::isfinite(pivot_abs) || pivot_abs <= Fuzz()) - return false; + if (!std::isfinite(pivot_abs) || pivot_abs <= Fuzz()) return false; if (pivot != col) { for (int j = col; j < 3; ++j) { const Real tmp = A[col][j]; @@ -124,37 +120,31 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingDense3x3( Real value = rhs[row]; for (int j = row + 1; j < 3; ++j) value -= A[row][j] * x[j]; - if (!std::isfinite(A[row][row]) || - std::abs(A[row][row]) <= Fuzz()) + if (!std::isfinite(A[row][row]) || std::abs(A[row][row]) <= Fuzz()) return false; x[row] = value / A[row][row]; - if (!std::isfinite(x[row])) - return false; + if (!std::isfinite(x[row])) return false; } return true; } template KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualBeta( - const std::array &beta, - const std::array &beta0, - const std::array &Fr0, const Real E0, const Real efloor, - const Real B, const Real dEg, const Real Q, const Real dens, - const Real eref, const Real c, const Real chat, const Real sigp, - const Real sigs, const Real energy_floor_scale, - std::array &F, Real &E, Real &dEk, - std::array &residual, Real &residual_norm) { + const std::array &beta, const std::array &beta0, + const std::array &Fr0, const Real E0, const Real efloor, const Real B, + const Real dEg, const Real Q, const Real dens, const Real eref, const Real c, + const Real chat, const Real sigp, const Real sigs, const Real energy_floor_scale, + std::array &F, Real &E, Real &dEk, std::array &residual, + Real &residual_norm) { // Enforce gas+radiation pseudo-momentum conservation algebraically: // rho c (beta-beta0) + (eref/chat) (F-Fr0) = 0, // where F is normalized by c*eref. This relation is the source of the // instability in a direct Picard update when eref/(rho*c*chat) is large. const Real mu = eref / (dens * c * chat); - if (!std::isfinite(mu) || mu <= 0.0) - return false; + if (!std::isfinite(mu) || mu <= 0.0) return false; const Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); - if (!std::isfinite(beta2) || beta2 >= 1.0) - return false; + if (!std::isfinite(beta2) || beta2 >= 1.0) return false; std::array delta_F{0.0, 0.0, 0.0}; Real delta_F2 = 0.0; @@ -169,27 +159,21 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualBeta( // From beta-beta0 = -mu*(F-Fr0), evaluate the kinetic-energy // change without subtracting two nearly equal beta^2 values: // dEk = (c/chat)[-beta0.dF + 0.5*mu*|dF|^2]. - dEk = c / chat * - (-beta0_dot_delta_F + 0.5 * mu * delta_F2); + dEk = c / chat * (-beta0_dot_delta_F + 0.5 * mu * delta_F2); E = E0 + chat / c * (Q - dEg - dEk); - const Real floor_slop = - 128.0 * std::numeric_limits::epsilon() * - std::max(1.0, std::max(std::abs(E0), std::abs(E))); - if (!std::isfinite(E) || E < efloor - floor_slop || E <= 0.0) - return false; + const Real floor_slop = 128.0 * std::numeric_limits::epsilon() * + std::max(1.0, std::max(std::abs(E0), std::abs(E))); + if (!std::isfinite(E) || E < efloor - floor_slop || E <= 0.0) return false; E = std::max(E, efloor); const Real fmag = CouplingFluxMagnitude(F); const Real realizability_slop = - 256.0 * std::numeric_limits::epsilon() * - std::max(1.0, E); - if (!std::isfinite(fmag) || fmag > E + realizability_slop) - return false; + 256.0 * std::numeric_limits::epsilon() * std::max(1.0, E); + if (!std::isfinite(fmag) || fmag > E + realizability_slop) return false; const Real g2 = 1.0 / (1.0 - beta2); const Real g = std::sqrt(g2); - const auto fedd = EddingtonTensor( - {F[0] / E, F[1] / E, F[2] / E}); + const auto fedd = EddingtonTensor({F[0] / E, F[1] / E, F[2] / E}); const std::array bdp{ beta[0] * fedd[TensIdx::X11] + beta[1] * fedd[TensIdx::X12] + beta[2] * fedd[TensIdx::X13], @@ -197,58 +181,49 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualBeta( beta[2] * fedd[TensIdx::X23], beta[0] * fedd[TensIdx::X13] + beta[1] * fedd[TensIdx::X23] + beta[2] * fedd[TensIdx::X33]}; - const Real bdbdp = - beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; + const Real bdbdp = beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; const Real sigf = sigp + sigs; const Real a = g * sigf; const Real bcoef = 2.0 * g2 * g * sigs; - const Real d1 = - g * (sigp * B + g2 * sigs * (1.0 + bdbdp) * E); + const Real d1 = g * (sigp * B + g2 * sigs * (1.0 + bdbdp) * E); const Real d2 = g * sigf * E; - const std::array rhs{ - Fr0[0] + d1 * beta[0] + d2 * bdp[0], - Fr0[1] + d1 * beta[1] + d2 * bdp[1], - Fr0[2] + d1 * beta[2] + d2 * bdp[2]}; + const std::array rhs{Fr0[0] + d1 * beta[0] + d2 * bdp[0], + Fr0[1] + d1 * beta[1] + d2 * bdp[1], + Fr0[2] + d1 * beta[2] + d2 * bdp[2]}; const auto Fsolve = SolveRadFlux(1.0 + a, bcoef, beta, rhs); - if (!std::isfinite(Fsolve[0]) || !std::isfinite(Fsolve[1]) || - !std::isfinite(Fsolve[2])) + if (!std::isfinite(Fsolve[0]) || !std::isfinite(Fsolve[1]) || !std::isfinite(Fsolve[2])) return false; const auto Ftarget = ProjectCouplingFluxToEnergy(Fsolve, E); - const Real scale = std::max( - energy_floor_scale, - std::max(E, std::max(CouplingFluxMagnitude(Fr0), - std::max(fmag, CouplingFluxMagnitude(Ftarget))))); + const Real scale = + std::max(energy_floor_scale, + std::max(E, std::max(CouplingFluxMagnitude(Fr0), + std::max(fmag, CouplingFluxMagnitude(Ftarget))))); residual_norm = 0.0; for (int d = 0; d < 3; ++d) { residual[d] = F[d] - Ftarget[d]; - residual_norm = - std::max(residual_norm, std::abs(residual[d]) / scale); + residual_norm = std::max(residual_norm, std::abs(residual[d]) / scale); } return std::isfinite(residual_norm); } template KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( - const std::array &beta0, - const std::array &Fr0, const Real E0, const Real efloor, - const Real B, const Real dEg, const Real Q, const Real dens, - const Real eref, const Real c, const Real chat, const Real sigp, - const Real sigs, const Real energy_floor_scale, const Real tolerance, - std::array &beta, std::array &F, Real &E, - Real &dEk, Real &momentum_error, int &iterations) { + const std::array &beta0, const std::array &Fr0, const Real E0, + const Real efloor, const Real B, const Real dEg, const Real Q, const Real dens, + const Real eref, const Real c, const Real chat, const Real sigp, const Real sigs, + const Real energy_floor_scale, const Real tolerance, std::array &beta, + std::array &F, Real &E, Real &dEk, Real &momentum_error, int &iterations) { const Real mu = eref / (dens * c * chat); - if (!std::isfinite(mu) || mu <= 0.0 || !std::isfinite(sigp) || - !std::isfinite(sigs) || sigp < 0.0 || sigs < 0.0) + if (!std::isfinite(mu) || mu <= 0.0 || !std::isfinite(sigp) || !std::isfinite(sigs) || + sigp < 0.0 || sigs < 0.0) return false; - const Real solve_tol = - std::max(tolerance, 64.0 * std::numeric_limits::epsilon()); + const Real solve_tol = std::max(tolerance, 64.0 * std::numeric_limits::epsilon()); constexpr int max_iterations = 32; constexpr int max_line_search = 20; - const Real fd_factor = - std::pow(std::numeric_limits::epsilon(), 1.0 / 3.0); + const Real fd_factor = std::pow(std::numeric_limits::epsilon(), 1.0 / 3.0); std::array residual{0.0, 0.0, 0.0}; std::array best_beta = beta; @@ -260,9 +235,8 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( for (iterations = 1; iterations <= max_iterations; ++iterations) { Real current_error = std::numeric_limits::max(); if (!EvaluateCouplingMomentumEnergyResidualBeta( - beta, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, - chat, sigp, sigs, energy_floor_scale, F, E, dEk, residual, - current_error)) + beta, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, F, E, dEk, residual, current_error)) break; if (current_error < best_error) { @@ -277,14 +251,11 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( return true; } - Real J[3][3]{{0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0}}; + Real J[3][3]{{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; bool jacobian_valid = true; for (int col = 0; col < 3; ++col) { - const Real h = std::max( - 32.0 * std::numeric_limits::epsilon(), - fd_factor * std::max(1.0e-3, std::abs(beta[col]))); + const Real h = std::max(32.0 * std::numeric_limits::epsilon(), + fd_factor * std::max(1.0e-3, std::abs(beta[col]))); auto beta_plus = beta; auto beta_minus = beta; beta_plus[col] += h; @@ -298,16 +269,12 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( Real dEkminus = dEk; Real err_plus = 0.0; Real err_minus = 0.0; - const bool plus_valid = - EvaluateCouplingMomentumEnergyResidualBeta( - beta_plus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, - eref, c, chat, sigp, sigs, energy_floor_scale, Fplus, - Eplus, dEkplus, Rplus, err_plus); - const bool minus_valid = - EvaluateCouplingMomentumEnergyResidualBeta( - beta_minus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, - eref, c, chat, sigp, sigs, energy_floor_scale, Fminus, - Eminus, dEkminus, Rminus, err_minus); + const bool plus_valid = EvaluateCouplingMomentumEnergyResidualBeta( + beta_plus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, Fplus, Eplus, dEkplus, Rplus, err_plus); + const bool minus_valid = EvaluateCouplingMomentumEnergyResidualBeta( + beta_minus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, Fminus, Eminus, dEkminus, Rminus, err_minus); if (plus_valid && minus_valid) { for (int row = 0; row < 3; ++row) @@ -335,16 +302,13 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( // gain is approximately -(4/3) mu E; the denominator below removes that // stiffness while retaining the correct fixed point. if (!jacobian_valid) { - const Real relax = 1.0 / (1.0 + (4.0 / 3.0) * mu * - std::max(E, energy_floor_scale)); + const Real relax = 1.0 / (1.0 + (4.0 / 3.0) * mu * std::max(E, energy_floor_scale)); for (int d = 0; d < 3; ++d) step[d] = relax * mu * residual[d]; } - const Real step_mag = - std::sqrt(SQR(step[0]) + SQR(step[1]) + SQR(step[2])); - if (!std::isfinite(step_mag) || step_mag <= Fuzz()) - break; + const Real step_mag = std::sqrt(SQR(step[0]) + SQR(step[1]) + SQR(step[2])); + if (!std::isfinite(step_mag) || step_mag <= Fuzz()) break; if (step_mag > 0.25) { const Real scale = 0.25 / step_mag; for (int d = 0; d < 3; ++d) @@ -354,20 +318,16 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( bool accepted = false; Real alpha = 1.0; for (int ls = 0; ls < max_line_search; ++ls) { - std::array beta_trial{ - beta[0] + alpha * step[0], beta[1] + alpha * step[1], - beta[2] + alpha * step[2]}; + std::array beta_trial{beta[0] + alpha * step[0], beta[1] + alpha * step[1], + beta[2] + alpha * step[2]}; std::array Ftrial, Rtrial; Real Etrial = E; Real dEktrial = dEk; Real trial_error = std::numeric_limits::max(); - const bool trial_valid = - EvaluateCouplingMomentumEnergyResidualBeta( - beta_trial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, - eref, c, chat, sigp, sigs, energy_floor_scale, Ftrial, - Etrial, dEktrial, Rtrial, trial_error); - if (trial_valid && - trial_error < current_error * (1.0 - 1.0e-4 * alpha)) { + const bool trial_valid = EvaluateCouplingMomentumEnergyResidualBeta( + beta_trial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, Ftrial, Etrial, dEktrial, Rtrial, trial_error); + if (trial_valid && trial_error < current_error * (1.0 - 1.0e-4 * alpha)) { beta = beta_trial; F = Ftrial; E = Etrial; @@ -381,24 +341,21 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( if (!accepted) { // Retry with the stiffness-aware Picard direction even when the Newton // Jacobian existed but produced a poor globalization step. - const Real relax = 1.0 / (1.0 + (4.0 / 3.0) * mu * - std::max(E, energy_floor_scale)); + const Real relax = 1.0 / (1.0 + (4.0 / 3.0) * mu * std::max(E, energy_floor_scale)); for (int d = 0; d < 3; ++d) step[d] = relax * mu * residual[d]; alpha = 1.0; for (int ls = 0; ls < max_line_search; ++ls) { - std::array beta_trial{ - beta[0] + alpha * step[0], beta[1] + alpha * step[1], - beta[2] + alpha * step[2]}; + std::array beta_trial{beta[0] + alpha * step[0], + beta[1] + alpha * step[1], + beta[2] + alpha * step[2]}; std::array Ftrial, Rtrial; Real Etrial = E; Real dEktrial = dEk; Real trial_error = std::numeric_limits::max(); - const bool trial_valid = - EvaluateCouplingMomentumEnergyResidualBeta( - beta_trial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, - eref, c, chat, sigp, sigs, energy_floor_scale, Ftrial, - Etrial, dEktrial, Rtrial, trial_error); + const bool trial_valid = EvaluateCouplingMomentumEnergyResidualBeta( + beta_trial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, + sigs, energy_floor_scale, Ftrial, Etrial, dEktrial, Rtrial, trial_error); if (trial_valid && trial_error < current_error) { beta = beta_trial; F = Ftrial; @@ -410,8 +367,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( alpha *= 0.5; } } - if (!accepted) - break; + if (!accepted) break; } beta = best_beta; @@ -422,7 +378,6 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( return best_error <= solve_tol; } - // In the high-matter-inertia limit mu = eref/(rho*c*chat) << 1, the // physically resolvable momentum unknown is the radiation-flux change, not // beta. A finite flux correction corresponds to a beta correction of order @@ -430,17 +385,14 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( // beta then cannot represent the root even though F remains well resolved. template KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualFlux( - const std::array &F, - const std::array &beta0, - const std::array &Fr0, const Real E0, const Real efloor, - const Real B, const Real dEg, const Real Q, const Real dens, - const Real eref, const Real c, const Real chat, const Real sigp, - const Real sigs, const Real energy_floor_scale, - std::array &beta, Real &E, Real &dEk, - std::array &residual, Real &residual_norm) { + const std::array &F, const std::array &beta0, + const std::array &Fr0, const Real E0, const Real efloor, const Real B, + const Real dEg, const Real Q, const Real dens, const Real eref, const Real c, + const Real chat, const Real sigp, const Real sigs, const Real energy_floor_scale, + std::array &beta, Real &E, Real &dEk, std::array &residual, + Real &residual_norm) { const Real mu = eref / (dens * c * chat); - if (!std::isfinite(mu) || mu <= 0.0) - return false; + if (!std::isfinite(mu) || mu <= 0.0) return false; std::array delta_F{0.0, 0.0, 0.0}; Real delta_F2 = 0.0; @@ -453,32 +405,25 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualFlux( } const Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); - if (!std::isfinite(beta2) || beta2 >= 1.0) - return false; + if (!std::isfinite(beta2) || beta2 >= 1.0) return false; // This form remains accurate even when beta-beta0 is below the resolution // of a stored velocity component. - dEk = c / chat * - (-beta0_dot_delta_F + 0.5 * mu * delta_F2); + dEk = c / chat * (-beta0_dot_delta_F + 0.5 * mu * delta_F2); E = E0 + chat / c * (Q - dEg - dEk); - const Real floor_slop = - 128.0 * std::numeric_limits::epsilon() * - std::max(1.0, std::max(std::abs(E0), std::abs(E))); - if (!std::isfinite(E) || E < efloor - floor_slop || E <= 0.0) - return false; + const Real floor_slop = 128.0 * std::numeric_limits::epsilon() * + std::max(1.0, std::max(std::abs(E0), std::abs(E))); + if (!std::isfinite(E) || E < efloor - floor_slop || E <= 0.0) return false; E = std::max(E, efloor); const Real fmag = CouplingFluxMagnitude(F); const Real realizability_slop = - 256.0 * std::numeric_limits::epsilon() * - std::max(1.0, E); - if (!std::isfinite(fmag) || fmag > E + realizability_slop) - return false; + 256.0 * std::numeric_limits::epsilon() * std::max(1.0, E); + if (!std::isfinite(fmag) || fmag > E + realizability_slop) return false; const Real g2 = 1.0 / (1.0 - beta2); const Real g = std::sqrt(g2); - const auto fedd = EddingtonTensor( - {F[0] / E, F[1] / E, F[2] / E}); + const auto fedd = EddingtonTensor({F[0] / E, F[1] / E, F[2] / E}); const std::array bdp{ beta[0] * fedd[TensIdx::X11] + beta[1] * fedd[TensIdx::X12] + beta[2] * fedd[TensIdx::X13], @@ -486,58 +431,49 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualFlux( beta[2] * fedd[TensIdx::X23], beta[0] * fedd[TensIdx::X13] + beta[1] * fedd[TensIdx::X23] + beta[2] * fedd[TensIdx::X33]}; - const Real bdbdp = - beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; + const Real bdbdp = beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; const Real sigf = sigp + sigs; const Real a = g * sigf; const Real bcoef = 2.0 * g2 * g * sigs; - const Real d1 = - g * (sigp * B + g2 * sigs * (1.0 + bdbdp) * E); + const Real d1 = g * (sigp * B + g2 * sigs * (1.0 + bdbdp) * E); const Real d2 = g * sigf * E; - const std::array rhs{ - Fr0[0] + d1 * beta[0] + d2 * bdp[0], - Fr0[1] + d1 * beta[1] + d2 * bdp[1], - Fr0[2] + d1 * beta[2] + d2 * bdp[2]}; + const std::array rhs{Fr0[0] + d1 * beta[0] + d2 * bdp[0], + Fr0[1] + d1 * beta[1] + d2 * bdp[1], + Fr0[2] + d1 * beta[2] + d2 * bdp[2]}; const auto Fsolve = SolveRadFlux(1.0 + a, bcoef, beta, rhs); - if (!std::isfinite(Fsolve[0]) || !std::isfinite(Fsolve[1]) || - !std::isfinite(Fsolve[2])) + if (!std::isfinite(Fsolve[0]) || !std::isfinite(Fsolve[1]) || !std::isfinite(Fsolve[2])) return false; const auto Ftarget = ProjectCouplingFluxToEnergy(Fsolve, E); - const Real scale = std::max( - energy_floor_scale, - std::max(E, std::max(CouplingFluxMagnitude(Fr0), - std::max(fmag, CouplingFluxMagnitude(Ftarget))))); + const Real scale = + std::max(energy_floor_scale, + std::max(E, std::max(CouplingFluxMagnitude(Fr0), + std::max(fmag, CouplingFluxMagnitude(Ftarget))))); residual_norm = 0.0; for (int d = 0; d < 3; ++d) { residual[d] = F[d] - Ftarget[d]; - residual_norm = - std::max(residual_norm, std::abs(residual[d]) / scale); + residual_norm = std::max(residual_norm, std::abs(residual[d]) / scale); } return std::isfinite(residual_norm); } template KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( - const std::array &beta0, - const std::array &Fr0, const Real E0, const Real efloor, - const Real B, const Real dEg, const Real Q, const Real dens, - const Real eref, const Real c, const Real chat, const Real sigp, - const Real sigs, const Real energy_floor_scale, const Real tolerance, - std::array &beta, std::array &F, Real &E, - Real &dEk, Real &momentum_error, int &iterations) { + const std::array &beta0, const std::array &Fr0, const Real E0, + const Real efloor, const Real B, const Real dEg, const Real Q, const Real dens, + const Real eref, const Real c, const Real chat, const Real sigp, const Real sigs, + const Real energy_floor_scale, const Real tolerance, std::array &beta, + std::array &F, Real &E, Real &dEk, Real &momentum_error, int &iterations) { const Real mu = eref / (dens * c * chat); - if (!std::isfinite(mu) || mu <= 0.0 || !std::isfinite(sigp) || - !std::isfinite(sigs) || sigp < 0.0 || sigs < 0.0) + if (!std::isfinite(mu) || mu <= 0.0 || !std::isfinite(sigp) || !std::isfinite(sigs) || + sigp < 0.0 || sigs < 0.0) return false; - const Real solve_tol = - std::max(tolerance, 64.0 * std::numeric_limits::epsilon()); + const Real solve_tol = std::max(tolerance, 64.0 * std::numeric_limits::epsilon()); constexpr int max_iterations = 32; constexpr int max_line_search = 20; - const Real fd_factor = - std::pow(std::numeric_limits::epsilon(), 1.0 / 3.0); + const Real fd_factor = std::pow(std::numeric_limits::epsilon(), 1.0 / 3.0); std::array residual{0.0, 0.0, 0.0}; std::array best_beta = beta; @@ -549,9 +485,8 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( for (iterations = 1; iterations <= max_iterations; ++iterations) { Real current_error = std::numeric_limits::max(); if (!EvaluateCouplingMomentumEnergyResidualFlux( - F, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, - chat, sigp, sigs, energy_floor_scale, beta, E, dEk, - residual, current_error)) + F, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, beta, E, dEk, residual, current_error)) break; if (current_error < best_error) { @@ -566,19 +501,18 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( return true; } - Real J[3][3]{{0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0}}; + Real J[3][3]{{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; bool jacobian_valid = true; for (int col = 0; col < 3; ++col) { - const Real variable_scale = std::max( - energy_floor_scale, - std::max(std::abs(F[col]), - std::max(std::abs(Fr0[col]), 1.0e-6 * E))); - const Real h = std::max( - 64.0 * std::numeric_limits::epsilon() * - std::max(1.0, std::abs(F[col])), - fd_factor * variable_scale); + const Real variable_scale = + std::max(energy_floor_scale, + std::max(std::abs(F[col]), std::max(std::abs(Fr0[col]), 1.0e-6 * E))); + // F is normalized by the local energy reference and can be much smaller + // than unity. An absolute O(epsilon) perturbation would then be larger + // than the correction implied by the requested relative residual. + const Real h = std::max(64.0 * std::numeric_limits::epsilon() * + variable_scale, + fd_factor * variable_scale); auto Fplus = F; auto Fminus = F; Fplus[col] += h; @@ -592,16 +526,12 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( Real dEkminus = dEk; Real err_plus = 0.0; Real err_minus = 0.0; - const bool plus_valid = - EvaluateCouplingMomentumEnergyResidualFlux( - Fplus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, - eref, c, chat, sigp, sigs, energy_floor_scale, - beta_plus, Eplus, dEkplus, Rplus, err_plus); - const bool minus_valid = - EvaluateCouplingMomentumEnergyResidualFlux( - Fminus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, - eref, c, chat, sigp, sigs, energy_floor_scale, - beta_minus, Eminus, dEkminus, Rminus, err_minus); + const bool plus_valid = EvaluateCouplingMomentumEnergyResidualFlux( + Fplus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, beta_plus, Eplus, dEkplus, Rplus, err_plus); + const bool minus_valid = EvaluateCouplingMomentumEnergyResidualFlux( + Fminus, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, beta_minus, Eminus, dEkminus, Rminus, err_minus); if (plus_valid && minus_valid) { for (int row = 0; row < 3; ++row) @@ -628,15 +558,16 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( step[d] = -residual[d]; } - const Real step_mag = - std::sqrt(SQR(step[0]) + SQR(step[1]) + SQR(step[2])); + const Real step_mag = std::sqrt(SQR(step[0]) + SQR(step[1]) + SQR(step[2])); const Real state_scale = std::max( energy_floor_scale, - std::max(E, std::max(CouplingFluxMagnitude(F), - CouplingFluxMagnitude(Fr0)))); - if (!std::isfinite(step_mag) || - step_mag <= 64.0 * std::numeric_limits::epsilon() * - std::max(1.0, state_scale)) + std::max(E, std::max(CouplingFluxMagnitude(F), CouplingFluxMagnitude(Fr0)))); + // state_scale includes the radiation floor. Do not impose a unit-scale + // absolute cutoff here: for E << 1 it can reject a resolvable correction + // before the normalized momentum residual reaches its tolerance. + if (!std::isfinite(step_mag) || step_mag <= 64.0 * + std::numeric_limits::epsilon() * + state_scale) break; if (step_mag > 0.25 * state_scale) { const Real scale = 0.25 * state_scale / step_mag; @@ -647,20 +578,16 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( bool accepted = false; Real alpha = 1.0; for (int ls = 0; ls < max_line_search; ++ls) { - std::array Ftrial{ - F[0] + alpha * step[0], F[1] + alpha * step[1], - F[2] + alpha * step[2]}; + std::array Ftrial{F[0] + alpha * step[0], F[1] + alpha * step[1], + F[2] + alpha * step[2]}; std::array beta_trial, Rtrial; Real Etrial = E; Real dEktrial = dEk; Real trial_error = std::numeric_limits::max(); - const bool trial_valid = - EvaluateCouplingMomentumEnergyResidualFlux( - Ftrial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, - eref, c, chat, sigp, sigs, energy_floor_scale, - beta_trial, Etrial, dEktrial, Rtrial, trial_error); - if (trial_valid && - trial_error < current_error * (1.0 - 1.0e-4 * alpha)) { + const bool trial_valid = EvaluateCouplingMomentumEnergyResidualFlux( + Ftrial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, beta_trial, Etrial, dEktrial, Rtrial, trial_error); + if (trial_valid && trial_error < current_error * (1.0 - 1.0e-4 * alpha)) { F = Ftrial; beta = beta_trial; E = Etrial; @@ -677,18 +604,15 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( step[d] = -residual[d]; alpha = 1.0; for (int ls = 0; ls < max_line_search; ++ls) { - std::array Ftrial{ - F[0] + alpha * step[0], F[1] + alpha * step[1], - F[2] + alpha * step[2]}; + std::array Ftrial{F[0] + alpha * step[0], F[1] + alpha * step[1], + F[2] + alpha * step[2]}; std::array beta_trial, Rtrial; Real Etrial = E; Real dEktrial = dEk; Real trial_error = std::numeric_limits::max(); - const bool trial_valid = - EvaluateCouplingMomentumEnergyResidualFlux( - Ftrial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, - eref, c, chat, sigp, sigs, energy_floor_scale, - beta_trial, Etrial, dEktrial, Rtrial, trial_error); + const bool trial_valid = EvaluateCouplingMomentumEnergyResidualFlux( + Ftrial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, beta_trial, Etrial, dEktrial, Rtrial, trial_error); if (trial_valid && trial_error < current_error) { F = Ftrial; beta = beta_trial; @@ -700,8 +624,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( alpha *= 0.5; } } - if (!accepted) - break; + if (!accepted) break; } beta = best_beta; @@ -714,13 +637,11 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( template KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergy( - const std::array &beta0, - const std::array &Fr0, const Real E0, const Real efloor, - const Real B, const Real dEg, const Real Q, const Real dens, - const Real eref, const Real c, const Real chat, const Real sigp, - const Real sigs, const Real energy_floor_scale, const Real tolerance, - std::array &beta, std::array &F, Real &E, - Real &dEk, Real &momentum_error, int &iterations) { + const std::array &beta0, const std::array &Fr0, const Real E0, + const Real efloor, const Real B, const Real dEg, const Real Q, const Real dens, + const Real eref, const Real c, const Real chat, const Real sigp, const Real sigs, + const Real energy_floor_scale, const Real tolerance, std::array &beta, + std::array &F, Real &E, Real &dEk, Real &momentum_error, int &iterations) { const Real mu = eref / (dens * c * chat); // Choose the nonlinear variable from the relative pseudo-inertia. For // mu <= 1 the gas velocity change can be much less resolvable than the @@ -728,41 +649,34 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergy( // magnifying flux perturbations into large velocity perturbations. if (mu <= 1.0) { return SolveCouplingMomentumEnergyFlux( - beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, - sigp, sigs, energy_floor_scale, tolerance, beta, F, E, dEk, - momentum_error, iterations); + beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, tolerance, beta, F, E, dEk, momentum_error, iterations); } return SolveCouplingMomentumEnergyBeta( - beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, - sigp, sigs, energy_floor_scale, tolerance, beta, F, E, dEk, - momentum_error, iterations); + beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, + energy_floor_scale, tolerance, beta, F, E, dEk, momentum_error, iterations); } template KOKKOS_INLINE_FUNCTION bool EvaluateCouplingInnerScalarResidual( - const Real eg, const Real eg0, const Real dEk, const Real Q, - const Real E0, const Real efloor, const Real Bfloor, const Real dens, - const Real eref, const Real arad, const Real tfloor, const Real dt, - const Real chat, const Real c, const Real g, const Real g2, - const Real beta2, const Real bdbdp, const Real bdf, const EOSType &eos, - const OpacityType &opacity, const ScatteringType &scattering, Real &E, - Real &B, Real &residual) { + const Real eg, const Real eg0, const Real dEk, const Real Q, const Real E0, + const Real efloor, const Real Bfloor, const Real dens, const Real eref, + const Real arad, const Real tfloor, const Real dt, const Real chat, const Real c, + const Real g, const Real g2, const Real beta2, const Real bdbdp, const Real bdf, + const EOSType &eos, const OpacityType &opacity, const ScatteringType &scattering, + Real &E, Real &B, Real &residual) { // Exact reduced-speed-of-light total-energy conservation at fixed velocity: // dEg + dEk + (c/chat) dEr = Q. E = E0 + chat / c * (Q - dEk - (eg - eg0)); - const Real floor_slop = - 64.0 * std::numeric_limits::epsilon() * - std::max(1.0, std::max(std::abs(E0), std::abs(E))); - if (!std::isfinite(E) || E < efloor - floor_slop) - return false; + const Real floor_slop = 64.0 * std::numeric_limits::epsilon() * + std::max(1.0, std::max(std::abs(E0), std::abs(E))); + if (!std::isfinite(E) || E < efloor - floor_slop) return false; E = std::max(E, efloor); const Real eint = eg * eref / dens; - const Real T = std::max( - tfloor, eos.TemperatureFromDensityInternalEnergy(dens, eint)); + const Real T = std::max(tfloor, eos.TemperatureFromDensityInternalEnergy(dens, eint)); B = std::max(Bfloor, arad * SQR(SQR(T)) / eref); - const Real sigp = - chat * dt * opacity.PlanckMeanAbsorptionCoefficient(dens, T); + const Real sigp = chat * dt * opacity.PlanckMeanAbsorptionCoefficient(dens, T); const Real sigs = chat * dt * scattering.RosselandMeanTotalScatteringCoefficient(dens, T); if (!std::isfinite(T) || !std::isfinite(B) || !std::isfinite(sigp) || @@ -773,8 +687,8 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingInnerScalarResidual( Real cb = 0.0; Real cd = 0.0; Real Eeq = E; - if (!ComputeCouplingEnergyCoefficients(sigp, sigs, g, g2, beta2, - bdbdp, bdf, ca, cb, cd) || + if (!ComputeCouplingEnergyCoefficients(sigp, sigs, g, g2, beta2, bdbdp, bdf, ca, cb, + cd) || !ComputeCouplingEquilibriumEnergy(E0, B, ca, cb, cd, Eeq)) return false; @@ -788,67 +702,57 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingInnerScalarResidual( return std::isfinite(residual); } -KOKKOS_INLINE_FUNCTION Real CouplingInnerScalarError( - const Real eg, const Real eg0, const Real dEk, const Real Q, - const Real E, const Real E0, const Real chat, const Real c, - const Real energy_floor_scale, const Real residual) { +KOKKOS_INLINE_FUNCTION Real CouplingInnerScalarError(const Real eg, const Real eg0, + const Real dEk, const Real Q, + const Real E, const Real E0, + const Real chat, const Real c, + const Real energy_floor_scale, + const Real residual) { const Real scale = std::max( energy_floor_scale, - std::max( - std::max(std::abs(eg0), std::abs(eg)), - std::max(std::abs(Q), - std::max(c / chat * (std::abs(E0) + std::abs(E)), - std::abs(dEk))))); + std::max(std::max(std::abs(eg0), std::abs(eg)), + std::max(std::abs(Q), std::max(c / chat * (std::abs(E0) + std::abs(E)), + std::abs(dEk))))); return c / chat * std::abs(residual) / scale; } -KOKKOS_INLINE_FUNCTION bool CouplingResidualChangesSign(const Real a, - const Real b) { +KOKKOS_INLINE_FUNCTION bool CouplingResidualChangesSign(const Real a, const Real b) { return (a <= 0.0 && b >= 0.0) || (a >= 0.0 && b <= 0.0); } template KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( - const Real eg0, const Real dEk, const Real Q, const Real E0, - const Real efloor, const Real Bfloor, const Real dens, const Real eref, - const Real arad, const Real tfloor, const Real dt, const Real chat, - const Real c, const Real g, const Real g2, const Real beta2, - const Real bdbdp, const Real bdf, const Real energy_floor_scale, - const Real tolerance, - const EOSType &eos, const OpacityType &opacity, - const ScatteringType &scattering, Real &E, Real &B, Real &inner_err, - int &iterations) { + const Real eg0, const Real dEk, const Real Q, const Real E0, const Real efloor, + const Real Bfloor, const Real dens, const Real eref, const Real arad, + const Real tfloor, const Real dt, const Real chat, const Real c, const Real g, + const Real g2, const Real beta2, const Real bdbdp, const Real bdf, + const Real energy_floor_scale, const Real tolerance, const EOSType &eos, + const OpacityType &opacity, const ScatteringType &scattering, Real &E, Real &B, + Real &inner_err, int &iterations) { const Real eg_floor = dens * eos.InternalEnergyFromDensityTemperature(dens, tfloor) / eref; - const Real eg_ceiling = - eg0 + Q - dEk + c / chat * (E0 - efloor); - if (!std::isfinite(eg_floor) || !std::isfinite(eg_ceiling) || - eg_ceiling < eg_floor) + const Real eg_ceiling = eg0 + Q - dEk + c / chat * (E0 - efloor); + if (!std::isfinite(eg_floor) || !std::isfinite(eg_ceiling) || eg_ceiling < eg_floor) return false; - const Real Tguess = std::max( - tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); - Real eg_guess = - dens * eos.InternalEnergyFromDensityTemperature(dens, Tguess) / eref; + const Real Tguess = std::max(tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); + Real eg_guess = dens * eos.InternalEnergyFromDensityTemperature(dens, Tguess) / eref; eg_guess = std::max(eg_floor, std::min(eg_ceiling, eg_guess)); Real Eguess = E; Real Bguess = B; Real Rguess = 0.0; bool guess_valid = EvaluateCouplingInnerScalarResidual( - eg_guess, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, - tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, - scattering, - Eguess, Bguess, Rguess); + eg_guess, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, tfloor, dt, chat, c, g, + g2, beta2, bdbdp, bdf, eos, opacity, scattering, Eguess, Bguess, Rguess); Real best_eg = eg_guess; Real best_E = Eguess; Real best_B = Bguess; - Real best_error = - guess_valid - ? CouplingInnerScalarError(eg_guess, eg0, dEk, Q, Eguess, E0, - chat, c, energy_floor_scale, Rguess) - : std::numeric_limits::max(); + Real best_error = guess_valid + ? CouplingInnerScalarError(eg_guess, eg0, dEk, Q, Eguess, E0, + chat, c, energy_floor_scale, Rguess) + : std::numeric_limits::max(); if (guess_valid && best_error <= tolerance) { E = Eguess; B = Bguess; @@ -861,13 +765,11 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( Real Blo = B; Real Rlo = 0.0; const bool lo_valid = EvaluateCouplingInnerScalarResidual( - eg_floor, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, - tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, - scattering, - Elo, Blo, Rlo); + eg_floor, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, tfloor, dt, chat, c, g, + g2, beta2, bdbdp, bdf, eos, opacity, scattering, Elo, Blo, Rlo); if (lo_valid) { const Real err = CouplingInnerScalarError(eg_floor, eg0, dEk, Q, Elo, E0, chat, c, - energy_floor_scale, Rlo); + energy_floor_scale, Rlo); if (err < best_error) { best_eg = eg_floor; best_E = Elo; @@ -880,13 +782,11 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( Real Bhi = B; Real Rhi = 0.0; const bool hi_valid = EvaluateCouplingInnerScalarResidual( - eg_ceiling, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, - tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, - scattering, - Ehi, Bhi, Rhi); + eg_ceiling, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, tfloor, dt, chat, c, + g, g2, beta2, bdbdp, bdf, eos, opacity, scattering, Ehi, Bhi, Rhi); if (hi_valid) { const Real err = CouplingInnerScalarError(eg_ceiling, eg0, dEk, Q, Ehi, E0, chat, c, - energy_floor_scale, Rhi); + energy_floor_scale, Rhi); if (err < best_error) { best_eg = eg_ceiling; best_E = Ehi; @@ -922,13 +822,11 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( Real Bleft = B; Real Rleft = 0.0; const bool left_valid = EvaluateCouplingInnerScalarResidual( - next_left, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, - tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, - scattering, Eleft, Bleft, Rleft); + next_left, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, tfloor, dt, chat, + c, g, g2, beta2, bdbdp, bdf, eos, opacity, scattering, Eleft, Bleft, Rleft); if (left_valid) { - const Real err = CouplingInnerScalarError( - next_left, eg0, dEk, Q, Eleft, E0, chat, c, - energy_floor_scale, Rleft); + const Real err = CouplingInnerScalarError(next_left, eg0, dEk, Q, Eleft, E0, chat, + c, energy_floor_scale, Rleft); if (err < best_error) { best_eg = next_left; best_E = Eleft; @@ -952,13 +850,11 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( Real Bright = B; Real Rright = 0.0; const bool right_valid = EvaluateCouplingInnerScalarResidual( - next_right, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, - tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, - scattering, Eright, Bright, Rright); + next_right, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, tfloor, dt, chat, + c, g, g2, beta2, bdbdp, bdf, eos, opacity, scattering, Eright, Bright, Rright); if (right_valid) { - const Real err = CouplingInnerScalarError( - next_right, eg0, dEk, Q, Eright, E0, chat, c, - energy_floor_scale, Rright); + const Real err = CouplingInnerScalarError(next_right, eg0, dEk, Q, Eright, E0, + chat, c, energy_floor_scale, Rright); if (err < best_error) { best_eg = next_right; best_E = Eright; @@ -999,8 +895,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( const Real eg_secant = (bracket_lo * bracket_Rhi - bracket_hi * bracket_Rlo) / denom; const Real width = bracket_hi - bracket_lo; - if (eg_secant > bracket_lo + 0.1 * width && - eg_secant < bracket_hi - 0.1 * width) + if (eg_secant > bracket_lo + 0.1 * width && eg_secant < bracket_hi - 0.1 * width) eg_trial = eg_secant; } @@ -1008,23 +903,19 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( Real Btrial = B; Real Rtrial = 0.0; const bool trial_valid = EvaluateCouplingInnerScalarResidual( - eg_trial, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, - tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, eos, opacity, - scattering, Etrial, Btrial, Rtrial); - if (!trial_valid) - break; + eg_trial, eg0, dEk, Q, E0, efloor, Bfloor, dens, eref, arad, tfloor, dt, chat, + c, g, g2, beta2, bdbdp, bdf, eos, opacity, scattering, Etrial, Btrial, Rtrial); + if (!trial_valid) break; - const Real err = CouplingInnerScalarError( - eg_trial, eg0, dEk, Q, Etrial, E0, chat, c, energy_floor_scale, - Rtrial); + const Real err = CouplingInnerScalarError(eg_trial, eg0, dEk, Q, Etrial, E0, chat, + c, energy_floor_scale, Rtrial); if (err < best_error) { best_eg = eg_trial; best_E = Etrial; best_B = Btrial; best_error = err; } - if (err <= tolerance) - break; + if (err <= tolerance) break; if (CouplingResidualChangesSign(bracket_Rlo, Rtrial)) { bracket_hi = eg_trial; @@ -1034,9 +925,8 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( bracket_Rlo = Rtrial; } - const Real width_scale = std::max( - energy_floor_scale, - std::abs(bracket_lo) + std::abs(bracket_hi)); + const Real width_scale = + std::max(energy_floor_scale, std::abs(bracket_lo) + std::abs(bracket_hi)); if ((bracket_hi - bracket_lo) / width_scale <= 8.0 * std::numeric_limits::epsilon()) break; @@ -1054,23 +944,18 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( template KOKKOS_INLINE_FUNCTION bool EvaluateCouplingEnergyOnlyResidual( - const Real eg, const Real eg0, const Real Q, const Real E0, - const Real Emin, const Real dens, const Real eref, const Real arad, - const Real tfloor, const Real dt, const Real chat, const Real c, - const EOSType &eos, const OpacityType &opacity, Real &E, Real &B, - Real &residual) { + const Real eg, const Real eg0, const Real Q, const Real E0, const Real Emin, + const Real dens, const Real eref, const Real arad, const Real tfloor, const Real dt, + const Real chat, const Real c, const EOSType &eos, const OpacityType &opacity, + Real &E, Real &B, Real &residual) { E = E0 + chat / c * (Q - (eg - eg0)); - if (!std::isfinite(E) || E < Emin) - return false; + if (!std::isfinite(E) || E < Emin) return false; const Real eint = eg * eref / dens; - const Real T = std::max( - tfloor, eos.TemperatureFromDensityInternalEnergy(dens, eint)); + const Real T = std::max(tfloor, eos.TemperatureFromDensityInternalEnergy(dens, eint)); B = arad * SQR(SQR(T)) / eref; - const Real sigp = - chat * dt * opacity.PlanckMeanAbsorptionCoefficient(dens, T); - if (!std::isfinite(T) || !std::isfinite(B) || !std::isfinite(sigp) || - sigp < 0.0) + const Real sigp = chat * dt * opacity.PlanckMeanAbsorptionCoefficient(dens, T); + if (!std::isfinite(T) || !std::isfinite(B) || !std::isfinite(sigp) || sigp < 0.0) return false; // This fallback intentionally retains only thermal absorption/emission. @@ -1110,16 +995,11 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const auto rad_efloor = moments_pkg->template Param("efloor"); const auto tfloor = moments_pkg->template Param("tfloor"); const auto Bfloor_phys = arad * SQR(SQR(tfloor)); - const auto outer_max = - moments_pkg->template Param("outer_iteration_max"); - const auto inner_max = - moments_pkg->template Param("inner_iteration_max"); - const auto outer_tol = - moments_pkg->template Param("outer_iteration_tol"); - const auto inner_tol = - moments_pkg->template Param("inner_iteration_tol"); - const Real nonlinear_roundoff_tol = - 64.0 * std::numeric_limits::epsilon(); + const auto outer_max = moments_pkg->template Param("outer_iteration_max"); + const auto inner_max = moments_pkg->template Param("inner_iteration_max"); + const auto outer_tol = moments_pkg->template Param("outer_iteration_tol"); + const auto inner_tol = moments_pkg->template Param("inner_iteration_tol"); + const Real nonlinear_roundoff_tol = 64.0 * std::numeric_limits::epsilon(); const auto fatal_if_unconverged = moments_pkg->template Param("fatal_if_unconverged"); @@ -1127,8 +1007,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { Real om0 = 0.0; Real qshear = 0.0; Real gm_bg = 0.0; - if (pm->packages.Get("artemis")->template Param( - "do_orbital_advection") || + if (pm->packages.Get("artemis")->template Param("do_orbital_advection") || pm->packages.Get("artemis")->template Param("do_rotating_frame")) { auto &rframe_pkg = pm->packages.Get("rotating_frame"); qshear = rframe_pkg->template Param("qshear"); @@ -1136,32 +1015,28 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { gm_bg = rframe_pkg->template Param("gm"); } const auto &cpars = - pm->packages.Get("artemis")->template Param( - "coord_params"); + pm->packages.Get("artemis")->template Param("coord_params"); const bool do_raytrace = pm->packages.Get("artemis")->template Param("do_raytrace"); // Packing and indexing - static auto desc = - parthenon::MakePackDescriptor( - resolved_pkgs.get()); + static auto desc = parthenon::MakePackDescriptor< + rad::cons::energy, rad::cons::flux, gas::cons::density, gas::cons::momentum, + gas::cons::internal_energy, gas::cons::total_energy, gas::src::energy>( + resolved_pkgs.get()); const auto v0 = desc.GetPack(u0); - static auto desc_g = - MakePackDescriptor(resolved_pkgs.get()); + static auto desc_g = MakePackDescriptor(resolved_pkgs.get()); auto vg = desc_g.GetPack(u0); const auto ib = u0->GetBoundsI(IndexDomain::interior); const auto jb = u0->GetBoundsJ(IndexDomain::interior); const auto kb = u0->GetBoundsK(IndexDomain::interior); parthenon::par_for( - DEFAULT_LOOP_PATTERN, "MatterCoupling", DevExecSpace(), 0, - u0->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + DEFAULT_LOOP_PATTERN, "MatterCoupling", DevExecSpace(), 0, u0->NumBlocks() - 1, + kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { geometry::Coords coords(cpars, v0.GetCoordinates(b), k, j, i); const auto &hx = coords.GetScaleFactors(vg, b, k, j, i); @@ -1169,11 +1044,9 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // U^(0) values const Real dens_raw = v0(b, gas::cons::density(), k, j, i); - const Real dens = - (std::isfinite(dens_raw) && dens_raw > dflr) ? dens_raw : dflr; + const Real dens = (std::isfinite(dens_raw) && dens_raw > dflr) ? dens_raw : dflr; Real Q = 0.0; - if (do_raytrace) - Q = FiniteOrZero(dt * v0(b, gas::src::energy(), k, j, i)); + if (do_raytrace) Q = FiniteOrZero(dt * v0(b, gas::src::energy(), k, j, i)); // In the numerical atmosphere, retain thermal emission/absorption but // suppress the momentum update. This avoids accelerating floor-density @@ -1181,8 +1054,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // the moment field through the conservative energy-only fallback. constexpr Real atmosphere_floor_factor = 100.0; const bool numerical_atmosphere = - !std::isfinite(dens_raw) || - dens_raw <= atmosphere_floor_factor * dflr; + !std::isfinite(dens_raw) || dens_raw <= atmosphere_floor_factor * dflr; // Note(AMD): There is some floating point difference between the // internal energy used to compute the temperature and the internal @@ -1197,8 +1069,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // increment. const Real eint0 = std::max(0.0, v0(b, gas::cons::internal_energy(), k, j, i) / dens); - Real T = std::max( - tfloor, eos_d.TemperatureFromDensityInternalEnergy(dens, eint0)); + Real T = + std::max(tfloor, eos_d.TemperatureFromDensityInternalEnergy(dens, eint0)); Real eg0 = dens * eos_d.InternalEnergyFromDensityTemperature(dens, T); Real B = std::max(Bfloor_phys, arad * SQR(SQR(T))); @@ -1214,18 +1086,15 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // have E_r << a_r T_floor^4. Promoting E_r to Bfloor_phys here creates // an artificial LTE radiation bath and can make the constrained // matter-coupling equations inconsistent at the first timestep. - const Real Er_state = - FiniteOrZero(v0(b, rad::cons::energy(), k, j, i)); + const Real Er_state = FiniteOrZero(v0(b, rad::cons::energy(), k, j, i)); Real E0 = std::max(rad_efloor, Er_state); // Use the arithmetic mean as the reference scale. Unlike the // geometric mean, this keeps both normalized energies bounded when the // gas and radiation temperatures are initially very different. const Real eref_max = std::max(E0, B); const Real eref_min = std::min(E0, B); - Real eref = - 0.5 * eref_max * (1.0 + eref_min / std::max(eref_max, Fuzz())); - if (!std::isfinite(eref) || eref <= 0.0) - eref = eref_max; + Real eref = 0.5 * eref_max * (1.0 + eref_min / std::max(eref_max, Fuzz())); + if (!std::isfinite(eref) || eref <= 0.0) eref = eref_max; if (!std::isfinite(eref) || eref <= 0.0) { v0(b, gas::cons::internal_energy(), k, j, i) += Q; v0(b, gas::cons::total_energy(), k, j, i) += Q; @@ -1247,14 +1116,10 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { B /= eref; const auto fred0 = NormalizeFlux( - FiniteOrZero(v0(b, rad::cons::flux(0), k, j, i) / - (hx[0] * fref * E0)), - FiniteOrZero(v0(b, rad::cons::flux(1), k, j, i) / - (hx[1] * fref * E0)), - FiniteOrZero(v0(b, rad::cons::flux(2), k, j, i) / - (hx[2] * fref * E0))); - const std::array Fr0{E0 * fred0[0], E0 * fred0[1], - E0 * fred0[2]}; + FiniteOrZero(v0(b, rad::cons::flux(0), k, j, i) / (hx[0] * fref * E0)), + FiniteOrZero(v0(b, rad::cons::flux(1), k, j, i) / (hx[1] * fref * E0)), + FiniteOrZero(v0(b, rad::cons::flux(2), k, j, i) / (hx[2] * fref * E0))); + const std::array Fr0{E0 * fred0[0], E0 * fred0[1], E0 * fred0[2]}; std::array v{p0[0] / dens, p0[1] / dens, p0[2] / dens}; const std::array beta0{v[0] / c, v[1] / c, v[2] / c}; @@ -1268,8 +1133,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { v0(b, rad::cons::flux(2), k, j, i) = Fr0[2] * hx[2] * fref; return; } - const Real ke0 = - 0.5 * dens * (SQR(v[0]) + SQR(v[1]) + SQR(v[2])) / eref; + const Real ke0 = 0.5 * dens * (SQR(v[0]) + SQR(v[1]) + SQR(v[2])) / eref; Real E = E0; auto F = Fr0; @@ -1284,12 +1148,10 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { Real dEk = 0.0; Real dEg = 0.0; Real dEr = 0.0; - const Real energy_floor_scale = - std::max(efloor + Bfloor, Fuzz()); + const Real energy_floor_scale = std::max(efloor + Bfloor, Fuzz()); Real escale = std::max( energy_floor_scale, - std::max(std::abs(eg0), - std::max(std::abs(Q), c / chat * std::abs(E0)))); + std::max(std::abs(eg0), std::max(std::abs(Q), c / chat * std::abs(E0)))); bool solve_valid = !numerical_atmosphere && std::isfinite(escale) && std::isfinite(ke0) && escale > 0.0; bool outer_converged = false; @@ -1306,11 +1168,9 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { Real momentum_error = std::numeric_limits::max(); for (outer_iter = 1; outer_iter <= outer_max; ++outer_iter) { - if (!solve_valid) - break; + if (!solve_valid) break; - const Real ke = - 0.5 * dens * (SQR(v[0]) + SQR(v[1]) + SQR(v[2])) / eref; + const Real ke = 0.5 * dens * (SQR(v[0]) + SQR(v[1]) + SQR(v[2])) / eref; std::array beta{v[0] / c, v[1] / c, v[2] / c}; const Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); if (!std::isfinite(beta2) || beta2 >= 1.0 || E <= 0.0) { @@ -1321,8 +1181,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const Real g = std::sqrt(g2); // F is normalized by c*eref, so F/E is the reduced flux. - const auto fedd = - EddingtonTensor({F[0] / E, F[1] / E, F[2] / E}); + const auto fedd = EddingtonTensor({F[0] / E, F[1] / E, F[2] / E}); const std::array bdp{ beta[0] * fedd[TensIdx::X11] + beta[1] * fedd[TensIdx::X12] + beta[2] * fedd[TensIdx::X13], @@ -1330,10 +1189,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { beta[2] * fedd[TensIdx::X23], beta[0] * fedd[TensIdx::X13] + beta[1] * fedd[TensIdx::X23] + beta[2] * fedd[TensIdx::X33]}; - const Real bdbdp = - beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; - const Real bdf = - beta[0] * F[0] + beta[1] * F[1] + beta[2] * F[2]; + const Real bdbdp = beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; + const Real bdf = beta[0] * F[0] + beta[1] * F[1] + beta[2] * F[2]; // Damped quasi-Newton solve for (B,E). Opacity derivatives are not // available, so accept only residual-decreasing trial steps. @@ -1341,20 +1198,15 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { int stalled_iterations = 0; bool inner_converged = false; for (inner_iter = 1; inner_iter <= inner_max; ++inner_iter) { - T = std::max( - tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); + T = std::max(tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); const Real eint = - dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / - eref; - const Real Cv = - dens * eos_d.SpecificHeatFromDensityTemperature(dens, T); + dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / eref; + const Real Cv = dens * eos_d.SpecificHeatFromDensityTemperature(dens, T); const Real fleck = FleckFactor(arad, T, Cv); - const Real sigp = - chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); + const Real sigp = chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); const Real sigs = - chat * dt * - scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); + chat * dt * scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); Real ca = 0.0; Real cb = 0.0; Real cd = 0.0; @@ -1362,15 +1214,14 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { sigp, sigs, g, g2, beta2, bdbdp, bdf, ca, cb, cd); Real Eeq_inner = E; const bool equilibrium_valid = - coeff_valid && ComputeCouplingEquilibriumEnergy( - E0, B, ca, cb, cd, Eeq_inner); + coeff_valid && + ComputeCouplingEquilibriumEnergy(E0, B, ca, cb, cd, Eeq_inner); // Retain the original residuals only to form the quasi-Newton // search direction. Use the conservation/equilibrium form below // for convergence and line-search acceptance. const Real G0 = ca * E - cb * B + cd; - const Real Fi = - (ke - ke0) + (eint - eg0) - c / chat * G0 - Q; + const Real Fi = (ke - ke0) + (eint - eg0) - c / chat * G0 - Q; const Real Fr = (E - E0) + G0; const Real conservation_inner = (ke - ke0) + (eint - eg0) + c / chat * (E - E0) - Q; @@ -1378,22 +1229,16 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const Real escale_inner = std::max( energy_floor_scale, - std::max( - std::max(std::abs(eg0), std::abs(eint)), - std::max(std::abs(Q), - std::max(c / chat * - (std::abs(E0) + std::abs(E)), - std::abs(ke - ke0))))); - inner_err = - std::max(std::abs(conservation_inner) / escale_inner, - c / chat * std::abs(source_inner) / escale_inner); - solve_valid = - solve_valid && equilibrium_valid && - std::isfinite(inner_err) && std::isfinite(fleck) && - std::isfinite(eint) && Cv > 0.0 && fleck >= 0.0 && - sigp >= 0.0 && sigs >= 0.0; - if (!solve_valid) - break; + std::max(std::max(std::abs(eg0), std::abs(eint)), + std::max(std::abs(Q), + std::max(c / chat * (std::abs(E0) + std::abs(E)), + std::abs(ke - ke0))))); + inner_err = std::max(std::abs(conservation_inner) / escale_inner, + c / chat * std::abs(source_inner) / escale_inner); + solve_valid = solve_valid && equilibrium_valid && std::isfinite(inner_err) && + std::isfinite(fleck) && std::isfinite(eint) && Cv > 0.0 && + fleck >= 0.0 && sigp >= 0.0 && sigs >= 0.0; + if (!solve_valid) break; if (inner_err <= inner_tol) { inner_converged = true; break; @@ -1404,8 +1249,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { else stalled_iterations = 0; previous_inner_err = inner_err; - if (stalled_iterations >= 8) - break; + if (stalled_iterations >= 8) break; const Real dfac = 1.0 + c / chat * fleck * cb; const Real denom = dfac + ca; @@ -1413,8 +1257,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { solve_valid = false; break; } - const Real dE = - dfac / denom * (-Fr) + fleck / denom * (-Fi * cb); + const Real dE = dfac / denom * (-Fr) + fleck / denom * (-Fi * cb); const Real dB = c / chat * fleck / denom * (-ca * Fr) + (1.0 + ca) * fleck / denom * (-Fi); @@ -1426,44 +1269,35 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { for (int ls = 0; ls < 12; ++ls) { const Real Etrial = std::max(efloor, E + alpha * dE); const Real Btrial = std::max(Bfloor, B + alpha * dB); - const Real Ttrial = std::max( - tfloor, - std::pow(std::max(eref * Btrial / arad, 0.0), 0.25)); + const Real Ttrial = + std::max(tfloor, std::pow(std::max(eref * Btrial / arad, 0.0), 0.25)); const Real eint_trial = - dens * - eos_d.InternalEnergyFromDensityTemperature(dens, Ttrial) / - eref; + dens * eos_d.InternalEnergyFromDensityTemperature(dens, Ttrial) / eref; const Real sigp_trial = - chat * dt * - opac_d.PlanckMeanAbsorptionCoefficient(dens, Ttrial); + chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, Ttrial); const Real sigs_trial = chat * dt * scat_d.RosselandMeanTotalScatteringCoefficient(dens, Ttrial); Real ca_trial = 0.0; Real cb_trial = 0.0; Real cd_trial = 0.0; - const bool coeff_trial_valid = - ComputeCouplingEnergyCoefficients( - sigp_trial, sigs_trial, g, g2, beta2, bdbdp, bdf, - ca_trial, cb_trial, cd_trial); + const bool coeff_trial_valid = ComputeCouplingEnergyCoefficients( + sigp_trial, sigs_trial, g, g2, beta2, bdbdp, bdf, ca_trial, cb_trial, + cd_trial); Real Eeq_trial = Etrial; const bool equilibrium_trial_valid = - coeff_trial_valid && ComputeCouplingEquilibriumEnergy( - E0, Btrial, ca_trial, cb_trial, - cd_trial, Eeq_trial); + coeff_trial_valid && + ComputeCouplingEquilibriumEnergy(E0, Btrial, ca_trial, cb_trial, + cd_trial, Eeq_trial); const Real conservation_trial = - (ke - ke0) + (eint_trial - eg0) + - c / chat * (Etrial - E0) - Q; + (ke - ke0) + (eint_trial - eg0) + c / chat * (Etrial - E0) - Q; const Real source_trial = Etrial - Eeq_trial; const Real escale_trial = std::max( energy_floor_scale, - std::max( - std::max(std::abs(eg0), std::abs(eint_trial)), - std::max( - std::abs(Q), - std::max(c / chat * - (std::abs(E0) + std::abs(Etrial)), - std::abs(ke - ke0))))); + std::max(std::max(std::abs(eg0), std::abs(eint_trial)), + std::max(std::abs(Q), + std::max(c / chat * (std::abs(E0) + std::abs(Etrial)), + std::abs(ke - ke0))))); const Real trial_err = std::max(std::abs(conservation_trial) / escale_trial, c / chat * std::abs(source_trial) / escale_trial); @@ -1477,9 +1311,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { best_trial_E = Etrial; best_trial_B = Btrial; } - if (trial_valid && - (trial_err <= inner_tol || - trial_err <= inner_err * (1.0 - 1.0e-4 * alpha))) { + if (trial_valid && (trial_err <= inner_tol || + trial_err <= inner_err * (1.0 - 1.0e-4 * alpha))) { E = Etrial; B = Btrial; accepted = true; @@ -1515,63 +1348,51 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { scalar_inner_used = true; int scalar_iterations = 0; const bool scalar_converged = SolveCouplingInnerScalar( - eg0, ke - ke0, Q, E0, efloor, Bfloor, dens, eref, arad, - tfloor, dt, chat, c, g, g2, beta2, bdbdp, bdf, - energy_floor_scale, inner_tol, eos_d, opac_d, scat_d, E, B, - inner_err, scalar_iterations); + eg0, ke - ke0, Q, E0, efloor, Bfloor, dens, eref, arad, tfloor, dt, chat, + c, g, g2, beta2, bdbdp, bdf, energy_floor_scale, inner_tol, eos_d, opac_d, + scat_d, E, B, inner_err, scalar_iterations); inner_iter += scalar_iterations; inner_converged = scalar_converged; } - if (!inner_converged && solve_valid && - momentum_predictor_count < outer_max) { + if (!inner_converged && solve_valid && momentum_predictor_count < outer_max) { // Bootstrap the kinetic-work term with the same coupled momentum // solve used below. Unlike the old full Picard predictor, this // cannot jump to a superluminal intermediate state when the // reduced-c radiation pseudo-inertia dominates the gas inertia. - const Real Tpred = std::max( - tfloor, - std::pow(std::max(eref * B / arad, 0.0), 0.25)); + const Real Tpred = + std::max(tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); const Real eg_pred = - dens * eos_d.InternalEnergyFromDensityTemperature(dens, - Tpred) / - eref; + dens * eos_d.InternalEnergyFromDensityTemperature(dens, Tpred) / eref; const Real dEg_pred = eg_pred - eg0; const Real sigp_pred = - chat * dt * - opac_d.RosselandMeanAbsorptionCoefficient(dens, Tpred); + chat * dt * opac_d.RosselandMeanAbsorptionCoefficient(dens, Tpred); const Real sigs_pred = - chat * dt * - scat_d.RosselandMeanTotalScatteringCoefficient(dens, Tpred); + chat * dt * scat_d.RosselandMeanTotalScatteringCoefficient(dens, Tpred); std::array beta_pred{v[0] / c, v[1] / c, v[2] / c}; std::array Fpred = F; Real Epred = E; Real dEkpred = dEk; Real predictor_error = std::numeric_limits::max(); int predictor_iterations = 0; - const bool predictor_converged = - SolveCouplingMomentumEnergy( - beta0, Fr0, E0, efloor, B, dEg_pred, Q, dens, eref, - c, chat, sigp_pred, sigs_pred, energy_floor_scale, - std::max(outer_tol, nonlinear_roundoff_tol), beta_pred, - Fpred, Epred, dEkpred, predictor_error, - predictor_iterations); + const bool predictor_converged = SolveCouplingMomentumEnergy( + beta0, Fr0, E0, efloor, B, dEg_pred, Q, dens, eref, c, chat, sigp_pred, + sigs_pred, energy_floor_scale, + std::max(outer_tol, nonlinear_roundoff_tol), beta_pred, Fpred, Epred, + dEkpred, predictor_error, predictor_iterations); momentum_iteration_count += predictor_iterations; if (predictor_converged) { - const Real predictor_scale = std::max( - energy_floor_scale, - std::max(E, std::max(CouplingFluxMagnitude(F), - CouplingFluxMagnitude(Fpred)))); + const Real predictor_scale = + std::max(energy_floor_scale, + std::max(E, std::max(CouplingFluxMagnitude(F), + CouplingFluxMagnitude(Fpred)))); Real predictor_change = std::abs(Epred - E) / predictor_scale; for (int d = 0; d < 3; ++d) { - predictor_change = std::max( - predictor_change, - std::abs(Fpred[d] - F[d]) / predictor_scale); - predictor_change = std::max( - predictor_change, - std::abs(beta_pred[d] - v[d] / c)); + predictor_change = std::max(predictor_change, + std::abs(Fpred[d] - F[d]) / predictor_scale); + predictor_change = + std::max(predictor_change, std::abs(beta_pred[d] - v[d] / c)); } - if (predictor_change > - 8.0 * std::numeric_limits::epsilon()) { + if (predictor_change > 8.0 * std::numeric_limits::epsilon()) { E = Epred; F = Fpred; dEk = dEkpred; @@ -1586,24 +1407,19 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { } } } - if (!inner_converged) - break; + if (!inner_converged) break; - T = std::max( - tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); + T = std::max(tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); const Real eg = - dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / - eref; + dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / eref; dEg = eg - eg0; const Real sigp_flux = - chat * dt * - opac_d.RosselandMeanAbsorptionCoefficient(dens, T); + chat * dt * opac_d.RosselandMeanAbsorptionCoefficient(dens, T); const Real sigs_flux = - chat * dt * - scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); - if (!std::isfinite(sigp_flux) || !std::isfinite(sigs_flux) || - sigp_flux < 0.0 || sigs_flux < 0.0) { + chat * dt * scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); + if (!std::isfinite(sigp_flux) || !std::isfinite(sigs_flux) || sigp_flux < 0.0 || + sigs_flux < 0.0) { solve_valid = false; break; } @@ -1614,16 +1430,12 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // mu = eref/(rho*c*chat) // is large; that is exactly the regime encountered just above a // density floor. The local 3D Newton solve removes this stiffness. - std::array beta_momentum{v[0] / c, v[1] / c, - v[2] / c}; + std::array beta_momentum{v[0] / c, v[1] / c, v[2] / c}; int momentum_iterations = 0; - const bool momentum_converged = - SolveCouplingMomentumEnergy( - beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, - chat, sigp_flux, sigs_flux, energy_floor_scale, - std::max(outer_tol, nonlinear_roundoff_tol), - beta_momentum, F, E, dEk, momentum_error, - momentum_iterations); + const bool momentum_converged = SolveCouplingMomentumEnergy( + beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp_flux, + sigs_flux, energy_floor_scale, std::max(outer_tol, nonlinear_roundoff_tol), + beta_momentum, F, E, dEk, momentum_error, momentum_iterations); momentum_iteration_count += momentum_iterations; if (!momentum_converged) { solve_valid = false; @@ -1640,121 +1452,97 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // outer iteration. This avoids declaring convergence merely because // the kinetic-energy correction changed little in a Keplerian flow. const std::array beta_out{v[0] / c, v[1] / c, v[2] / c}; - const Real beta2_out = SQR(beta_out[0]) + SQR(beta_out[1]) + - SQR(beta_out[2]); + const Real beta2_out = SQR(beta_out[0]) + SQR(beta_out[1]) + SQR(beta_out[2]); if (!std::isfinite(beta2_out) || beta2_out >= 1.0 || E <= 0.0) { solve_valid = false; break; } const Real g2_out = 1.0 / (1.0 - beta2_out); const Real g_out = std::sqrt(g2_out); - const auto fedd_out = EddingtonTensor( - {F[0] / E, F[1] / E, F[2] / E}); - const std::array bdp_out{ - beta_out[0] * fedd_out[TensIdx::X11] + - beta_out[1] * fedd_out[TensIdx::X12] + - beta_out[2] * fedd_out[TensIdx::X13], - beta_out[0] * fedd_out[TensIdx::X12] + - beta_out[1] * fedd_out[TensIdx::X22] + - beta_out[2] * fedd_out[TensIdx::X23], - beta_out[0] * fedd_out[TensIdx::X13] + - beta_out[1] * fedd_out[TensIdx::X23] + - beta_out[2] * fedd_out[TensIdx::X33]}; - const Real bdbdp_out = beta_out[0] * bdp_out[0] + - beta_out[1] * bdp_out[1] + - beta_out[2] * bdp_out[2]; - const Real bdf_out = beta_out[0] * F[0] + beta_out[1] * F[1] + - beta_out[2] * F[2]; - - T = std::max( - tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); + const auto fedd_out = EddingtonTensor({F[0] / E, F[1] / E, F[2] / E}); + const std::array bdp_out{beta_out[0] * fedd_out[TensIdx::X11] + + beta_out[1] * fedd_out[TensIdx::X12] + + beta_out[2] * fedd_out[TensIdx::X13], + beta_out[0] * fedd_out[TensIdx::X12] + + beta_out[1] * fedd_out[TensIdx::X22] + + beta_out[2] * fedd_out[TensIdx::X23], + beta_out[0] * fedd_out[TensIdx::X13] + + beta_out[1] * fedd_out[TensIdx::X23] + + beta_out[2] * fedd_out[TensIdx::X33]}; + const Real bdbdp_out = beta_out[0] * bdp_out[0] + beta_out[1] * bdp_out[1] + + beta_out[2] * bdp_out[2]; + const Real bdf_out = + beta_out[0] * F[0] + beta_out[1] * F[1] + beta_out[2] * F[2]; + + T = std::max(tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); const Real eg_out = - dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / - eref; + dens * eos_d.InternalEnergyFromDensityTemperature(dens, T) / eref; dEg = eg_out - eg0; const Real sigp_energy = chat * dt * opac_d.PlanckMeanAbsorptionCoefficient(dens, T); const Real sigs_out = - chat * dt * - scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); + chat * dt * scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); Real ca_out = 0.0; Real cb_out = 0.0; Real cd_out = 0.0; Real Eeq_out = E; const bool energy_coeff_valid = ComputeCouplingEnergyCoefficients( - sigp_energy, sigs_out, g_out, g2_out, beta2_out, - bdbdp_out, bdf_out, ca_out, cb_out, cd_out); + sigp_energy, sigs_out, g_out, g2_out, beta2_out, bdbdp_out, bdf_out, ca_out, + cb_out, cd_out); const bool equilibrium_valid = - energy_coeff_valid && ComputeCouplingEquilibriumEnergy( - E0, B, ca_out, cb_out, cd_out, - Eeq_out); - const Real conservation_residual = - dEg + dEk + c / chat * dEr - Q; + energy_coeff_valid && + ComputeCouplingEquilibriumEnergy(E0, B, ca_out, cb_out, cd_out, Eeq_out); + const Real conservation_residual = dEg + dEk + c / chat * dEr - Q; const Real source_residual = E - Eeq_out; const Real sigp_flux_out = - chat * dt * - opac_d.RosselandMeanAbsorptionCoefficient(dens, T); + chat * dt * opac_d.RosselandMeanAbsorptionCoefficient(dens, T); const Real sigf_flux_out = sigp_flux_out + sigs_out; const Real a_out = g_out * sigf_flux_out; const Real b_out = 2.0 * g2_out * g_out * sigs_out; const Real d1_out = - g_out * - (sigp_flux_out * B + - g2_out * sigs_out * (1.0 + bdbdp_out) * E); + g_out * (sigp_flux_out * B + g2_out * sigs_out * (1.0 + bdbdp_out) * E); const Real d2_out = g_out * sigf_flux_out * E; const std::array rhs_out{ Fr0[0] + d1_out * beta_out[0] + d2_out * bdp_out[0], Fr0[1] + d1_out * beta_out[1] + d2_out * bdp_out[1], Fr0[2] + d1_out * beta_out[2] + d2_out * bdp_out[2]}; - const auto Fsolve_out = - SolveRadFlux(1.0 + a_out, b_out, beta_out, rhs_out); + const auto Fsolve_out = SolveRadFlux(1.0 + a_out, b_out, beta_out, rhs_out); const auto Ftarget_out = ProjectCouplingFluxToEnergy(Fsolve_out, E); - escale = std::max( - energy_floor_scale, - std::max( - std::max(std::abs(eg0), std::abs(eg_out)), - std::max( - std::abs(Q), - std::max(c / chat * - (std::abs(E0) + std::abs(E)), - std::abs(dEk))))); + escale = + std::max(energy_floor_scale, + std::max(std::max(std::abs(eg0), std::abs(eg_out)), + std::max(std::abs(Q), + std::max(c / chat * (std::abs(E0) + std::abs(E)), + std::abs(dEk))))); const Real energy_residual = std::max(std::abs(conservation_residual) / escale, c / chat * std::abs(source_residual) / escale); const Real flux_scale = std::max( energy_floor_scale, - std::max( - E, std::max(CouplingFluxMagnitude(Fr0), - std::max(CouplingFluxMagnitude(F), - CouplingFluxMagnitude(Ftarget_out))))); + std::max(E, std::max(CouplingFluxMagnitude(Fr0), + std::max(CouplingFluxMagnitude(F), + CouplingFluxMagnitude(Ftarget_out))))); Real flux_residual = 0.0; for (int d = 0; d < 3; ++d) { flux_residual = - std::max(flux_residual, - std::abs(F[d] - Ftarget_out[d]) / flux_scale); + std::max(flux_residual, std::abs(F[d] - Ftarget_out[d]) / flux_scale); } - const Real realizability_residual = - std::max(0.0, (CouplingFluxMagnitude(F) - E) / - std::max(E, energy_floor_scale)); + const Real realizability_residual = std::max( + 0.0, (CouplingFluxMagnitude(F) - E) / std::max(E, energy_floor_scale)); // These are residuals of the coupled equations evaluated at the // current state. Iterate-to-iterate changes are deliberately not an // acceptance criterion: a converged nonlinear root need not move by // less than a user tolerance tighter than the local arithmetic can // resolve. - outer_err = std::max( - energy_residual, - std::max(flux_residual, realizability_residual)); - solve_valid = - solve_valid && equilibrium_valid && - std::isfinite(outer_err) && std::isfinite(dEr) && - std::isfinite(sigp_energy) && std::isfinite(sigs_out) && - std::isfinite(sigp_flux_out) && - sigp_energy >= 0.0 && sigs_out >= 0.0 && - sigp_flux_out >= 0.0; - if (!solve_valid) - break; + outer_err = + std::max(energy_residual, std::max(flux_residual, realizability_residual)); + solve_valid = solve_valid && equilibrium_valid && std::isfinite(outer_err) && + std::isfinite(dEr) && std::isfinite(sigp_energy) && + std::isfinite(sigs_out) && std::isfinite(sigp_flux_out) && + sigp_energy >= 0.0 && sigs_out >= 0.0 && sigp_flux_out >= 0.0; + if (!solve_valid) break; have_complete_iterate = true; if (outer_err < best_outer_err) { @@ -1778,8 +1566,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // energy-only solve below with momentum frozen, rather than repeatedly // applying a radiation drag whose kinetic-energy work cannot be paid by // either the gas or radiation field. - const Real floor_active_tol = - 256.0 * std::numeric_limits::epsilon(); + const Real floor_active_tol = 256.0 * std::numeric_limits::epsilon(); const bool material_floor_active = B <= Bfloor * (1.0 + floor_active_tol) + Fuzz(); const bool radiation_floor_active = @@ -1810,14 +1597,11 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { dEk = 0.0; const Real Emin = std::max(efloor, CouplingFluxMagnitude(Fr0)); const Real eg_floor = - dens * eos_d.InternalEnergyFromDensityTemperature(dens, - tfloor) / - eref; - const Real eg_max = - eg0 + Q + c / chat * (E0 - Emin); - - bool fallback_valid = std::isfinite(eg_floor) && - std::isfinite(eg_max) && eg_max >= eg_floor; + dens * eos_d.InternalEnergyFromDensityTemperature(dens, tfloor) / eref; + const Real eg_max = eg0 + Q + c / chat * (E0 - Emin); + + bool fallback_valid = + std::isfinite(eg_floor) && std::isfinite(eg_max) && eg_max >= eg_floor; bool fallback_have = false; bool bracketed = false; Real bracket_lo = eg_floor; @@ -1834,18 +1618,15 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { if (fallback_valid) { constexpr int fallback_scan_points = 64; for (int n = 0; n <= fallback_scan_points; ++n) { - const Real frac = - static_cast(n) / fallback_scan_points; - const Real eg_trial = - eg_floor + frac * (eg_max - eg_floor); + const Real frac = static_cast(n) / fallback_scan_points; + const Real eg_trial = eg_floor + frac * (eg_max - eg_floor); Real E_trial = E0; Real B_trial = B; Real residual = 0.0; const bool valid_trial = EvaluateCouplingEnergyOnlyResidual( - eg_trial, eg0, Q, E0, Emin, dens, eref, arad, tfloor, dt, - chat, c, eos_d, opac_d, E_trial, B_trial, residual); - if (!valid_trial) - continue; + eg_trial, eg0, Q, E0, Emin, dens, eref, arad, tfloor, dt, chat, c, + eos_d, opac_d, E_trial, B_trial, residual); + if (!valid_trial) continue; fallback_have = true; if (std::abs(residual) < best_residual) { @@ -1875,10 +1656,9 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { Real B_mid = B; Real residual_mid = 0.0; const bool valid_mid = EvaluateCouplingEnergyOnlyResidual( - eg_mid, eg0, Q, E0, Emin, dens, eref, arad, tfloor, dt, - chat, c, eos_d, opac_d, E_mid, B_mid, residual_mid); - if (!valid_mid) - break; + eg_mid, eg0, Q, E0, Emin, dens, eref, arad, tfloor, dt, chat, c, eos_d, + opac_d, E_mid, B_mid, residual_mid); + if (!valid_mid) break; if (std::abs(residual_mid) < best_residual) { best_residual = std::abs(residual_mid); best_eg = eg_mid; @@ -1894,8 +1674,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { } if ((bracket_hi - bracket_lo) / std::max(energy_floor_scale, - std::abs(bracket_hi) + - std::abs(bracket_lo)) <= + std::abs(bracket_hi) + std::abs(bracket_lo)) <= nonlinear_roundoff_tol) break; } @@ -1908,10 +1687,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { dEr = E - E0; const Real fallback_scale = std::max( energy_floor_scale, - std::max(std::abs(Q), - std::max(std::abs(eg0) + std::abs(best_eg), - c / chat * - (std::abs(E0) + std::abs(E))))); + std::max(std::abs(Q), std::max(std::abs(eg0) + std::abs(best_eg), + c / chat * (std::abs(E0) + std::abs(E))))); fallback_error = c / chat * best_residual / std::max(fallback_scale, Fuzz()); fallback_success = std::isfinite(fallback_error); @@ -1928,15 +1705,14 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const bool accepted_floor_fallback = floor_constrained_failure && fallback_success; - if (!outer_converged && fatal_if_unconverged && - !numerical_atmosphere && !accepted_floor_fallback) { - const Real beta_fail = - std::sqrt(SQR(v[0] / c) + SQR(v[1] / c) + SQR(v[2] / c)); + if (!outer_converged && fatal_if_unconverged && !numerical_atmosphere && + !accepted_floor_fallback) { + const Real beta_fail = std::sqrt(SQR(v[0] / c) + SQR(v[1] / c) + SQR(v[2] / c)); const Real fred_fail = CouplingFluxMagnitude(F) / std::max(E, energy_floor_scale); const Real mu_fail = eref / (dens * c * chat); - const Real delta_f_fail = CouplingFluxMagnitude( - {F[0] - Fr0[0], F[1] - Fr0[1], F[2] - Fr0[2]}); + const Real delta_f_fail = + CouplingFluxMagnitude({F[0] - Fr0[0], F[1] - Fr0[1], F[2] - Fr0[2]}); printf("MatterCoupling full fail (%d,%d,%d,%d): outer=%.17e " "(best=%.17e, tol=%.17e), inner=%.17e (tol=%.17e), " "fallback=%.17e, E=%.17e, E0=%.17e, efloor=%.17e, " @@ -1947,14 +1723,13 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { "inner_iter=%d, momentum_iter=%d, " "momentum_err=%.17e, valid=%d, scalar=%d, " "predictors=%d, complete=%d, floor_constrained=%d\n", - b, k, j, i, outer_err, best_outer_err, outer_tol, inner_err, - inner_tol, fallback_error, E, E0, efloor, B, Bfloor, Q, - dEk, dens, dens / dflr, std::sqrt(beta20), beta_fail, - fred_fail, delta_f_fail, chat / c, mu_fail, + b, k, j, i, outer_err, best_outer_err, outer_tol, inner_err, inner_tol, + fallback_error, E, E0, efloor, B, Bfloor, Q, dEk, dens, dens / dflr, + std::sqrt(beta20), beta_fail, fred_fail, delta_f_fail, chat / c, mu_fail, static_cast(mu_fail <= 1.0), outer_iter, inner_iter, - momentum_iteration_count, momentum_error, - solve_valid, scalar_inner_used, momentum_predictor_count, - have_complete_iterate, floor_constrained_failure); + momentum_iteration_count, momentum_error, solve_valid, scalar_inner_used, + momentum_predictor_count, have_complete_iterate, + floor_constrained_failure); PARTHENON_FAIL("Outer not converged"); } From 747d992e4beb5f032cdee2f95594b7284da61db2 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Mon, 13 Jul 2026 11:59:54 -0600 Subject: [PATCH 06/14] Split and Unsplit moments. --- src/artemis_driver.cpp | 43 +++++++++++----- src/artemis_driver.hpp | 2 +- src/radiation/moments/moments.cpp | 85 ++++++++++++++++++++++++++----- 3 files changed, 101 insertions(+), 29 deletions(-) diff --git a/src/artemis_driver.cpp b/src/artemis_driver.cpp index a4af03db..c3d55cf6 100644 --- a/src/artemis_driver.cpp +++ b/src/artemis_driver.cpp @@ -83,11 +83,17 @@ ArtemisDriver::ArtemisDriver(ParameterInput *pin, ApplicationInput *app_in ndim = pm->ndim; // Moments integrator + do_moment_unsplit = false; + do_moment_split = false; if (do_moment) { auto rad_int = pin->GetOrAddString("radiation/moment", "integrator", "rk2"); PARTHENON_REQUIRE(((rad_int == "rk1") || (rad_int == "rk2") || (rad_int == "rk3")), "radiation/integrator must be rk1,rk2, or rk3.") - rad_integrator = std::make_unique(rad_int); + do_moment_split = pin->GetOrAddBoolean("radiation/moment", "substep", false); + do_moment_unsplit = !do_moment_split; + if (do_moment_split) { + rad_integrator = std::make_unique(rad_int); + } } // NBody integrator and initialization @@ -153,7 +159,7 @@ TaskListStatus ArtemisDriver::Step() { } // Operator split, moments subcycling (M1 or P1) - if (do_moment) { + if (do_moment_split) { status = Moments::MomentsDriver(pmesh, tm, rad_integrator.get()); if (status != TaskListStatus::complete) return status; } @@ -180,7 +186,7 @@ void ArtemisDriver::PreStepTasks() { // set the integration timestep integrator->dt = tm.dt; if (do_nbody) nbody_integrator->dt = tm.dt; - if (do_moment) rad_integrator->dt = tm.dt; + if (do_moment_split) rad_integrator->dt = tm.dt; // Extract Base MeshData Registers auto &base = pmesh->mesh_data.Get(); @@ -193,7 +199,7 @@ void ArtemisDriver::PreStepTasks() { auto &u1 = pmesh->mesh_data.Add("u1", u0); // Assign registers with fields required for moments - if (do_moment) { + if (do_moment_split) { parthenon::Metadata::FlagCollection moments_flags, geom_flags; moments_flags.TakeUnion(pmesh->packages.Get("moments")->GetMetadataFlag()); geom_flags.TakeUnion(pmesh->packages.Get("geometry")->GetMetadataFlag()); @@ -258,10 +264,13 @@ TaskCollection ArtemisDriver::StepTasks() { // Compute hydrodynamic fluxes // NOTE(@adempsey): 1st stage of VL2 uses piecewise constant reconstruction const bool do_pcm = ((stage == 1) && (integrator->GetName() == "vl2")); - TaskID gas_flx = none, dust_flx = none; + TaskID gas_flx = none, dust_flx = none, rad_flx = none; if (do_gas && update_fluxes) gas_flx = tl.AddTask(none, Gas::CalculateFluxes, u0.get(), do_pcm); if (do_dust) dust_flx = tl.AddTask(none, Dust::CalculateFluxes, u0.get(), do_pcm); + if (do_moment_unsplit) { + rad_flx = tl.AddTask(none, Moments::CalculateFluxes, u0.get()); + } // Compute (gas) diffusive fluxes TaskID diff_flx = none; @@ -275,20 +284,23 @@ TaskCollection ArtemisDriver::StepTasks() { // Communicate and set fluxes auto send_flx = - tl.AddTask(gas_flx | dust_flx | diff_flx, + tl.AddTask(gas_flx | dust_flx | rad_flx | diff_flx, parthenon::SendBoundBufs, u0); auto recv_flx = tl.AddTask(start_flx_recv, parthenon::ReceiveFluxCorrections, u0); auto set_flx = tl.AddTask(recv_flx, parthenon::SetFluxCorrections, u0); // Apply flux divergence auto update = - tl.AddTask(gas_flx | dust_flx | set_flx, ArtemisUtils::ApplyUpdate, - u0.get(), u1.get(), g0, g1, bdt); + tl.AddTask(gas_flx | dust_flx | rad_flx | set_flx, + ArtemisUtils::ApplyUpdate, u0.get(), u1.get(), g0, g1, bdt); // Apply "coordinate source terms" - TaskID gas_coord_src = update, dust_coord_src = update; + TaskID gas_coord_src = update, dust_coord_src = update, rad_coord_src; if (do_gas) gas_coord_src = tl.AddTask(update, Gas::FluxSource, u0.get(), bdt); if (do_dust) dust_coord_src = tl.AddTask(update, Dust::FluxSource, u0.get(), bdt); + if (do_moment_unsplit) { + rad_coord_src = tl.AddTask(update, Moments::FluxSource, u0.get(), bdt); + } // Apply (gas) diffusion sources TaskID gas_diff_src = gas_coord_src | diff_flx | set_flx; @@ -311,17 +323,20 @@ TaskCollection ArtemisDriver::StepTasks() { tl.AddTask(gravity_src, SelfGravity::SelfGravity, u0.get(), time, bdt); } - TaskID rt_src = self_gravity_src; + TaskID rad_src = self_gravity_src; // Note that radiation moments will handle this source term if active - if (do_raytrace && !do_moment) { - rt_src = tl.AddTask(self_gravity_src, Gas::DepositEnergy, u0.get(), bdt); + if (do_moment_unsplit) { + rad_src = + tl.AddTask(self_gravity_src, Moments::MatterCoupling, u0.get(), bdt); + } else if (do_raytrace && !do_moment) { + rad_src = tl.AddTask(self_gravity_src, Gas::DepositEnergy, u0.get(), bdt); } // Apply rotating frame source term - TaskID rframe_src = rt_src; + TaskID rframe_src = rad_src; if (do_rotating_frame || do_orbital_advection) { rframe_src = - tl.AddTask(rt_src, RotatingFrame::RotatingFrameForce, u0.get(), time, bdt); + tl.AddTask(rad_src, RotatingFrame::RotatingFrameForce, u0.get(), time, bdt); } // Apply drag source term diff --git a/src/artemis_driver.hpp b/src/artemis_driver.hpp index 0ac79ce4..e08edeeb 100644 --- a/src/artemis_driver.hpp +++ b/src/artemis_driver.hpp @@ -53,7 +53,7 @@ class ArtemisDriver : public EvolutionDriver { IntegratorPtr_t integrator, nbody_integrator, rad_integrator; StateDescriptor *artemis_pkg; int ndim; - bool do_gas, do_dust, do_moment, do_imc; + bool do_gas, do_dust, do_imc, do_moment, do_moment_split, do_moment_unsplit; bool do_gravity, do_self_gravity, do_nbody, do_rotating_frame; bool do_orbital_advection; bool do_cooling, do_drag, do_viscosity, do_conduction, do_diffusion; diff --git a/src/radiation/moments/moments.cpp b/src/radiation/moments/moments.cpp index 1e7d7a66..e8f99b91 100644 --- a/src/radiation/moments/moments.cpp +++ b/src/radiation/moments/moments.cpp @@ -127,7 +127,7 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("inner_iteration_tol", pin->GetOrAddReal("radiation/moment", "inner_iteration_tol", 1e-10)); - const bool substep = pin->GetOrAddBoolean("radiation/moment", "substep", true); + const bool substep = pin->GetOrAddBoolean("radiation/moment", "substep", false); if (!substep) { if (coords == Coordinates::cartesian) { moments->EstimateTimestepMesh = EstimateTimeStepMesh; @@ -162,27 +162,85 @@ std::shared_ptr Initialize(ParameterInput *pin, // Control field for sparse radiation fields std::string control_field = rad::cons::energy::name(); - // Conserved Energy Density - Metadata m = Metadata({Metadata::Cell, Metadata::Conserved, Metadata::Independent, + auto mflags_cons = [&MetadataMoments, &MetadataOperatorSplit, + &substep](const int size) { + if (size == 1) { + if (substep) { + return Metadata({Metadata::Cell, Metadata::Conserved, Metadata::Independent, Metadata::WithFluxes, Metadata::Sparse, MetadataMoments, MetadataOperatorSplit}); + } else { + + return Metadata({Metadata::Cell, Metadata::Conserved, Metadata::Independent, + Metadata::WithFluxes, Metadata::Sparse}); + } + + } else { + if (substep) { + return Metadata({Metadata::Cell, Metadata::Vector, Metadata::Conserved, + Metadata::Independent, Metadata::WithFluxes, Metadata::Sparse, + MetadataMoments, MetadataOperatorSplit}, + std::vector({size})); + } else { + + return Metadata({Metadata::Cell, Metadata::Vector, Metadata::Conserved, + Metadata::Independent, Metadata::WithFluxes, Metadata::Sparse}, + std::vector({size})); + } + } + }; + + auto mflags_prim = [&MetadataMoments, &MetadataOperatorSplit, + &substep](const int size) { + if (size == 1) { + if (substep) { + return Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, + Metadata::OneCopy, Metadata::FillGhost, Metadata::Sparse, + MetadataMoments, MetadataOperatorSplit}); + } else { + return Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, + Metadata::OneCopy, Metadata::FillGhost, Metadata::Sparse}); + } + + } else { + if (substep) { + return Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, + Metadata::OneCopy, Metadata::FillGhost, Metadata::Sparse, + MetadataMoments, MetadataOperatorSplit}, + std::vector({size})); + } else { + return Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, + Metadata::OneCopy, Metadata::FillGhost, Metadata::Sparse}, + std::vector({size})); + } + } + }; + auto mflags_prim_withflux = [&MetadataMoments, &MetadataOperatorSplit, &substep]() { + if (substep) { + return Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, + Metadata::WithFluxes, Metadata::Sparse, MetadataMoments, + MetadataOperatorSplit}); + } else { + return Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, + Metadata::WithFluxes, Metadata::Sparse}); + } + }; + + // Conserved Energy Density + Metadata m = mflags_cons(1); + ArtemisUtils::EnrollArtemisRefinementOps(m, coords, log); m.SetSparseThresholds(0.0, 0.0, 0.0); moments->AddSparsePool(m, control_field, fluidids); // Conserved Flux - m = Metadata({Metadata::Cell, Metadata::Vector, Metadata::Conserved, - Metadata::Independent, Metadata::WithFluxes, Metadata::Sparse, - MetadataMoments, MetadataOperatorSplit}, - std::vector({3})); + m = mflags_cons(3); ArtemisUtils::EnrollArtemisRefinementOps(m, coords, log); m.SetSparseThresholds(0.0, 0.0, 0.0); moments->AddSparsePool(m, control_field, fluidids); // Primitive Energy Density - m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, Metadata::OneCopy, - Metadata::FillGhost, Metadata::Sparse, MetadataMoments, - MetadataOperatorSplit}); + m = mflags_prim(1); ArtemisUtils::EnrollArtemisRefinementOps(m, coords, log); m.SetSparseThresholds(0.0, 0.0, 0.0); moments->AddSparsePool(m, control_field, fluidids); @@ -191,15 +249,14 @@ std::shared_ptr Initialize(ParameterInput *pin, m = Metadata({Metadata::Cell, Metadata::Derived, Metadata::Intensive, Metadata::OneCopy, Metadata::WithFluxes, Metadata::Sparse, MetadataMoments, MetadataOperatorSplit}); + m = mflags_prim_withflux(); + ArtemisUtils::EnrollArtemisRefinementOps(m, coords, log); m.SetSparseThresholds(0.0, 0.0, 0.0); moments->AddSparsePool(m, control_field, fluidids); // Primitive Reduced Flux - m = Metadata({Metadata::Cell, Metadata::Vector, Metadata::Derived, Metadata::Intensive, - Metadata::OneCopy, Metadata::FillGhost, Metadata::Sparse, MetadataMoments, - MetadataOperatorSplit}, - std::vector({3})); + m = mflags_prim(3); ArtemisUtils::EnrollArtemisRefinementOps(m, coords, log); m.SetSparseThresholds(0.0, 0.0, 0.0); moments->AddSparsePool(m, control_field, fluidids); From 67dd42eb8c679fce9e9d5117b88363ad65c9e05e Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Tue, 14 Jul 2026 08:40:15 -0600 Subject: [PATCH 07/14] Add Eps function --- src/artemis.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/artemis.hpp b/src/artemis.hpp index 20e13e93..42be71d4 100644 --- a/src/artemis.hpp +++ b/src/artemis.hpp @@ -191,6 +191,10 @@ KOKKOS_FORCEINLINE_FUNCTION constexpr auto Big() { return std::numeric_limits::max(); } template +KOKKOS_FORCEINLINE_FUNCTION constexpr auto Eps() { + return std::numeric_limits::epsilon(); +} +template KOKKOS_FORCEINLINE_FUNCTION constexpr auto Tiny() { return std::numeric_limits::lowest(); } From 898c3d34cee9ee6b5303f8f9058a6d516e30f5dd Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Tue, 14 Jul 2026 08:40:26 -0600 Subject: [PATCH 08/14] Add VNorm --- src/utils/artemis_utils.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/utils/artemis_utils.hpp b/src/utils/artemis_utils.hpp index c81a797b..7090e7ae 100644 --- a/src/utils/artemis_utils.hpp +++ b/src/utils/artemis_utils.hpp @@ -43,6 +43,13 @@ KOKKOS_FORCEINLINE_FUNCTION Real VDot(const V1 &a, const V2 &b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } +KOKKOS_FORCEINLINE_FUNCTION Real VNorm(const std::array &v) { + return std::sqrt(SQR(v[0]) + SQR(v[1])); +} +KOKKOS_FORCEINLINE_FUNCTION Real VNorm(const std::array &v) { + return std::sqrt(SQR(v[0]) + SQR(v[1]) + SQR(v[2])); +} + //---------------------------------------------------------------------------------------- //! \fn Real ArtemisUtils::DualEnergySIE(vmesh, const int b, const int n, const int k, //! const int j, const int i, const Real de_switch, From 1fb23695a9b97016844b8a4196a8dc29c800d420 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Tue, 14 Jul 2026 08:40:55 -0600 Subject: [PATCH 09/14] Cleanup and removal of unneeded code --- src/radiation/moments/matter_coupling.hpp | 359 +++++++++------------- 1 file changed, 149 insertions(+), 210 deletions(-) diff --git a/src/radiation/moments/matter_coupling.hpp b/src/radiation/moments/matter_coupling.hpp index a7e8660e..f23002aa 100644 --- a/src/radiation/moments/matter_coupling.hpp +++ b/src/radiation/moments/matter_coupling.hpp @@ -13,6 +13,7 @@ // the public, perform publicly and display publicly, and to permit others to do // so. //======================================================================================== + // This file was modified in part with the assistance of generative AI. #ifndef RADIATION_MOMENTS_MATTER_COUPLING_HPP_ #define RADIATION_MOMENTS_MATTER_COUPLING_HPP_ @@ -29,18 +30,11 @@ using ArtemisUtils::EOS; using ArtemisUtils::MeanOpacity; using ArtemisUtils::MeanScattering; using ArtemisUtils::VI; +using ArtemisUtils::VNorm; namespace Moments { -KOKKOS_INLINE_FUNCTION Real FiniteOrZero(const Real x) { - return std::isfinite(x) ? x : 0.0; -} - -KOKKOS_INLINE_FUNCTION Real CouplingFluxMagnitude(const std::array &F) { - return std::sqrt(SQR(F[0]) + SQR(F[1]) + SQR(F[2])); -} - -KOKKOS_INLINE_FUNCTION bool +KOKKOS_INLINE_FUNCTION void ComputeCouplingEnergyCoefficients(const Real sigp, const Real sigs, const Real g, const Real g2, const Real beta2, const Real bdbdp, const Real bdf, Real &ca, Real &cb, Real &cd) { @@ -52,7 +46,6 @@ ComputeCouplingEnergyCoefficients(const Real sigp, const Real sigs, const Real g ca = g * (sigp - g2 * sigs * (beta2 + bdbdp)); cb = g * sigp; cd = g * bdf * (g2 * sigs * (1.0 + beta2) - sigp); - return std::isfinite(ca) && std::isfinite(cb) && std::isfinite(cd); } KOKKOS_INLINE_FUNCTION bool ComputeCouplingEquilibriumEnergy(const Real E0, const Real B, @@ -64,17 +57,17 @@ KOKKOS_INLINE_FUNCTION bool ComputeCouplingEquilibriumEnergy(const Real E0, cons // the optically thick limit, where ca and cb may be O(1e10) and directly // evaluating ca*E - cb*B loses all useful digits near LTE. const Real denom = 1.0 + ca; - if (!std::isfinite(denom) || std::abs(denom) <= Fuzz()) return false; + if (std::abs(denom) <= Fuzz()) return false; const Real inv_denom = 1.0 / denom; Eeq = E0 * inv_denom + (cb * inv_denom) * B - cd * inv_denom; - return std::isfinite(Eeq); + return true; } KOKKOS_INLINE_FUNCTION std::array ProjectCouplingFluxToEnergy(const std::array &Fin, const Real E) { - const Real fmag = CouplingFluxMagnitude(Fin); - if (!std::isfinite(fmag) || !std::isfinite(E) || E <= 0.0) return {0.0, 0.0, 0.0}; - if (fmag <= E || fmag <= Fuzz()) return Fin; + const Real fmag = VNorm(Fin); + if (E <= 0.0) return {0.0, 0.0, 0.0}; + if (fmag <= E || fmag <= Fuzz()) return Fin; const Real scale = E / fmag; return {scale * Fin[0], scale * Fin[1], scale * Fin[2]}; } @@ -94,7 +87,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingDense3x3(Real A[3][3], Real rhs[3], pivot_abs = candidate; } } - if (!std::isfinite(pivot_abs) || pivot_abs <= Fuzz()) return false; + if (pivot_abs <= Fuzz()) return false; if (pivot != col) { for (int j = col; j < 3; ++j) { const Real tmp = A[col][j]; @@ -120,10 +113,8 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingDense3x3(Real A[3][3], Real rhs[3], Real value = rhs[row]; for (int j = row + 1; j < 3; ++j) value -= A[row][j] * x[j]; - if (!std::isfinite(A[row][row]) || std::abs(A[row][row]) <= Fuzz()) - return false; + if (std::abs(A[row][row]) <= Fuzz()) return false; x[row] = value / A[row][row]; - if (!std::isfinite(x[row])) return false; } return true; } @@ -141,10 +132,9 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualBeta( // where F is normalized by c*eref. This relation is the source of the // instability in a direct Picard update when eref/(rho*c*chat) is large. const Real mu = eref / (dens * c * chat); - if (!std::isfinite(mu) || mu <= 0.0) return false; const Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); - if (!std::isfinite(beta2) || beta2 >= 1.0) return false; + if (beta2 >= 1.0) return false; std::array delta_F{0.0, 0.0, 0.0}; Real delta_F2 = 0.0; @@ -161,15 +151,14 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualBeta( // dEk = (c/chat)[-beta0.dF + 0.5*mu*|dF|^2]. dEk = c / chat * (-beta0_dot_delta_F + 0.5 * mu * delta_F2); E = E0 + chat / c * (Q - dEg - dEk); - const Real floor_slop = 128.0 * std::numeric_limits::epsilon() * - std::max(1.0, std::max(std::abs(E0), std::abs(E))); - if (!std::isfinite(E) || E < efloor - floor_slop || E <= 0.0) return false; + const Real floor_slop = + 128.0 * Eps() * std::max(1.0, std::max(std::abs(E0), std::abs(E))); + if (E < efloor - floor_slop || E <= 0.0) return false; E = std::max(E, efloor); - const Real fmag = CouplingFluxMagnitude(F); - const Real realizability_slop = - 256.0 * std::numeric_limits::epsilon() * std::max(1.0, E); - if (!std::isfinite(fmag) || fmag > E + realizability_slop) return false; + const Real fmag = VNorm(F); + const Real realizability_slop = 256.0 * Eps() * std::max(1.0, E); + if (fmag > E + realizability_slop) return false; const Real g2 = 1.0 / (1.0 - beta2); const Real g = std::sqrt(g2); @@ -192,20 +181,17 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualBeta( Fr0[1] + d1 * beta[1] + d2 * bdp[1], Fr0[2] + d1 * beta[2] + d2 * bdp[2]}; const auto Fsolve = SolveRadFlux(1.0 + a, bcoef, beta, rhs); - if (!std::isfinite(Fsolve[0]) || !std::isfinite(Fsolve[1]) || !std::isfinite(Fsolve[2])) - return false; const auto Ftarget = ProjectCouplingFluxToEnergy(Fsolve, E); const Real scale = std::max(energy_floor_scale, - std::max(E, std::max(CouplingFluxMagnitude(Fr0), - std::max(fmag, CouplingFluxMagnitude(Ftarget))))); + std::max(E, std::max(VNorm(Fr0), std::max(fmag, VNorm(Ftarget))))); residual_norm = 0.0; for (int d = 0; d < 3; ++d) { residual[d] = F[d] - Ftarget[d]; residual_norm = std::max(residual_norm, std::abs(residual[d]) / scale); } - return std::isfinite(residual_norm); + return true; } template @@ -216,24 +202,21 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( const Real energy_floor_scale, const Real tolerance, std::array &beta, std::array &F, Real &E, Real &dEk, Real &momentum_error, int &iterations) { const Real mu = eref / (dens * c * chat); - if (!std::isfinite(mu) || mu <= 0.0 || !std::isfinite(sigp) || !std::isfinite(sigs) || - sigp < 0.0 || sigs < 0.0) - return false; - const Real solve_tol = std::max(tolerance, 64.0 * std::numeric_limits::epsilon()); + const Real solve_tol = std::max(tolerance, 64.0 * Eps()); constexpr int max_iterations = 32; constexpr int max_line_search = 20; - const Real fd_factor = std::pow(std::numeric_limits::epsilon(), 1.0 / 3.0); + const Real fd_factor = std::pow(Eps(), 1.0 / 3.0); std::array residual{0.0, 0.0, 0.0}; std::array best_beta = beta; std::array best_F = F; Real best_E = E; Real best_dEk = dEk; - Real best_error = std::numeric_limits::max(); + Real best_error = Big(); for (iterations = 1; iterations <= max_iterations; ++iterations) { - Real current_error = std::numeric_limits::max(); + Real current_error = Big(); if (!EvaluateCouplingMomentumEnergyResidualBeta( beta, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, energy_floor_scale, F, E, dEk, residual, current_error)) @@ -254,8 +237,8 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( Real J[3][3]{{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; bool jacobian_valid = true; for (int col = 0; col < 3; ++col) { - const Real h = std::max(32.0 * std::numeric_limits::epsilon(), - fd_factor * std::max(1.0e-3, std::abs(beta[col]))); + const Real h = + std::max(32.0 * Eps(), fd_factor * std::max(1.0e-3, std::abs(beta[col]))); auto beta_plus = beta; auto beta_minus = beta; beta_plus[col] += h; @@ -308,7 +291,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( } const Real step_mag = std::sqrt(SQR(step[0]) + SQR(step[1]) + SQR(step[2])); - if (!std::isfinite(step_mag) || step_mag <= Fuzz()) break; + if (step_mag <= Fuzz()) break; if (step_mag > 0.25) { const Real scale = 0.25 / step_mag; for (int d = 0; d < 3; ++d) @@ -323,7 +306,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( std::array Ftrial, Rtrial; Real Etrial = E; Real dEktrial = dEk; - Real trial_error = std::numeric_limits::max(); + Real trial_error = Big(); const bool trial_valid = EvaluateCouplingMomentumEnergyResidualBeta( beta_trial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, energy_floor_scale, Ftrial, Etrial, dEktrial, Rtrial, trial_error); @@ -352,7 +335,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyBeta( std::array Ftrial, Rtrial; Real Etrial = E; Real dEktrial = dEk; - Real trial_error = std::numeric_limits::max(); + Real trial_error = Big(); const bool trial_valid = EvaluateCouplingMomentumEnergyResidualBeta( beta_trial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, energy_floor_scale, Ftrial, Etrial, dEktrial, Rtrial, trial_error); @@ -392,8 +375,6 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualFlux( std::array &beta, Real &E, Real &dEk, std::array &residual, Real &residual_norm) { const Real mu = eref / (dens * c * chat); - if (!std::isfinite(mu) || mu <= 0.0) return false; - std::array delta_F{0.0, 0.0, 0.0}; Real delta_F2 = 0.0; Real beta0_dot_delta_F = 0.0; @@ -405,21 +386,20 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualFlux( } const Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); - if (!std::isfinite(beta2) || beta2 >= 1.0) return false; + if (beta2 >= 1.0) return false; // This form remains accurate even when beta-beta0 is below the resolution // of a stored velocity component. dEk = c / chat * (-beta0_dot_delta_F + 0.5 * mu * delta_F2); E = E0 + chat / c * (Q - dEg - dEk); - const Real floor_slop = 128.0 * std::numeric_limits::epsilon() * - std::max(1.0, std::max(std::abs(E0), std::abs(E))); - if (!std::isfinite(E) || E < efloor - floor_slop || E <= 0.0) return false; + const Real floor_slop = + 128.0 * Eps() * std::max(1.0, std::max(std::abs(E0), std::abs(E))); + if (E < efloor - floor_slop || E <= 0.0) return false; E = std::max(E, efloor); - const Real fmag = CouplingFluxMagnitude(F); - const Real realizability_slop = - 256.0 * std::numeric_limits::epsilon() * std::max(1.0, E); - if (!std::isfinite(fmag) || fmag > E + realizability_slop) return false; + const Real fmag = VNorm(F); + const Real realizability_slop = 256.0 * Eps() * std::max(1.0, E); + if (fmag > E + realizability_slop) return false; const Real g2 = 1.0 / (1.0 - beta2); const Real g = std::sqrt(g2); @@ -442,20 +422,17 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingMomentumEnergyResidualFlux( Fr0[1] + d1 * beta[1] + d2 * bdp[1], Fr0[2] + d1 * beta[2] + d2 * bdp[2]}; const auto Fsolve = SolveRadFlux(1.0 + a, bcoef, beta, rhs); - if (!std::isfinite(Fsolve[0]) || !std::isfinite(Fsolve[1]) || !std::isfinite(Fsolve[2])) - return false; const auto Ftarget = ProjectCouplingFluxToEnergy(Fsolve, E); const Real scale = std::max(energy_floor_scale, - std::max(E, std::max(CouplingFluxMagnitude(Fr0), - std::max(fmag, CouplingFluxMagnitude(Ftarget))))); + std::max(E, std::max(VNorm(Fr0), std::max(fmag, VNorm(Ftarget))))); residual_norm = 0.0; for (int d = 0; d < 3; ++d) { residual[d] = F[d] - Ftarget[d]; residual_norm = std::max(residual_norm, std::abs(residual[d]) / scale); } - return std::isfinite(residual_norm); + return true; } template @@ -466,24 +443,20 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( const Real energy_floor_scale, const Real tolerance, std::array &beta, std::array &F, Real &E, Real &dEk, Real &momentum_error, int &iterations) { const Real mu = eref / (dens * c * chat); - if (!std::isfinite(mu) || mu <= 0.0 || !std::isfinite(sigp) || !std::isfinite(sigs) || - sigp < 0.0 || sigs < 0.0) - return false; - - const Real solve_tol = std::max(tolerance, 64.0 * std::numeric_limits::epsilon()); + const Real solve_tol = std::max(tolerance, 64.0 * Eps()); constexpr int max_iterations = 32; constexpr int max_line_search = 20; - const Real fd_factor = std::pow(std::numeric_limits::epsilon(), 1.0 / 3.0); + const Real fd_factor = std::pow(Eps(), 1.0 / 3.0); std::array residual{0.0, 0.0, 0.0}; std::array best_beta = beta; std::array best_F = F; Real best_E = E; Real best_dEk = dEk; - Real best_error = std::numeric_limits::max(); + Real best_error = Big(); for (iterations = 1; iterations <= max_iterations; ++iterations) { - Real current_error = std::numeric_limits::max(); + Real current_error = Big(); if (!EvaluateCouplingMomentumEnergyResidualFlux( F, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, energy_floor_scale, beta, E, dEk, residual, current_error)) @@ -510,9 +483,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( // F is normalized by the local energy reference and can be much smaller // than unity. An absolute O(epsilon) perturbation would then be larger // than the correction implied by the requested relative residual. - const Real h = std::max(64.0 * std::numeric_limits::epsilon() * - variable_scale, - fd_factor * variable_scale); + const Real h = std::max(64.0 * Eps() * variable_scale, fd_factor * variable_scale); auto Fplus = F; auto Fminus = F; Fplus[col] += h; @@ -559,16 +530,12 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( } const Real step_mag = std::sqrt(SQR(step[0]) + SQR(step[1]) + SQR(step[2])); - const Real state_scale = std::max( - energy_floor_scale, - std::max(E, std::max(CouplingFluxMagnitude(F), CouplingFluxMagnitude(Fr0)))); + const Real state_scale = + std::max(energy_floor_scale, std::max(E, std::max(VNorm(F), VNorm(Fr0)))); // state_scale includes the radiation floor. Do not impose a unit-scale // absolute cutoff here: for E << 1 it can reject a resolvable correction // before the normalized momentum residual reaches its tolerance. - if (!std::isfinite(step_mag) || step_mag <= 64.0 * - std::numeric_limits::epsilon() * - state_scale) - break; + if (step_mag <= 64.0 * Eps() * state_scale) break; if (step_mag > 0.25 * state_scale) { const Real scale = 0.25 * state_scale / step_mag; for (int d = 0; d < 3; ++d) @@ -583,7 +550,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( std::array beta_trial, Rtrial; Real Etrial = E; Real dEktrial = dEk; - Real trial_error = std::numeric_limits::max(); + Real trial_error = Big(); const bool trial_valid = EvaluateCouplingMomentumEnergyResidualFlux( Ftrial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, energy_floor_scale, beta_trial, Etrial, dEktrial, Rtrial, trial_error); @@ -609,7 +576,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergyFlux( std::array beta_trial, Rtrial; Real Etrial = E; Real dEktrial = dEk; - Real trial_error = std::numeric_limits::max(); + Real trial_error = Big(); const bool trial_valid = EvaluateCouplingMomentumEnergyResidualFlux( Ftrial, beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, energy_floor_scale, beta_trial, Etrial, dEktrial, Rtrial, trial_error); @@ -642,11 +609,11 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingMomentumEnergy( const Real eref, const Real c, const Real chat, const Real sigp, const Real sigs, const Real energy_floor_scale, const Real tolerance, std::array &beta, std::array &F, Real &E, Real &dEk, Real &momentum_error, int &iterations) { - const Real mu = eref / (dens * c * chat); // Choose the nonlinear variable from the relative pseudo-inertia. For // mu <= 1 the gas velocity change can be much less resolvable than the // radiation-flux change, so solve in F. For mu > 1 solve in beta to avoid // magnifying flux perturbations into large velocity perturbations. + const Real mu = eref / (dens * c * chat); if (mu <= 1.0) { return SolveCouplingMomentumEnergyFlux( beta0, Fr0, E0, efloor, B, dEg, Q, dens, eref, c, chat, sigp, sigs, @@ -668,9 +635,9 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingInnerScalarResidual( // Exact reduced-speed-of-light total-energy conservation at fixed velocity: // dEg + dEk + (c/chat) dEr = Q. E = E0 + chat / c * (Q - dEk - (eg - eg0)); - const Real floor_slop = 64.0 * std::numeric_limits::epsilon() * - std::max(1.0, std::max(std::abs(E0), std::abs(E))); - if (!std::isfinite(E) || E < efloor - floor_slop) return false; + const Real floor_slop = + 64.0 * Eps() * std::max(1.0, std::max(std::abs(E0), std::abs(E))); + if (E < efloor - floor_slop) return false; E = std::max(E, efloor); const Real eint = eg * eref / dens; @@ -679,27 +646,18 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingInnerScalarResidual( const Real sigp = chat * dt * opacity.PlanckMeanAbsorptionCoefficient(dens, T); const Real sigs = chat * dt * scattering.RosselandMeanTotalScatteringCoefficient(dens, T); - if (!std::isfinite(T) || !std::isfinite(B) || !std::isfinite(sigp) || - !std::isfinite(sigs) || sigp < 0.0 || sigs < 0.0) - return false; - Real ca = 0.0; Real cb = 0.0; Real cd = 0.0; Real Eeq = E; - if (!ComputeCouplingEnergyCoefficients(sigp, sigs, g, g2, beta2, bdbdp, bdf, ca, cb, - cd) || - !ComputeCouplingEquilibriumEnergy(E0, B, ca, cb, cd, Eeq)) - return false; + ComputeCouplingEnergyCoefficients(sigp, sigs, g, g2, beta2, bdbdp, bdf, ca, cb, cd); + ComputeCouplingEquilibriumEnergy(E0, B, ca, cb, cd, Eeq); // Compare two independently constructed radiation energies: // E: exact total-energy conservation, // Eeq: the implicit radiation source equation. - // Do not evaluate (1+ca)E - (E0+cbB-cd) here. In an optically thick - // near-LTE cell that expression subtracts O(sigp) numbers and its absolute - // roundoff can be orders of magnitude larger than the requested tolerance. residual = E - Eeq; - return std::isfinite(residual); + return true; } KOKKOS_INLINE_FUNCTION Real CouplingInnerScalarError(const Real eg, const Real eg0, @@ -732,8 +690,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( const Real eg_floor = dens * eos.InternalEnergyFromDensityTemperature(dens, tfloor) / eref; const Real eg_ceiling = eg0 + Q - dEk + c / chat * (E0 - efloor); - if (!std::isfinite(eg_floor) || !std::isfinite(eg_ceiling) || eg_ceiling < eg_floor) - return false; + if (eg_ceiling < eg_floor) return false; const Real Tguess = std::max(tfloor, std::pow(std::max(eref * B / arad, 0.0), 0.25)); Real eg_guess = dens * eos.InternalEnergyFromDensityTemperature(dens, Tguess) / eref; @@ -752,7 +709,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( Real best_error = guess_valid ? CouplingInnerScalarError(eg_guess, eg0, dEk, Q, Eguess, E0, chat, c, energy_floor_scale, Rguess) - : std::numeric_limits::max(); + : Big(); if (guess_valid && best_error <= tolerance) { E = Eguess; B = Bguess; @@ -891,7 +848,7 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( // cases. Real eg_trial = 0.5 * (bracket_lo + bracket_hi); const Real denom = bracket_Rhi - bracket_Rlo; - if (std::isfinite(denom) && std::abs(denom) > Fuzz()) { + if (std::abs(denom) > Fuzz()) { const Real eg_secant = (bracket_lo * bracket_Rhi - bracket_hi * bracket_Rlo) / denom; const Real width = bracket_hi - bracket_lo; @@ -927,19 +884,15 @@ KOKKOS_INLINE_FUNCTION bool SolveCouplingInnerScalar( const Real width_scale = std::max(energy_floor_scale, std::abs(bracket_lo) + std::abs(bracket_hi)); - if ((bracket_hi - bracket_lo) / width_scale <= - 8.0 * std::numeric_limits::epsilon()) - break; + if ((bracket_hi - bracket_lo) / width_scale <= 8.0 * Eps()) break; } } E = best_E; B = best_B; inner_err = best_error; - const Real attainable_tolerance = - std::max(tolerance, 32.0 * std::numeric_limits::epsilon()); - return std::isfinite(best_eg) && std::isfinite(E) && std::isfinite(B) && - std::isfinite(inner_err) && inner_err <= attainable_tolerance; + const Real attainable_tolerance = std::max(tolerance, 32.0 * Eps()); + return inner_err <= attainable_tolerance; } template @@ -949,14 +902,12 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingEnergyOnlyResidual( const Real chat, const Real c, const EOSType &eos, const OpacityType &opacity, Real &E, Real &B, Real &residual) { E = E0 + chat / c * (Q - (eg - eg0)); - if (!std::isfinite(E) || E < Emin) return false; + if (E < Emin) return false; const Real eint = eg * eref / dens; const Real T = std::max(tfloor, eos.TemperatureFromDensityInternalEnergy(dens, eint)); B = arad * SQR(SQR(T)) / eref; const Real sigp = chat * dt * opacity.PlanckMeanAbsorptionCoefficient(dens, T); - if (!std::isfinite(T) || !std::isfinite(B) || !std::isfinite(sigp) || sigp < 0.0) - return false; // This fallback intentionally retains only thermal absorption/emission. // Scattering energy exchange is work associated with the momentum source, @@ -966,7 +917,7 @@ KOKKOS_INLINE_FUNCTION bool EvaluateCouplingEnergyOnlyResidual( const Real inv_denom = 1.0 / (1.0 + sigp); const Real Eeq = E0 * inv_denom + (sigp * inv_denom) * B; residual = E - Eeq; - return std::isfinite(Eeq) && std::isfinite(residual); + return true; } //---------------------------------------------------------------------------------------- @@ -999,7 +950,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const auto inner_max = moments_pkg->template Param("inner_iteration_max"); const auto outer_tol = moments_pkg->template Param("outer_iteration_tol"); const auto inner_tol = moments_pkg->template Param("inner_iteration_tol"); - const Real nonlinear_roundoff_tol = 64.0 * std::numeric_limits::epsilon(); + const Real nonlinear_roundoff_tol = 64.0 * Eps(); const auto fatal_if_unconverged = moments_pkg->template Param("fatal_if_unconverged"); @@ -1034,6 +985,31 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const auto jb = u0->GetBoundsJ(IndexDomain::interior); const auto kb = u0->GetBoundsK(IndexDomain::interior); + /* + * This is a multi-step solve with some fallbacks. Energies and fluxes are normalized + * to keep numbers around unity. + * + * If the cells are close to the density floor, we skip the full momentum solve and + * only do energy coupling. This prevents large accelerations. + * + * For the outer loop we: + * - Hold the current gas velocity and radiation flux estimate fixed. + * - It chooses flux or velocity (beta) as the Newton variable based on mu = eref / + * (rho * c * chat) + * - Run an inner nonlinear solve for material emission B and radiation energy E. + * - It uses a damped quasi-Newton update with line search. + * - If that stalls, it switches to a bracketed scalar solve in gas internal + * energy. + * - If needed, it can use a momentum-energy predictor to seed the next outer + * iteration. + * - With the converged thermal state, solve the coupled momentum system: + * + * If no full iterate is available, it runs the conservative energy-only fallback: + * - freeze gas momentum and radiation flux; + * - bracket a thermal solution in gas internal energy; + * - use it automatically for atmosphere cells and floor-constrained failures. + * + */ parthenon::par_for( DEFAULT_LOOP_PATTERN, "MatterCoupling", DevExecSpace(), 0, u0->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, @@ -1044,17 +1020,16 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // U^(0) values const Real dens_raw = v0(b, gas::cons::density(), k, j, i); - const Real dens = (std::isfinite(dens_raw) && dens_raw > dflr) ? dens_raw : dflr; + const Real dens = (dens_raw > dflr) ? dens_raw : dflr; Real Q = 0.0; - if (do_raytrace) Q = FiniteOrZero(dt * v0(b, gas::src::energy(), k, j, i)); + if (do_raytrace) Q = dt * v0(b, gas::src::energy(), k, j, i); // In the numerical atmosphere, retain thermal emission/absorption but // suppress the momentum update. This avoids accelerating floor-density // gas while still allowing absorbed raytraced energy to reradiate into // the moment field through the conservative energy-only fallback. constexpr Real atmosphere_floor_factor = 100.0; - const bool numerical_atmosphere = - !std::isfinite(dens_raw) || dens_raw <= atmosphere_floor_factor * dflr; + const bool numerical_atmosphere = dens_raw <= atmosphere_floor_factor * dflr; // Note(AMD): There is some floating point difference between the // internal energy used to compute the temperature and the internal @@ -1086,24 +1061,14 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // have E_r << a_r T_floor^4. Promoting E_r to Bfloor_phys here creates // an artificial LTE radiation bath and can make the constrained // matter-coupling equations inconsistent at the first timestep. - const Real Er_state = FiniteOrZero(v0(b, rad::cons::energy(), k, j, i)); + const Real Er_state = v0(b, rad::cons::energy(), k, j, i); Real E0 = std::max(rad_efloor, Er_state); // Use the arithmetic mean as the reference scale. Unlike the // geometric mean, this keeps both normalized energies bounded when the // gas and radiation temperatures are initially very different. const Real eref_max = std::max(E0, B); const Real eref_min = std::min(E0, B); - Real eref = 0.5 * eref_max * (1.0 + eref_min / std::max(eref_max, Fuzz())); - if (!std::isfinite(eref) || eref <= 0.0) eref = eref_max; - if (!std::isfinite(eref) || eref <= 0.0) { - v0(b, gas::cons::internal_energy(), k, j, i) += Q; - v0(b, gas::cons::total_energy(), k, j, i) += Q; - v0(b, rad::cons::energy(), k, j, i) = rad_efloor; - v0(b, rad::cons::flux(0), k, j, i) = 0.0; - v0(b, rad::cons::flux(1), k, j, i) = 0.0; - v0(b, rad::cons::flux(2), k, j, i) = 0.0; - return; - } + Real eref = 0.5 * eref_max * (1.0 + eref_min / std::max(eref_max, Fuzz())); const Real fref = c * eref; // Keep the radiation floor independent of the material temperature // floor. Bfloor constrains T; efloor constrains E_r. @@ -1115,16 +1080,16 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { eg0 /= eref; B /= eref; - const auto fred0 = NormalizeFlux( - FiniteOrZero(v0(b, rad::cons::flux(0), k, j, i) / (hx[0] * fref * E0)), - FiniteOrZero(v0(b, rad::cons::flux(1), k, j, i) / (hx[1] * fref * E0)), - FiniteOrZero(v0(b, rad::cons::flux(2), k, j, i) / (hx[2] * fref * E0))); + const auto fred0 = + NormalizeFlux(v0(b, rad::cons::flux(0), k, j, i) / (hx[0] * fref * E0), + v0(b, rad::cons::flux(1), k, j, i) / (hx[1] * fref * E0), + v0(b, rad::cons::flux(2), k, j, i) / (hx[2] * fref * E0)); const std::array Fr0{E0 * fred0[0], E0 * fred0[1], E0 * fred0[2]}; std::array v{p0[0] / dens, p0[1] / dens, p0[2] / dens}; const std::array beta0{v[0] / c, v[1] / c, v[2] / c}; const Real beta20 = SQR(beta0[0]) + SQR(beta0[1]) + SQR(beta0[2]); - if (!std::isfinite(beta20) || beta20 >= 1.0) { + if (beta20 >= 1.0) { v0(b, gas::cons::internal_energy(), k, j, i) += Q * eref; v0(b, gas::cons::total_energy(), k, j, i) += Q * eref; v0(b, rad::cons::energy(), k, j, i) = E0 * eref; @@ -1141,22 +1106,21 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // Start outer iteration int outer_iter = 0; int inner_iter = 0; - Real outer_err = std::numeric_limits::max(); - Real inner_err = std::numeric_limits::max(); + Real outer_err = Big(); + Real inner_err = Big(); std::array dv{0.0, 0.0, 0.0}; Real dEk = 0.0; Real dEg = 0.0; Real dEr = 0.0; - const Real energy_floor_scale = std::max(efloor + Bfloor, Fuzz()); + const Real energy_floor_scale = std::max(efloor + Bfloor, Fuzz()); Real escale = std::max( energy_floor_scale, std::max(std::abs(eg0), std::max(std::abs(Q), c / chat * std::abs(E0)))); - bool solve_valid = !numerical_atmosphere && std::isfinite(escale) && - std::isfinite(ke0) && escale > 0.0; + bool solve_valid = !numerical_atmosphere && escale > 0.0; bool outer_converged = false; bool have_complete_iterate = false; - Real best_outer_err = std::numeric_limits::max(); + Real best_outer_err = Big(); std::array best_F = Fr0; std::array best_dv{0.0, 0.0, 0.0}; Real best_dEg = 0.0; @@ -1165,7 +1129,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { bool scalar_inner_used = false; int momentum_predictor_count = 0; int momentum_iteration_count = 0; - Real momentum_error = std::numeric_limits::max(); + Real momentum_error = Big(); for (outer_iter = 1; outer_iter <= outer_max; ++outer_iter) { if (!solve_valid) break; @@ -1173,7 +1137,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const Real ke = 0.5 * dens * (SQR(v[0]) + SQR(v[1]) + SQR(v[2])) / eref; std::array beta{v[0] / c, v[1] / c, v[2] / c}; const Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); - if (!std::isfinite(beta2) || beta2 >= 1.0 || E <= 0.0) { + if (beta2 >= 1.0 || E <= 0.0) { solve_valid = false; break; } @@ -1192,9 +1156,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const Real bdbdp = beta[0] * bdp[0] + beta[1] * bdp[1] + beta[2] * bdp[2]; const Real bdf = beta[0] * F[0] + beta[1] * F[1] + beta[2] * F[2]; - // Damped quasi-Newton solve for (B,E). Opacity derivatives are not - // available, so accept only residual-decreasing trial steps. - Real previous_inner_err = std::numeric_limits::max(); + // Damped quasi-Newton solve for (B,E). + Real previous_inner_err = Big(); int stalled_iterations = 0; bool inner_converged = false; for (inner_iter = 1; inner_iter <= inner_max; ++inner_iter) { @@ -1210,12 +1173,10 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { Real ca = 0.0; Real cb = 0.0; Real cd = 0.0; - const bool coeff_valid = ComputeCouplingEnergyCoefficients( - sigp, sigs, g, g2, beta2, bdbdp, bdf, ca, cb, cd); + ComputeCouplingEnergyCoefficients(sigp, sigs, g, g2, beta2, bdbdp, bdf, ca, + cb, cd); Real Eeq_inner = E; - const bool equilibrium_valid = - coeff_valid && - ComputeCouplingEquilibriumEnergy(E0, B, ca, cb, cd, Eeq_inner); + ComputeCouplingEquilibriumEnergy(E0, B, ca, cb, cd, Eeq_inner); // Retain the original residuals only to form the quasi-Newton // search direction. Use the conservation/equilibrium form below @@ -1235,9 +1196,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { std::abs(ke - ke0))))); inner_err = std::max(std::abs(conservation_inner) / escale_inner, c / chat * std::abs(source_inner) / escale_inner); - solve_valid = solve_valid && equilibrium_valid && std::isfinite(inner_err) && - std::isfinite(fleck) && std::isfinite(eint) && Cv > 0.0 && - fleck >= 0.0 && sigp >= 0.0 && sigs >= 0.0; + solve_valid = + solve_valid && Cv > 0.0 && fleck >= 0.0 && sigp >= 0.0 && sigs >= 0.0; if (!solve_valid) break; if (inner_err <= inner_tol) { inner_converged = true; @@ -1253,7 +1213,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const Real dfac = 1.0 + c / chat * fleck * cb; const Real denom = dfac + ca; - if (!std::isfinite(denom) || std::abs(denom) <= Fuzz()) { + if (std::abs(denom) <= Fuzz()) { solve_valid = false; break; } @@ -1281,14 +1241,11 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { Real ca_trial = 0.0; Real cb_trial = 0.0; Real cd_trial = 0.0; - const bool coeff_trial_valid = ComputeCouplingEnergyCoefficients( - sigp_trial, sigs_trial, g, g2, beta2, bdbdp, bdf, ca_trial, cb_trial, - cd_trial); + ComputeCouplingEnergyCoefficients(sigp_trial, sigs_trial, g, g2, beta2, + bdbdp, bdf, ca_trial, cb_trial, cd_trial); Real Eeq_trial = Etrial; - const bool equilibrium_trial_valid = - coeff_trial_valid && - ComputeCouplingEquilibriumEnergy(E0, Btrial, ca_trial, cb_trial, - cd_trial, Eeq_trial); + ComputeCouplingEquilibriumEnergy(E0, Btrial, ca_trial, cb_trial, cd_trial, + Eeq_trial); const Real conservation_trial = (ke - ke0) + (eint_trial - eg0) + c / chat * (Etrial - E0) - Q; const Real source_trial = Etrial - Eeq_trial; @@ -1301,18 +1258,13 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { const Real trial_err = std::max(std::abs(conservation_trial) / escale_trial, c / chat * std::abs(source_trial) / escale_trial); - const bool trial_valid = - std::isfinite(Ttrial) && std::isfinite(eint_trial) && - std::isfinite(sigp_trial) && std::isfinite(sigs_trial) && - equilibrium_trial_valid && std::isfinite(trial_err) && - sigp_trial >= 0.0 && sigs_trial >= 0.0; - if (trial_valid && trial_err < best_trial_err) { + if (trial_err < best_trial_err) { best_trial_err = trial_err; best_trial_E = Etrial; best_trial_B = Btrial; } - if (trial_valid && (trial_err <= inner_tol || - trial_err <= inner_err * (1.0 - 1.0e-4 * alpha))) { + if ((trial_err <= inner_tol || + trial_err <= inner_err * (1.0 - 1.0e-4 * alpha))) { E = Etrial; B = Btrial; accepted = true; @@ -1372,7 +1324,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { std::array Fpred = F; Real Epred = E; Real dEkpred = dEk; - Real predictor_error = std::numeric_limits::max(); + Real predictor_error = Big(); int predictor_iterations = 0; const bool predictor_converged = SolveCouplingMomentumEnergy( beta0, Fr0, E0, efloor, B, dEg_pred, Q, dens, eref, c, chat, sigp_pred, @@ -1381,10 +1333,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { dEkpred, predictor_error, predictor_iterations); momentum_iteration_count += predictor_iterations; if (predictor_converged) { - const Real predictor_scale = - std::max(energy_floor_scale, - std::max(E, std::max(CouplingFluxMagnitude(F), - CouplingFluxMagnitude(Fpred)))); + const Real predictor_scale = std::max( + energy_floor_scale, std::max(E, std::max(VNorm(F), VNorm(Fpred)))); Real predictor_change = std::abs(Epred - E) / predictor_scale; for (int d = 0; d < 3; ++d) { predictor_change = std::max(predictor_change, @@ -1392,7 +1342,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { predictor_change = std::max(predictor_change, std::abs(beta_pred[d] - v[d] / c)); } - if (predictor_change > 8.0 * std::numeric_limits::epsilon()) { + if (predictor_change > 8.0 * Eps()) { E = Epred; F = Fpred; dEk = dEkpred; @@ -1418,11 +1368,6 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { chat * dt * opac_d.RosselandMeanAbsorptionCoefficient(dens, T); const Real sigs_flux = chat * dt * scat_d.RosselandMeanTotalScatteringCoefficient(dens, T); - if (!std::isfinite(sigp_flux) || !std::isfinite(sigs_flux) || sigp_flux < 0.0 || - sigs_flux < 0.0) { - solve_valid = false; - break; - } // Solve the radiation-flux source equation together with gas // pseudo-momentum conservation and reduced-c total-energy @@ -1453,7 +1398,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // the kinetic-energy correction changed little in a Keplerian flow. const std::array beta_out{v[0] / c, v[1] / c, v[2] / c}; const Real beta2_out = SQR(beta_out[0]) + SQR(beta_out[1]) + SQR(beta_out[2]); - if (!std::isfinite(beta2_out) || beta2_out >= 1.0 || E <= 0.0) { + if (beta2_out >= 1.0 || E <= 0.0) { solve_valid = false; break; } @@ -1486,12 +1431,10 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { Real cb_out = 0.0; Real cd_out = 0.0; Real Eeq_out = E; - const bool energy_coeff_valid = ComputeCouplingEnergyCoefficients( - sigp_energy, sigs_out, g_out, g2_out, beta2_out, bdbdp_out, bdf_out, ca_out, - cb_out, cd_out); - const bool equilibrium_valid = - energy_coeff_valid && - ComputeCouplingEquilibriumEnergy(E0, B, ca_out, cb_out, cd_out, Eeq_out); + ComputeCouplingEnergyCoefficients(sigp_energy, sigs_out, g_out, g2_out, + beta2_out, bdbdp_out, bdf_out, ca_out, cb_out, + cd_out); + ComputeCouplingEquilibriumEnergy(E0, B, ca_out, cb_out, cd_out, Eeq_out); const Real conservation_residual = dEg + dEk + c / chat * dEr - Q; const Real source_residual = E - Eeq_out; @@ -1521,16 +1464,14 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { c / chat * std::abs(source_residual) / escale); const Real flux_scale = std::max( energy_floor_scale, - std::max(E, std::max(CouplingFluxMagnitude(Fr0), - std::max(CouplingFluxMagnitude(F), - CouplingFluxMagnitude(Ftarget_out))))); + std::max(E, std::max(VNorm(Fr0), std::max(VNorm(F), VNorm(Ftarget_out))))); Real flux_residual = 0.0; for (int d = 0; d < 3; ++d) { flux_residual = std::max(flux_residual, std::abs(F[d] - Ftarget_out[d]) / flux_scale); } - const Real realizability_residual = std::max( - 0.0, (CouplingFluxMagnitude(F) - E) / std::max(E, energy_floor_scale)); + const Real realizability_residual = + std::max(0.0, (VNorm(F) - E) / std::max(E, energy_floor_scale)); // These are residuals of the coupled equations evaluated at the // current state. Iterate-to-iterate changes are deliberately not an // acceptance criterion: a converged nonlinear root need not move by @@ -1538,10 +1479,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // resolve. outer_err = std::max(energy_residual, std::max(flux_residual, realizability_residual)); - solve_valid = solve_valid && equilibrium_valid && std::isfinite(outer_err) && - std::isfinite(dEr) && std::isfinite(sigp_energy) && - std::isfinite(sigs_out) && std::isfinite(sigp_flux_out) && - sigp_energy >= 0.0 && sigs_out >= 0.0 && sigp_flux_out >= 0.0; + solve_valid = solve_valid && sigp_energy >= 0.0 && sigs_out >= 0.0 && + sigp_flux_out >= 0.0; if (!solve_valid) break; have_complete_iterate = true; @@ -1566,16 +1505,16 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { // energy-only solve below with momentum frozen, rather than repeatedly // applying a radiation drag whose kinetic-energy work cannot be paid by // either the gas or radiation field. - const Real floor_active_tol = 256.0 * std::numeric_limits::epsilon(); + const Real floor_active_tol = 256.0 * Eps(); const bool material_floor_active = - B <= Bfloor * (1.0 + floor_active_tol) + Fuzz(); + B <= Bfloor * (1.0 + floor_active_tol) + Fuzz(); const bool radiation_floor_active = - E <= efloor * (1.0 + floor_active_tol) + Fuzz(); + E <= efloor * (1.0 + floor_active_tol) + Fuzz(); const bool floor_constrained_failure = solve_valid && !have_complete_iterate && (material_floor_active || radiation_floor_active); bool fallback_success = false; - Real fallback_error = std::numeric_limits::max(); + Real fallback_error = Big(); if (have_complete_iterate) { F = best_F; @@ -1595,13 +1534,12 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { dv = {0.0, 0.0, 0.0}; v = {p0[0] / dens, p0[1] / dens, p0[2] / dens}; dEk = 0.0; - const Real Emin = std::max(efloor, CouplingFluxMagnitude(Fr0)); + const Real Emin = std::max(efloor, VNorm(Fr0)); const Real eg_floor = dens * eos_d.InternalEnergyFromDensityTemperature(dens, tfloor) / eref; const Real eg_max = eg0 + Q + c / chat * (E0 - Emin); - bool fallback_valid = - std::isfinite(eg_floor) && std::isfinite(eg_max) && eg_max >= eg_floor; + bool fallback_valid = eg_max >= eg_floor; bool fallback_have = false; bool bracketed = false; Real bracket_lo = eg_floor; @@ -1610,11 +1548,12 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { Real best_eg = eg0; Real best_E_fallback = E0; Real best_B_fallback = B; - Real best_residual = std::numeric_limits::max(); + Real best_residual = Big(); Real previous_eg = eg_floor; Real previous_residual = 0.0; bool have_previous = false; + // bisection fallback #1 if (fallback_valid) { constexpr int fallback_scan_points = 64; for (int n = 0; n <= fallback_scan_points; ++n) { @@ -1689,9 +1628,8 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { energy_floor_scale, std::max(std::abs(Q), std::max(std::abs(eg0) + std::abs(best_eg), c / chat * (std::abs(E0) + std::abs(E))))); - fallback_error = - c / chat * best_residual / std::max(fallback_scale, Fuzz()); - fallback_success = std::isfinite(fallback_error); + fallback_error = c / chat * best_residual / std::max(fallback_scale, Fuzz()); + fallback_success = true; } else { // Last-resort finite update. This is intentionally limited to the // external ray source and leaves the moment state unchanged. @@ -1707,12 +1645,12 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { floor_constrained_failure && fallback_success; if (!outer_converged && fatal_if_unconverged && !numerical_atmosphere && !accepted_floor_fallback) { +#ifndef NDEBUG + // printf with args is slow on device, so hide this at compile time const Real beta_fail = std::sqrt(SQR(v[0] / c) + SQR(v[1] / c) + SQR(v[2] / c)); - const Real fred_fail = - CouplingFluxMagnitude(F) / std::max(E, energy_floor_scale); + const Real fred_fail = VNorm(F) / std::max(E, energy_floor_scale); const Real mu_fail = eref / (dens * c * chat); - const Real delta_f_fail = - CouplingFluxMagnitude({F[0] - Fr0[0], F[1] - Fr0[1], F[2] - Fr0[2]}); + const Real delta_f_fail = VNorm({F[0] - Fr0[0], F[1] - Fr0[1], F[2] - Fr0[2]}); printf("MatterCoupling full fail (%d,%d,%d,%d): outer=%.17e " "(best=%.17e, tol=%.17e), inner=%.17e (tol=%.17e), " "fallback=%.17e, E=%.17e, E0=%.17e, efloor=%.17e, " @@ -1730,6 +1668,7 @@ TaskStatus MatterCouplingFullSingleImpl(MeshData *u0, const Real dt) { momentum_iteration_count, momentum_error, solve_valid, scalar_inner_used, momentum_predictor_count, have_complete_iterate, floor_constrained_failure); +#endif PARTHENON_FAIL("Outer not converged"); } From 3484005046658a8d6478746141701366e240bcd0 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Tue, 14 Jul 2026 08:48:23 -0600 Subject: [PATCH 10/14] Update chicoma detection --- env/bash | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/env/bash b/env/bash index ab3ee9b7..54fc07a9 100644 --- a/env/bash +++ b/env/bash @@ -21,7 +21,7 @@ export MAKE_PROGRAM=${MAKE_PROGRAM:-make} # Identify partitions based on SLURM_JOB_PARTITION and HOSTNAME variables PARTITION="unknown" -if [[ $HOSTNAME == ch-fe* || $CLUSTER == "chicoma" ]]; then +if [[ $HOSTNAME == ch* || $CLUSTER == "chicoma" ]]; then if [[ $SLURM_JOB_PARTITION == *gpu* ]]; then PARTITION="chicoma-gpu" else From cdb6b4114099dc77efabf3585791009f8529af42 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Tue, 14 Jul 2026 10:07:40 -0600 Subject: [PATCH 11/14] Remove "simple" matter coupling --- src/radiation/moments/moments.cpp | 20 ++------------------ src/radiation/params.yaml | 4 ---- 2 files changed, 2 insertions(+), 22 deletions(-) diff --git a/src/radiation/moments/moments.cpp b/src/radiation/moments/moments.cpp index e8f99b91..9c8f4b8c 100644 --- a/src/radiation/moments/moments.cpp +++ b/src/radiation/moments/moments.cpp @@ -18,7 +18,6 @@ #include "artemis.hpp" #include "geometry/geometry.hpp" #include "matter_coupling.hpp" -#include "matter_coupling_simple.hpp" #include "moments.hpp" #include "utils/artemis_utils.hpp" #include "utils/fluxes/fluid_fluxes.hpp" @@ -87,12 +86,6 @@ std::shared_ptr Initialize(ParameterInput *pin, params.Add("fatal_if_unconverged", pin->GetOrAddBoolean("radiation/moment", "fatal_if_unconverged", true)); - // how to handle the matter coupling: - // full_coupling = false only does a loop over energy coupling - // full_coupling = true also does an outer loop over momentum coupling - params.Add("full_coupling", - pin->GetOrAddBoolean("radiation/moment", "full_coupling", true)); - // Radiation constants (including chat for Moments) // NOTE(@pdmullen): These are also stored in top level radiation package... const Real light = constants.GetCCode(); @@ -439,21 +432,12 @@ TaskStatus MatterCoupling(MeshData *u0, const Real dt) { // Extract moments package and params auto &moments_pkg = pm->packages.Get("moments"); auto closure_type = moments_pkg->template Param("closure_type"); - auto full_coupling = moments_pkg->template Param("full_coupling"); // Call MatterCoupling with appropriate GEOM, Fluid, and Closure type given coupling if (closure_type == Closure::m1) { - if (full_coupling) { - return MatterCouplingFullSingleImpl(u0, dt); - } else { - // return MatterCouplingSimpleImpl(u0, dt); - } + return MatterCouplingFullSingleImpl(u0, dt); } else if (closure_type == Closure::p1) { - if (full_coupling) { - return MatterCouplingFullSingleImpl(u0, dt); - } else { - // return MatterCouplingSimpleImpl(u0, dt); - } + return MatterCouplingFullSingleImpl(u0, dt); } return TaskStatus::complete; } diff --git a/src/radiation/params.yaml b/src/radiation/params.yaml index e2a98ddd..0d34b426 100644 --- a/src/radiation/params.yaml +++ b/src/radiation/params.yaml @@ -45,10 +45,6 @@ radiation: _type: Real _description: "CFL number for radiation" _default: 0.8 - full_coupling: - _type: bool - _description: "If set to true, this will do a full outer iteration of momentum exchange with the material. If set to false, artemis will take one outer iteration of momentum coupling." - _default: true creduc: _type: Real _description: "Speed of light reduction factor, c_hat = c / creduc." From 5265146805a5fb076ca433d6bb5e679b96d9e1ed Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Tue, 14 Jul 2026 10:32:35 -0600 Subject: [PATCH 12/14] Add substep to docs --- src/radiation/params.yaml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/radiation/params.yaml b/src/radiation/params.yaml index 0d34b426..1b2dd9b8 100644 --- a/src/radiation/params.yaml +++ b/src/radiation/params.yaml @@ -112,6 +112,10 @@ radiation: _type: bool _description: "When dynamically initializing the radiation field, scale a T^4 by min(1, tau), where tau is the optical depth of the cell" _default: "true" + substep: + _type: bool + _description: "Enable radiation subcycling" + _default: false imc: _description: "IMC radiation model from the Jaybenne library" _type: node From 108a6b963bf0a91b1930a05fb4c04a48217c15c3 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Tue, 14 Jul 2026 10:34:05 -0600 Subject: [PATCH 13/14] No need to check opacities being finite or negative --- src/radiation/raytrace/raytrace.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/radiation/raytrace/raytrace.cpp b/src/radiation/raytrace/raytrace.cpp index 0ee3502b..cb95c180 100644 --- a/src/radiation/raytrace/raytrace.cpp +++ b/src/radiation/raytrace/raytrace.cpp @@ -307,9 +307,7 @@ TaskStatus EvalOpac(MeshData *md) { // Evaluated at T* //%%%%%%%%%%%%%%%% const Real temp = eos_d.TemperatureFromDensityInternalEnergy(rho, sie); - const Real alpha = opacity_d.AbsorptionCoefficient(rho, temp, 1.0); - vmesh(b, rad::star::absorption(), k, j, i) = - (std::isfinite(alpha) && alpha > 0.0) ? alpha : 0.0; + vmesh(b, rad::star::absorption(), k, j, i) = opacity_d.AbsorptionCoefficient(rho, temp, 1.0); }); return TaskStatus::complete; From e26ebb2a14e8102bbfa2293384232299fe6f0d43 Mon Sep 17 00:00:00 2001 From: Adam Dempsey Date: Mon, 20 Jul 2026 07:20:56 -0600 Subject: [PATCH 14/14] format --- src/radiation/raytrace/raytrace.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/radiation/raytrace/raytrace.cpp b/src/radiation/raytrace/raytrace.cpp index cb95c180..04548a0a 100644 --- a/src/radiation/raytrace/raytrace.cpp +++ b/src/radiation/raytrace/raytrace.cpp @@ -307,7 +307,8 @@ TaskStatus EvalOpac(MeshData *md) { // Evaluated at T* //%%%%%%%%%%%%%%%% const Real temp = eos_d.TemperatureFromDensityInternalEnergy(rho, sie); - vmesh(b, rad::star::absorption(), k, j, i) = opacity_d.AbsorptionCoefficient(rho, temp, 1.0); + vmesh(b, rad::star::absorption(), k, j, i) = + opacity_d.AbsorptionCoefficient(rho, temp, 1.0); }); return TaskStatus::complete;