Improve moments radiation - #137
Conversation
| // 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 |
There was a problem hiding this comment.
Does your text editor or claude need to be aware of 90 character limit?
| using ArtemisUtils::VI; | ||
| using ArtemisUtils::VNorm; | ||
|
|
||
| namespace Moments { |
There was a problem hiding this comment.
Just leaving this comment at the top of this...
This solver is now substantially more complex than the previous coupling path. Even if just by private communication (as opposed to an additional CI test), can you show me how this large changeset impacts the solution? Is it avoiding crashes... or providing a higher fidelity solution? Both?
If you really did want to go the route of an added unit test... one where the old solver fails/new solver passes might be worthwhile.
I don’t want to hold this PR up unnecessarily, especially since I understand users need these improvements for active runs. That said, I’m having a hard time understanding which parts of this large change set are actually responsible for the improved behavior.
From the diff, it looks like much of the substantive work is in the matter-coupling solver, so I’m not convinced the improvement is simply from moving moments radiation from an operator-split/subcycled update into the unsplit integrator. Could you help identify the main “difference maker” in this file?
If the robustness improvement is mostly isolated to one part of the rewrite, I’d like us to consider whether the PR can be narrowed to that piece, or at least documented in a way that makes the necessity of the other changes clear. Right now, without piece-by-piece evidence or a focused reproducer showing which changes improve which failures, I’m being asked to take a fairly large leap of faith that the full change set is needed.
| } | ||
| } // outer iteration | ||
|
|
||
| // A hard floor turns the nonlinear equality problem into a |
There was a problem hiding this comment.
My goodness.... this is starting to feel like HARM's fixup.cpp. How often are we hitting this? Everywhere in atmosphere? Did you ever set a boolean to see where these fixups are hit in practice?
Or in a realistic problem how bad conservation gets?
| void InitMesh(parthenon::Mesh *pmesh); | ||
| //---------------------------------------------------------------------------------------- | ||
| //! \fn Real Moments::EstimateTimeStepMesh | ||
| //! \brief Not enrolled in parthenon's determination for global dt |
There was a problem hiding this comment.
huh? this comment isn't true anymore if we are enrolling moments in an unsplit manner, right?
On that note---is taking the radiation dt for hydro + other physics packages even a viable path for us? I would have thought that even with reduced speed of light that this would be a showstopper compared to just subcycling the radiatino physics in an operator split update.
What is the recommendation for users then? Do we want to continue supporting both paths. It does make our driver a bit more unwieldy.
| auto &rt_pkg = pm->packages.Get("raytrace"); | ||
| auto x1max = rt_pkg->Param<Real>("x1max"); | ||
| const Real x1tol = | ||
| 32.0 * std::numeric_limits<Real>::epsilon() * std::max(1.0, std::abs(x1max)); |
| _type: Real | ||
| _description: "CFL number for radiation" | ||
| _default: 0.8 | ||
| full_coupling: |
There was a problem hiding this comment.
does this need to be eliminated from any input files or physics.rst?
| "radiation/integrator must be rk1,rk2, or rk3.") | ||
| rad_integrator = std::make_unique<Integrator_t>(rad_int); | ||
| do_moment_split = pin->GetOrAddBoolean("radiation/moment", "substep", false); | ||
| do_moment_unsplit = !do_moment_split; |
There was a problem hiding this comment.
why do we want rad_int inside pin when substep==false? That is confusing because rad_int != parthenon/time/integrator is not allowed for unsplit physics.
| @@ -83,11 +83,17 @@ ArtemisDriver<GEOM>::ArtemisDriver(ParameterInput *pin, ApplicationInput *app_in | |||
| ndim = pm->ndim; | |||
There was a problem hiding this comment.
Does everything work correctly with package turn-on and turn-off if using moments unsplit mode? In that event, do we still take the radiation dt?
2 participants
Background
This does 2 things:
The cost of this new approach is definitely more, but I think the robustness might make it worth it.
This change may not play nice with restarts.
Description of Changes
Checklist
// This file was created in part or in whole by generative AI