官方建议的ARKODE编程步骤

1. Initialize parallel or multi-threaded environment, if appropriate.
   For example, call MPI_Init to initialize MPI if used, or set num_threads, the number of threads to use
   within the threaded vector functions, if used.

---

2. Set problem dimensions, etc.
   This generally includes the problem size, N, and may include the local vector length Nlocal.
   Note: The variables N and Nlocal should be of type sunindextype.

---

3. Set vector of initial values
   To set the vector y0 of initial values, use the appropriate functions defined by the particular NVECTOR implementation.
   For native SUNDIALS vector implementations (except the CUDA and RAJA based ones), use a call of the form
   y0 = N_VMake_***(…, ydata);
   if the realtype array ydata containing the initial values of 𝑦 already exists. Otherwise, create a new vector
   by making a call of the form
   y0 = N_VNew_***(…);
   and then set its elements by accessing the underlying data where it is located with a call of the form
   ydata = N_VGetArrayPointer_***(y0);
   See the sections The NVECTOR_SERIAL Module through The NVECTOR_PTHREADS Module for details.
   For the HYPRE and PETSc vector wrappers, first create and initialize the underlying vector, and then create the
   NVECTOR wrapper with a call of the form
   y0 = N_VMake_***(yvec);
   where yvec is a HYPRE or PETSc vector. Note that calls like N_VNew_***(…) and
   N_VGetArrayPointer_***(…) are not available for these vector wrappers. See the sections The
   NVECTOR_PARHYP Module and The NVECTOR_PETSC Module for details.
   If using either the CUDA- or RAJA-based vector implementations use a call of the form
   y0 = N_VMake_***(…, c);
   where c is a pointer to a suncudavec or sunrajavec vector class if this class already exists. Otherwise,
   create a new vector by making a call of the form
   N_VGetDeviceArrayPointer_***
   or
   N_VGetHostArrayPointer_***
   Note that the vector class will allocate memory on both the host and device when instantiated. See the sections
   The NVECTOR_CUDA Module and The NVECTOR_RAJA Module for details.

---

4. Create an ARKStep object for the fast (inner) integration
   Call inner_arkode_mem = ARKStepCreate(…) to create the ARKStep memory block.
   ARKStepCreate() returns a void* pointer to this memory structure. See the section ARKStep initialization
   and deallocation functions for details.

---

5. Configure the fast integrator
   Specify tolerances, create and attach matrix and/or solver objects, or call ARKStepSet* functions to configure
   the fast integrator as desired. See sections A skeleton of the user’s main program and Optional input functions
   for details on configuring ARKStep.
   Notes on using ARKStep as a fast integrator:
   If the inner method is not explicitly specified then the default method in the ARKStep module will be used. If
   a particular fast method is desired it should be set in this step. The slow method can be set when configuring the
   slow integrator in the following steps.
   By default the fast integrator will use adaptive step sizes. To use a fixed fast step a call to
   ARKStepSetFixedStep() should be made in this step otherwise fast integration tolerances should be set
   in this step as described in A skeleton of the user’s main program.
   If attaching a user_data pointer, it should be attached to the slow integrator in the following steps with
   MRIStepSetUserData(). This pointer will subsequently be passed to user-provided functions during the
   fast integration.
   Specifying a rootfinding problem for the fast integration is not supported. Rootfinding problems should be
   created and initialized with the slow integrator. See the steps below and MRIStepRootInit() for more
   details.
   The ARKStep module used for the fast time scale must be configured with an identity mass matrix.

---

6. Create an MRIStep object for the slow (outer) integration
   Call arkode_mem = MRIStepCreate(…) to create the MRIStep memory block.
   MRIStepCreate() returns a void* pointer to this memory structure. See the section MRIStep initialization and deallocation functions for details.

---

7. Set the slow step size
   Call MRIStepSetFixedStep() to specify the slow time step size.

---

8. Set optional inputs
   Call MRIStepSet* functions to change any optional inputs that control the behavior of MRIStep from their
   default values. See the section Optional input functions for details.

---

9. Specify rootfinding problem
   Optionally, call MRIStepRootInit() to initialize a rootfinding problem to be solved during the integration
   of the ODE system. See the section Rootfinding initialization function for general details, and the section
   Optional input functions for relevant optional input calls.

---

10. Advance solution in time
    For each point at which output is desired, call
    ier = MRIStepEvolve(arkode_mem, tout, yout, &tret, itask);
    Here, itask specifies the return mode. The vector yout (which can be the same as the vector y0 above) will
    contain 𝑦(𝑡out). See the section MRIStep solver function for details.

---

11. Get optional outputs
    Call MRIStepGet* and/or ARKStepGet* functions to obtain optional output from the slow or fast integrators
    respectively. See the section Optional output functions and Optional output functions for details.

---

12. Deallocate memory for solution vector
    Upon completion of the integration, deallocate memory for the vector y (or yout) by calling the NVECTOR
    destructor function:
    N_VDestroy(y);

---

13. Free solver memory
    Call ARKStepFree(&inner_arkode_mem) and MRIStepFree(&arkode_mem) to free
    the memory allocated for fast and slow integration modules respectively.

---

14. Free linear solver and matrix memory
    Call SUNLinSolFree() and (possibly) SUNMatDestroy() to free any memory allocated for
    any linear solver and/or matrix objects created above for the fast integrator.

---

15. Finalize MPI, if used
    Call MPI_Finalize to terminate MPI.

---