Improvements to matter coupling solver#128
Conversation
| std::max(B, static_cast<Real>(0.0))); | ||
| const Real eref_avg = 0.5 * (std::abs(E0) + std::abs(B)); | ||
| const Real eref = | ||
| std::max(erad_floor, std::max(eref_geom, std::max(eref_avg, Fuzz<Real>()))); |
There was a problem hiding this comment.
Basic math question here - eref_geom would always be less than eref_avg right? Could eref_geom be removed then, and this eref be made a max over erad_floor, eref_avg and Fuzz<Real>()?
There was a problem hiding this comment.
(This is the AM-GM inequality, basically, and the entries for either average are always positive here - in fact the absolute values will be greater in general than the max with 0... sorry if missing something here.)
| // choose a strictly positive reference scale to keep the normalized solve finite | ||
| const Real erad_floor = arad * SQR(SQR(tfloor)); | ||
| const Real eref_geom = std::sqrt(std::max(E0, static_cast<Real>(0.0)) * | ||
| std::max(B, static_cast<Real>(0.0))); |
There was a problem hiding this comment.
Are the static_cast<Real>(0.0) necessary here? (Assuming eref_geom is necessary - see comment on eref_geom.)
| const Real icc = 1. / (c * chat * dens); | ||
| Real escale = et0 + c / chat * E; | ||
| const Real escale_floor = std::max(efloor, Fuzz<Real>()); | ||
| const Real beta2_max = 1.0 - 1.0e-12; |
There was a problem hiding this comment.
| const Real beta2_max = 1.0 - 1.0e-12; | |
| constexpr Real beta2_max = 1.0 - 1.0e-12; |
?
| std::array<Real, 3> beta{v[0] / c, v[1] / c, v[2] / c}; | ||
| const Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); | ||
| Real beta2 = SQR(beta[0]) + SQR(beta[1]) + SQR(beta[2]); | ||
| if (beta2 > beta2_max) { |
There was a problem hiding this comment.
If beta2 is greater than beta2_max, that presumably means this is ultra-relativistic, right? Is it worth adding a PARTHENON_DEBUG check here that beta2 < 1.0?
| flux_change = std::max(flux_change, std::abs(F[d] - F_prev[d]) / flux_scale); | ||
| const Real vel_scale = std::abs(v[d]) + std::abs(v_prev[d]) + Fuzz<Real>(); | ||
| vel_change = std::max(vel_change, std::abs(v[d] - v_prev[d]) / vel_scale); | ||
| } |
There was a problem hiding this comment.
Is it worth moving the computation of flux_scale and vel_scale out of this loop, in to a separate loop over dimension? (Furthermore, is it worth restricting this loop to the number of problem dimensions?) It seems like this could make the subsequent outer_error calculation below too conservative.
|
|
||
| // This needs to be something else related to the change from the last iteration | ||
| outer_err = std::abs(dEk - dEk_prev) / escale; | ||
| auto fedd_outer = |
There was a problem hiding this comment.
Is all this _outer computation necessary (e.g. the Eddington tensor stuff)? It seems to be recomputing the whole set of equations for error terms... it may be worth adding some comments (or having Claude add some comments) on these lines.
2 participants
Background
This explores some AI suggested improvements to the radiation-matter coupling.
Description of Changes
Checklist
// This file was created in part or in whole by generative AI