Using SUNDIALS without ifdefs for CPU/GPU?

[BenWibking]

I'd like to use SUNDIALS to integrate cooling source terms, very similar to what is done in Nyx (but I need to add additional cooling terms for our applications). I see a lot of ifdefs in the Nyx code that calls SUNDIALS depending on whether CUDA is enabled or not. The example codes here don't have explicit ifdefs, but maybe these examples are CPU-only.

Is there a wrapper available that hides the GPU details when using SUNDIALS?

[asalmgren]

Hey Ben, We're actually actively working on wrapping SUNDIALS time integrators in amrex to make it easier for applications to call them. Don Willcox is the one who can probably tell you the most ... we're all just coming off our holiday break so hopefully he'll be able to answer in the next day or two.

…

On Mon, Jan 3, 2022 at 3:27 PM Ben Wibking ***@***.***> wrote: I'd like to use SUNDIALS to integrate cooling source terms, very similar to what is done in Nyx (but I need to add additional cooling terms for our applications). I see a lot of ifdefs in the Nyx code that calls SUNDIALS depending on whether CUDA is enabled or not. The example codes here <https://github.com/AMReX-Codes/ATPESC-codes/tree/main/SUNDIALS%2BAMReX> don't have explicit ifdefs, but maybe these examples are CPU-only. Is there a wrapper available that hides the GPU details when using SUNDIALS? — Reply to this email directly, view it on GitHub <#2557>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACRE6YXCQ6MUTNJHOBJVFHTUUIWFHANCNFSM5LGGHCBA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

-- Ann Almgren Senior Scientist; Acting Dept Head, Applied Mathematics Pronouns: she/her/hers

[BenWibking]

Hi Ann,

That's good to hear. I would be very glad to try this out as soon as it's available, even in beta/preview form.

Cheers,
Ben

[dwillcox]

Hi @BenWibking, Like @asalmgren said, I'm working on such a wrapper. Currently it supports only the ARKODE integrators in SUNDIALS and is designed to take a user-supplied function that computes a MultiFab to hold the right-hand-side for an ODE system. I haven't implemented the same sort of wrapper for CVODE, which is the implicit integrator Nyx uses, but I don't think that would be terribly difficult.

So a couple of questions:

• Do you require CVODE to integrate your ODE system (and if so, I assume you need its implicit BDF method?)
• Is your code designed to (or easily able to) calculate a MultiFab containing the right-hand side of the ODE system?

[BenWibking]

I haven't written the code yet, but I was thinking of CVODE with BDF. Occasionally, it might be worthwhile to use the explicit integrator. I could calculate this and put it in a MultiFab - currently all of the hydro integration is done with local Fabs without extra MultiFabs.

[zingale]

also @BenWibking if you need fully implicit, and don't mind a zone-by-zone interface for the integrator, we have a C++ port of VODE that runs on GPUs that we use for reaction nets which are very stiff. It is in our Microphysics repo:
https://github.com/starkiller-astro/Microphysics/tree/main/integration/VODE

[BenWibking]

This also looks like a good option. In some cases, the reaction terms are not very stiff. Is there a way to switch to an explicit method, or is this only implicit?

[zingale]

we stripped out the explicit option. Our stuff is very stiff, and the number of RHS evaluations can vary a lot from one zone to the next. In your case, you might do best with the sundials option.

[drummerdoc]

Mike types faster than I do....I was going to offer a similar capability via the combustion codes via https://github.com/AMReX-Combustion/PelePhysics The version here packs a bunch of cells from a MultiFab into a SUNDIALS vector and integrates over time using either CVODE or ARKCODE/ERKSTEP (suggesting a lot of code duplication between Castro, Nyx and Pele ..:) ) [templated interface].

Incidentally, I see that the MFEM group has nicely wrapped all the SUNDIALS functionality (at least to some fairly recent release date). I'll bet that their approach could be duplicated for AMReX applications with data from MultiFabs, then Pele/Castro/Nyx and this other app could use all the integrator options without duplicating even more code. Note however that the implicit solvers do introduce a bit more complexity associated with the Jacobians...might be too much and not worth it in your app.

[BenWibking]

This looks like a good option. I'll take a look in more detail next week. @dwillcox is this similar to your wrapper?

[drummerdoc]

One of the distinctions between these different implementations might be in how the state is passed to and maintained by SUNDIALS. In the Pele codes, we come in with a standard multicomponent Fab from the MultiFab and "flatten" the initial state data into one long 1D array with all the cells and components from the fab stacked up. SUNDIALS does its thing and returns the result as a 1D array, and we "unflatten" that advanced state back into the Fab. There probably should rather be an opaque pointer way to pass this data around which may allow us to avoid the flattening bit (but there are strides data access issues anyway that we need to avoid....).

As Zingale mentioned above, this is one approach...giving SUNDIALS a batch of cells (e.g. all the cells in the Fab) to work with all at once. Another approach is to give SUNDIALS a smaller number of cells to work with at a time. On GPUs, we take the approach that each GPU thread will carry out all the work associated with a cell, and we advance many cells together in parallel. Where that really happens though is that our RHS function is called by SUNDIALS via a kernel launch, and each GPU thread in that launch goes into the SUNDIALS vector to grab and set the data associated with his cell. One could actually implement the parallelism in a number of ways - we did the simplest thing.

In this GPU approach, where lots of cells are integrated together in batch, they all advance at a sub cycled dt constrained by the worst case (stiffest) cell. There is an opportunity to exploit that by, say, partitioning cells by their difficulty or something. We are not doing that. We used to do that on CPU and it made a huge difference, but here it's not so simple to do that an get good GPU parallelism.

Anyway, in Pele ARKODE and CVODE are both compiled and available, and the decision is templated in the SUNDIALS solver constructor used to advance each Fab. We may, for example, decide at runtime that we want a batch of cells integrated with one or the other based on the cfl dt or something. Have a look at the "Reactions" folder in PelePhysics, and let me know if things don't make sense.

[BenWibking]

Cell-by-cell might be ideal. For other source terms (not the ones I want to use SUNDIALS for), I've written a Newton solver that does independent iterations cell-by-cell inside in a ParallelFor loop. This seems to give good performance for those terms, which can have a wide range of stiffnesses, without any kind of explicit partitioning into separate batches. There is thread divergence, but it seems to not really matter. Is this the kind of partitioning that you have tried, or are you thinking of something slightly different?

[BenWibking]

For these terms, though, it might not make as big of a difference, since they're usually not as stiff. I'll take a look, thanks!