Use of AMReX_TimeIntegrator.H with AmrLevel

[bezmi]

Are there any examples that show use of the TimeIntegrator class, particularly for the SUNDIALS adaptive methods when using AmrLevel with subcycling? I would like to experiment with integrating electron and ion dynamics at different rates in our multi fluid plasma code, where we have separate states for each fluid species and the field (same distribution map, however).

For example, see here for our RK2 implementation:

https://github.com/PlasmaSimUQ/cerberus/blob/a9930678ecd2788d6d6a4b9a2f91d6679d80413a/Source/MFP_advance.cpp#L58

The call to get_data is a FillPatch operation into RK_step[data_idx].U and RK_step[data_idx].dU stores the update. data_idx is the state index.

As a simple first attempt, I have tried this with only a single fluid state, with the aim of getting it working with the TimeIntegrator class before moving on to coupled states.

I load the state data into a Vector<MultiFab>, initialize the time integrator and create a rhs function that calls the calc_time_derivative (source terms) and calc_spatial_derivative (flux divergence) update functions. Using this with the amrex native RKIntegrator gives correct output for 'ForwardEuler' and first order 'RungeKutta', but only the first step in the higher order RK methods appears to do anything, with the data remaining unchanged for remaining steps of the RK update. Trying it with the integrator type set to SUNDIALS with appropriate method (ERK) doesn't run at all. I can post this code when I get back to the office if necessary.

Is this why AmrLevel::RK exists (relevant discussion topic maybe: #3319)? I am interested in the multi-rate methods provided in ARKODE, so using that isn't really an option in this case.

All advice is greatly appreciated :)

[asalmgren]

@ajnonaka @gardner48 -- can you address this one?

[bezmi]

Some further clarification and discussion after getting sleep:

but only the first step in the higher order RK methods appears to do anything, with the data remaining unchanged for remaining steps of the RK update.

This is not strictly true. Interior data is being updated during each step of the RKIntegrator update, however ghost data (which is what I was looking at when I made the post) remains unchanged throughout the RK update, and only changes when FillPatch is called before the next simulation time-step. This is to be expected, I believe. I would have to implement something like you've done in STvAR: https://github.com/AMReX-Astro/STvAR/blob/6fa69d6382165d3f04c2093e3ba5c4c0eebcb925/Source/AmrCoreAdv.cpp#L716 by supplying a function to set_post_update that fills the ghost data by interpolating in time when in between simulation time-steps.

I am still debugging why SUNDIALS won't work. My code runs with the behaviour described above when I use:

integrator.type = RungeKutta
integrator.method = 2

But when I run with:

integrator.type = SUNDIALS
integrator.sundials.type = ERK

# or any other method
integration.sundials.method = ARKODE_FORWARD_EULER_1_1

It crashes when trying to access ghost cells (negative indices) of the data Array4. Interior cells can be accessed just fine. I will update if I get it working.

[ajnonaka]

As far as I know there are no example codes that presently incorporate SUNDIALS into a multilevel hierarchy, even without subcycling. The STvAR example you cite looks like the closest thing - it was an attempt to do this but I was unaware that code existed until just now (and looks like nobody has touched that code in 2 years).

In the past year there have been quite a few updates to how SUNDIALS/AMReX works in terms of ghost cells, post_update, and other interface calls so I'd be surprised if the STvAR code even works or compiles at this point. We are certainly interested in adding multilevel support (libraries and examples) it's just not explicitly there yet.

[bezmi]

Thank you for your response, I'd assumed that may be the case. I'll keep chipping away at it and see where I get. I have detailed some more debugging progress below - could this be a bug in the wrap_data or unpack_vector routines, or am I just doing something obviously wrong? EDIT: a FillPatch call within rhs_fun fixes the issue below. The existing TimeIntegrator tutorials have this (FillBoundary), but I missed it. We have to make sure that U gets filled with the appropriate ghost cell data (interpolating in time if necessary) before calculating dU.

I use the following function as close to a minimum working example. It calls the ARKODE routines, but nothing is really done with the data. I am simply printing the values inside Vector<MultiFab> U:

void MFP::advance_arkode(Real time, Real dt, int iteration, int ncycle)
{
    BL_PROFILE("MFP::advance_arkode");
    Print() << "Begin advance_arkode\n";

    EulerianState& istate = EulerianState::get_state(0);

    // +2 for num_grow_eb
    int ng = istate.get_num_grow() + 2;

    int nc = istate.n_cons();

    Vector<MultiFab> old_data(1);
    Vector<MultiFab> new_data(1);

    old_data[0]
      .define(this->boxArray(), this->DistributionMap(), nc, ng, MFInfo(), this->Factory());

    new_data[0]
      .define(this->boxArray(), this->DistributionMap(), nc, ng, MFInfo(), this->Factory());

    // fill the ghost cells
    FillPatch(*this, old_data[0], ng, time, 0, 0, nc);

    // copy old_data into new_data
    MultiFab::Copy(new_data[0], old_data[0], 0, 0, nc, ng);

    Print() << "old_data[0].boxArray().numPts(): " << old_data[0].boxArray().numPts()
            << ", old_data[0].nGrow(): " << old_data[0].nGrow()
            << ", old_data[0].nComp(): " << old_data[0].nComp() << "\n";
    for (MFIter mfi(old_data[0]); mfi.isValid(); ++mfi) {
        const Box& box = mfi.tilebox();
        FArrayBox& old_data_fab = old_data[0][mfi];
        FArrayBox& new_data_fab = new_data[0][mfi];

        Array4<Real> const& old_data4 = old_data_fab.array();
        Array4<Real> const& new_data4 = new_data_fab.array();
        Print() << "old_data4.begin: " << old_data4.begin << ", old_data4.end: " << old_data4.end
                << "\n";
        Print() << "new_data4.begin: " << new_data4.begin << ", new_data4.end: " << new_data4.end
                << "\n";

        // try accessing some interior cell data (non-ghost)
        Print() << "old_data4(0,0,0,0): " << old_data4(0, 0, 0, 0) << "\n";

        // try accessing ghost cells
        Print() << "old_data4(-1, 0, 0, 0): " << old_data4(-1, 0, 0, 0) << "\n";
        Print() << "old_data4(-2, 0, 0, 0): " << old_data4(-2, 0, 0, 0) << "\n";
        Print() << "old_data4(-3, 0, 0, 0): " << old_data4(-3, 0, 0, 0) << "\n\n";

        Print() << "new_data4(0,0,0,0): " << new_data4(0, 0, 0, 0) << "\n";
        Print() << "new_data4(-1, 0, 0, 0): " << new_data4(-1, 0, 0, 0) << "\n";
        Print() << "new_data4(-2, 0, 0, 0): " << new_data4(-2, 0, 0, 0) << "\n";
        Print() << "new_data4(-3, 0, 0, 0): " << new_data4(-3, 0, 0, 0) << "\n\n";
    }

    auto rhs_fun = [&](Vector<MultiFab>& dU, Vector<MultiFab>& U, const Real time) {
        Print() << "rhs_fun() called\n";
        Print() << "dU[0].boxArray().numPts(): " << dU[0].boxArray().numPts()
                << ", dU[0].nGrow(): " << dU[0].nGrow() << ", dU[0].nComp(): " << dU[0].nComp()
                << "\n";
        Print() << "U[0].boxArray().numPts: " << U[0].boxArray().numPts()
                << ", U[0].nGrow(): " << U[0].nGrow() << ", U[0].nComp(): " << U[0].nComp() << "\n";

        for (MFIter mfi(U[0]); mfi.isValid(); ++mfi) {
            const Box& box = mfi.tilebox();
            FArrayBox& Ufab = U[0][mfi];

            Array4<Real> const& U4 = Ufab.array();
            Print() << "U4.begin: " << U4.begin << ", U4.end: " << U4.end << "\n";

            // try accessing some interior cell data (non-ghost)
            Print() << "U4(0, 0, 0, 0): " << U4(0, 0, 0, 0) << "\n";

            // try accessing ghost cells
            Print() << "U4(-1, 0, 0, 0): " << U4(-1, 0, 0, 0) << "\n";
            Print() << "U4(-2, 0, 0, 0): " << U4(-2, 0, 0, 0) << "\n";
            Print() << "U4(-3, 0, 0, 0): " << U4(-3, 0, 0, 0) << "\n\n";
        }
    };

    TimeIntegrator<Vector<MultiFab>> time_integrator(old_data);

    time_integrator.set_rhs(rhs_fun);
    time_integrator.set_time_step(dt);

    time_integrator.advance(old_data, new_data, time, dt);
}

And here is the console output for a single call to advance_arkode with amrex.fpe_trap_invalid=0 and integration.type = SUNDIALS, integration.sundials.type = ERK, integration.sundials.method = ARKODE_FORWARD_EULER_1_1.

Begin advance_arkode
old_data[0].boxArray().numPts(): 1024, old_data[0].nGrow(): 3, old_data[0].nComp(): 6
old_data4.begin: (-3,-3,0), old_data4.end: (35,35,1)
new_data4.begin: (-3,-3,0), new_data4.end: (35,35,1)
old_data4(0,0,0,0): 1
old_data4(-1, 0, 0, 0): 1
old_data4(-2, 0, 0, 0): 1
old_data4(-3, 0, 0, 0): 1

new_data4(0,0,0,0): 1
new_data4(-1, 0, 0, 0): 1
new_data4(-2, 0, 0, 0): 1
new_data4(-3, 0, 0, 0): 1

SUNDIALS ERK time integrator
SUNDIALS ERK method ARKODE_FORWARD_EULER_1_1
rhs_fun() called
dU[0].boxArray().numPts(): 1024, dU[0].nGrow(): 3, dU[0].nComp(): 6
U[0].boxArray().numPts: 1024, U[0].nGrow(): 3, U[0].nComp(): 6
U4.begin: (-3,-3,0), U4.end: (35,35,1)
U4(0, 0, 0, 0): 1
U4(-1, 0, 0, 0): nan
U4(-2, 0, 0, 0): nan
U4(-3, 0, 0, 0): nan

As you can see, the nGrow, begin and end for the multifabs passed through to rhs_fun are correct, however, when accessing ghost cell values, they are nan.

I have also tried this with MultiFab instead of Vector<MultiFab> in the relevant places. Ghost cells still come out as nan.

With the native RKIntegrator, the ghost cell values are printed correctly as 1.

[ajnonaka]

Thanks for the update. What you have done should fill ghost cells to something sensible but not formally correct since, as Weiqun put it (#3319), intermediate RK stages do not exactly have a time associated with it.

[APTX-4869-core]

Sorry, I want to ask a question that is not related to your issue.
Do you use both RK integrator and subcycling in your program?
Subcycling in the finer level has splited the Δt into ref_ratio parts, won't the addtional RK integrator casuse extra computation load?
Thank you in advance!

[bezmi]

We do use subcycling, but I don't understand why that would matter when it comes to choosing a higher order time integrator. The relative increase in the number of RHS function evaluations remains identical (eg, if you go from 1st order Euler to RK2) whether the code subcycles or not.

In our multi-fluid plasma equations, we solve the compressible Euler equations for separate electron and ion fluids. The ions are on the order of 1e3 times more massive. Due to this, the electrons have a much faster sound speed, so the timestep for a given dx will be limited by them (ignoring the electomagnetic fields and coupling for the sake of brevity). So, if we want to explicitly resolve the electron dynamics, we have to step the ion state at a much smaller dt than necessary for stability. With the ARKODE multi-rate integrators, I am hoping to be able to step electrons at a faster rate and ions at the slower rate to reduce the number of RHS function evaluations whilst maintaining approximately 2nd order time accuracy. The relative gain in effiiciency for this should be the same whether subcycling or not.

[APTX-4869-core]

Thanks for your reply!
Yes, using the explicit scheme needs small step in time.
At first, I think both the RK scheme and subcycling may cause extra conputation expense.
But after that, I re-think this problem, and I thought the RK won't result in much cost, since the RK doesn't contain complex step such as equation solving...
So I'm not doubting your method, I just want to know if the cost resulting from RK should be considered.
It seems that your solver must be high-order, so the RK step may only account for a small fraction of the computation cost and I know your purpose about this now.
Thank you for your reply again, and sorry for sparing your time!