Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generalize Integrator Infrastructure #840

Merged
merged 23 commits into from
Apr 16, 2023
Merged

Generalize Integrator Infrastructure #840

merged 23 commits into from
Apr 16, 2023

Conversation

Yurlungur
Copy link
Collaborator

@Yurlungur Yurlungur commented Feb 24, 2023

PR Summary

I am pushing on high-order methods on parthenon downstream. To prepare for this, I needed some higher-order integrators. For that reason, I've implemented several improvements to our integrator machinery, without breaking any backwards compatibility:

Class hierarchy

I have split the StagedIntegrator class into a class hierarchy. The base class contains a little shared functionality. There are currently two subclasses:

  • LowStorageIntegrator, which contains the family of low-storage integrators used by Athena++ expressed in Shu Osher form. (This is what @pgrete updated our integrators to some time ago.)
  • ButcherIntegrator which contains a few commonly-used Butcher tableaus plus a 10th order one for fun. The butcher tableau machinery is introduced with IMEX and multi-rate integrators in mind down the line.

I also have templated the MultiStageDriver class on integrator type, with the default being LowStorageIntegrator, which returns to the old machinery. (To prevent people from having to type template<> when inheriting from the MultiStageDriver, I name the generic one MultiStageDriverGeneric and use a type alias to produce MultiStageDriver.)

The intent here is that if one knows they will be only a specific integration scheme, they should inherit from an iterative driver specialized to that specific integrator. However, the class hierarchy opens the possibility of switching between integration schemes at runtime.

In a similar vein, I have added some functions to update.hpp that can take the integrator object and do the update by reading it. I have not updated our examples or any downstream codes to use these functions, so they may have bugs. Nevertheless I'd like to merge them in if possible and fix them later as I work on my high-order machinery. If people object, I can remove them. I have some discussion points about this below.

Additional changes:

  • There was an ambiguity in our old machinery about stages vs. storage buffers. The integrator now reports not only how many stages it needs, but how many storage buffers.
  • I have added a low-storage version of SSPRK4(5) to the LowStorageIntegrator. For parthenon-hydro and AthenaPK this should offer immediate benefits if you want to try it out.
  • I have commented where exactly our equations for integration come from in more detail. I found the Athena++ paper wasn't enough for me to fully understand what was going on, so I dig into the low storage paper that is referenced there.
  • Our integration coefficients were not tested in isolation. I added a set of unit tests for all our integration schemes.

Questions

  1. I do not believe our methods ever worked in "true" shu osher form, only in minimal storage form. We were lucky that for RK1 and RK2 the minimal storage and true shu osher forms are the same. Should we also implement non-minimal storage SSP integrators with shu osher alpha/beta matrices? (I am probably not going to do this now, but food for thought.)
  2. The integration scheme used in our examples and the burgers' benchmark are not totally inline with the minimal storage approach (see point 1). They're a hybrid between shu osher notation and minimal storage 2S* form. I am considering updating them to use the more general integrators, but this will imply significant changes to the task lists for those examples. How do people feel about this?
  3. I added some extra update functionality to leverage the integrators more easily. However, it's untested. Should we leave it out until I have a test case? Or should we add it now and fix it later after testing?
  4. One possibility here is we merge the integrators + the update functionality in this MR and I update one of the examples in a subsequent MR. Thoughts?

@pgrete your thoughts on this PR in particular will be valuable given you're using the low-storage integrators downstream.

PR Checklist

  • Code passes cpplint
  • New features are documented.
  • Adds a test for any bugs fixed. Adds tests for new features.
  • Code is formatted
  • Changes are summarized in CHANGELOG.md
  • CI has been triggered on Darwin for performance regression tests.
  • Docs build
  • (@lanl.gov employees) Update copyright on changed files

@Yurlungur Yurlungur added the enhancement New feature or request label Feb 24, 2023
@Yurlungur Yurlungur self-assigned this Feb 24, 2023
@Yurlungur
Copy link
Collaborator Author

I still need to write updated docs, which I have not done yet.

@Yurlungur
Copy link
Collaborator Author

This now pretty much ready for review.

@Yurlungur Yurlungur enabled auto-merge February 28, 2023 18:18
@Yurlungur Yurlungur disabled auto-merge February 28, 2023 18:18
@brryan
Copy link
Collaborator

brryan commented Mar 1, 2023

Kudos for making this change non-breaking. Only one comment now: isn't the template<> thing no longer needed in C++17?

@Yurlungur
Copy link
Collaborator Author

Oh yeah that's true---so maybe the Generic thing isn't needed in that case...

@Yurlungur Yurlungur requested a review from pdmullen March 9, 2023 17:14
@Yurlungur Yurlungur requested a review from BenWibking March 9, 2023 17:48
@bensworth
Copy link

Simple naming convention, but I would not call classes integrators if they don't do any integration and only contain data

@Yurlungur
Copy link
Collaborator Author

Simple naming convention, but I would not call classes integrators if they don't do any integration and only contain data

This is a fair point. I left the name because that's what it was called before I touched it. If others are open to a different naming convention, I can change the name.

@Yurlungur
Copy link
Collaborator Author

ping @pdmullen and @pgrete :)

Copy link
Collaborator

@pdmullen pdmullen left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks great, @jonahm-LANL 😸 I don't have too much to contribute... It looks like you mostly covered your bases with unit tests and CI is passing for this PR. More generally, here are some of the thoughts I had regarding your infrastructure and higher order, in general:

(1) Round-off error conservation gets a bit weird with higher order integrators according to the Athena++ community. See this issue: PrincetonUniversity/athena#300. Parthenon may want to address this in the future. As implemented, I predict that higher order integrators with non-exact coefficients will see some departures from conservation.
(2) I was initially thinking about how this PR might influence how one would add a sub-cycling integrator for split physics (e.g., super time stepping). However, I am now not entirely convinced that this infrastructure is what one would use. Instead, I wonder if the solvers infrastructure might instead be the better place for this kind of thing, e.g., for STS:

    auto &solver = tr[i].AddSuperTimeStep("RKL1", ...);
    auto sts = CreateTaskList(..., i, pmesh, tr, solver, mbase, miter, mout, dt);

doc/sphinx/src/integrators.rst Outdated Show resolved Hide resolved
doc/sphinx/src/integrators.rst Outdated Show resolved Hide resolved
doc/sphinx/src/integrators.rst Show resolved Hide resolved
tst/unit/test_unit_integrators.cpp Show resolved Hide resolved
@Yurlungur
Copy link
Collaborator Author

Ping one more person: @BenWibking @bprather @brryan @forrestglines @pgrete this PR needs one more approving review to go in. It has been waiting for review for almost six weeks. If I do not see any signs of life by next Monday I will force merge.

@Yurlungur
Copy link
Collaborator Author

Thanks for the review, @pdmullen ! Regarding your comments:

(1) Round-off error conservation gets a bit weird with higher order integrators according to the Athena++ community. See this issue: PrincetonUniversity/athena#300. Parthenon may want to address this in the future. As implemented, I predict that higher order integrators with non-exact coefficients will see some departures from conservation.

Yes---I'm sure this is correct. Something to keep an eye on for sure.

(2) I was initially thinking about how this PR might influence how one would add a sub-cycling integrator for split physics (e.g., super time stepping). However, I am now not entirely convinced that this infrastructure is what one would use. Instead, I wonder if the solvers infrastructure might instead be the better place for this kind of thing, e.g., for STS:

Yeah... I haven't thought carefully about this. I tried to make this infrastructure "ready" to extend for such things. But there's definitely questions to answer still as far as that goes.

Copy link
Collaborator

@bprather bprather left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using this in KHARMA and broke nothing. I like making MultiStageDriver a specialization of a new generic version.
I have no way to test the higher-order methods implemented here independently, but the unit tests seem complete.

@Yurlungur
Copy link
Collaborator Author

Thanks @pdmullen @bprather ! We now have two approvals. AthenaPK folks: @forrestglines @pgrete speak now-ish or forever hold your peace on the integrator MR. I'll enable auto-merge on Monday morning.

Copy link
Collaborator

@pgrete pgrete left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good and bonus points for the additional documentation.
I just found couple of typos that I just fixed directly.

@Yurlungur
Copy link
Collaborator Author

Thanks @pgrete ! In that case I'm going to pull the trigger.

@Yurlungur Yurlungur enabled auto-merge April 16, 2023 19:48
@Yurlungur Yurlungur merged commit b473118 into develop Apr 16, 2023
@Yurlungur Yurlungur deleted the jmm/rk4 branch April 16, 2023 20:57
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

8 participants