A skeleton of the user’s main program

The following is a skeleton of the user’s main program (or calling program) for the integration of
an ODE IVP. Most of the steps are independent of the nvector, sunmatrix, sunlinsol, and
sunnonlinsol implementations used. For the steps that are not, refer to Chapters 7, 8, 9, and 10
for the specific name of the function to be called or macro to be referenced.

配合cvRoberts_dns.c来讲解

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

设置问题的维度大小

#define NEQ   3                /* number of equations  */

/* Create serial vector of length NEQ for I.C. and abstol */
y = N_VNew_Serial(NEQ);
if (check_retval((void *)y, "N_VNew_Serial", 0)) return(1);
abstol = N_VNew_Serial(NEQ); 
if (check_retval((void *)abstol, "N_VNew_Serial", 0)) return(1);

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 y 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 with a call of the
form ydata = N VGetArrayPointer(y0). See §7.3-7.6 for details.
For the hypre and petsc vector wrappers, first create and initialize the underlying vector, and
then create an 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 §7.7 and §7.8 for details.

初始化向量的值

#define Y1    RCONST(1.0)      /* initial y components */
#define Y2    RCONST(0.0)
#define Y3    RCONST(0.0)

 /* Create serial vector of length NEQ for I.C. and abstol */
  y = N_VNew_Serial(NEQ);
  if (check_retval((void *)y, "N_VNew_Serial", 0)) return(1);
  abstol = N_VNew_Serial(NEQ); 
  if (check_retval((void *)abstol, "N_VNew_Serial", 0)) return(1);

  /* Initialize y */
  Ith(y,1) = Y1;
  Ith(y,2) = Y2;
  Ith(y,3) = Y3;

4. Create cvode object

Call cvode mem = CVodeCreate(lmm) to create the cvode memory block and to specify the linear
multistep method. CVodeCreate returns a pointer to the cvode memory structure. See §4.5.1
for details.

申请求解器使用的内存空间

  /* Call CVodeCreate to create the solver memory and specify the 
   * Backward Differentiation Formula */
  cvode_mem = CVodeCreate(CV_BDF);
  if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);

5. Initialize cvode solver

Call CVodeInit(…) to provide required problem specifications, allocate internal memory for
cvode, and initialize cvode. CVodeInit returns a flag, the value of which indicates either
success or an illegal argument value. See §4.5.1 for details.

初始化cvode求解器

  /* Call CVodeInit to initialize the integrator memory and specify the
   * user's right hand side function in y'=f(t,y), the inital time T0, and
   * the initial dependent variable vector y. */
  retval = CVodeInit(cvode_mem, f, T0, y);
  if (check_retval(&retval, "CVodeInit", 1)) return(1);

6. Specify integration tolerances

Call CVodeSStolerances(…) or CVodeSVtolerances(…) to specify either a scalar relative
tolerance and scalar absolute tolerance, or a scalar relative tolerance and a vector of absolute
tolerances, respectively. Alternatively, call CVodeWFtolerances to specify a function which sets
directly the weights used in evaluating WRMS vector norms. See §4.5.2 for details.

设置积分的精度

  /* Set the scalar relative tolerance */
  reltol = RTOL;
  /* Set the vector absolute tolerance */
  Ith(abstol,1) = ATOL1;
  Ith(abstol,2) = ATOL2;
  Ith(abstol,3) = ATOL3;


  
  /* Call CVodeSVtolerances to specify the scalar relative tolerance
   * and vector absolute tolerances */
  retval = CVodeSVtolerances(cvode_mem, reltol, abstol);
  if (check_retval(&retval, "CVodeSVtolerances", 1)) return(1);

7. Create matrix object

If a nonlinear solver requiring a linear solve will be used (e.g., the default Newton iteration)
and the linear solver will be a matrix-based linear solver, then a template Jacobian matrix must
be created by calling the appropriate constructor function defined by the particular sunmatrix
implementation.
For the sundials-supplied sunmatrix implementations, the matrix object may be created using
a call of the form
SUNMatrix J = SUNBandMatrix(…);
or
SUNMatrix J = SUNDenseMatrix(…);
or
SUNMatrix J = SUNSparseMatrix(…);
NOTE: The dense, banded, and sparse matrix objects are usable only in a serial or threaded
environment.

创建计算使用到的矩阵

  /* Create dense SUNMatrix for use in linear solves */
  A = SUNDenseMatrix(NEQ, NEQ);
  if(check_retval((void *)A, "SUNDenseMatrix", 0)) return(1);

8. Create linear solver object

If a nonlinear solver requiring a linear solver is chosen (e.g., the default Newton iteration), then
the desired linear solver object must be created by calling the appropriate constructor function
defined by the particular sunlinsol implementation.
For any of the sundials-supplied sunlinsol implementations, the linear solver object may be
created using a call of the form
SUNLinearSolver LS = SUNLinSol *(…);
where * can be replaced with “Dense”, “SPGMR”, or other options, as discussed in §4.5.3 and
Chapter 9.

创建求解器

  /* Create dense SUNLinearSolver object for use by CVode */
  LS = SUNLinSol_Dense(y, A);
  if(check_retval((void *)LS, "SUNLinSol_Dense", 0)) return(1);

9. Set linear solver optional inputs

Call Set functions from the selected linear solver module to change optional inputs specific to
that linear solver. See the documentation for each sunlinsol module in Chapter 9 for details.

设置求解器的一些配置

 /* Call CVodeSetLinearSolver to attach the matrix and linear solver to CVode */
  retval = CVodeSetLinearSolver(cvode_mem, LS, A);
  if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);

10. Attach linear solver module

If a nonlinear solver requiring a linear solver is chosen (e.g., the default Newton iteration), then
initialize the cvls linear solver interface by attaching the linear solver object (and matrix object,
if applicable) with the call (for details see §4.5.3):
ier = CVodeSetLinearSolver(…);
Alternately, if the cvode-specific diagonal linear solver module, cvdiag, is desired, initialize the
linear solver module and attach it to cvode with the call
ier = CVDiag(…);

  /* Call CVodeSetLinearSolver to attach the matrix and linear solver to CVode */
  retval = CVodeSetLinearSolver(cvode_mem, LS, A);
  if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);

  /* Set the user-supplied Jacobian routine Jac */
  retval = CVodeSetJacFn(cvode_mem, Jac);
  if(check_retval(&retval, "CVodeSetJacFn", 1)) return(1);

11. Set optional inputs

Call CVodeSet* functions to change any optional inputs that control the behavior of cvode from
their default values. See §4.5.8.1 and §4.5.8 for details.

12-15是非线性的步骤，跟线性的一样

12. Create nonlinear solver object (optional )

If using a non-default nonlinear solver (see §4.5.4), then create the desired nonlinear solver object
by calling the appropriate constructor function defined by the particular sunnonlinsol imple-
mentation (e.g., NLS = SUNNonlinSol ***(…); where *** is the name of the nonlinear solver
(see Chapter 10 for details).

13. Attach nonlinear solver module (optional )

If using a non-default nonlinear solver, then initialize the nonlinear solver interface by attaching
the nonlinear solver object by calling ier = CVodeSetNonlinearSolver(cvode mem, NLS); (see
§4.5.4 for details).

14. Set nonlinear solver optional inputs (optional )

Call the appropriate set functions for the selected nonlinear solver module to change optional
inputs specific to that nonlinear solver. These must be called after CVodeInit if using the default
nonlinear solver or after attaching a new nonlinear solver to cvode, otherwise the optional inputs
will be overridden by cvode defaults. See Chapter 10 for more information on optional inputs.

15. Specify rootfinding problem (optional )

Call CVodeRootInit to initialize a rootfinding problem to be solved during the integration of the
ODE system. See §4.5.5, and see §4.5.8.3 for relevant optional input calls.

  /* Call CVodeRootInit to specify the root function g with 2 components */
  retval = CVodeRootInit(cvode_mem, 2, g);
  if (check_retval(&retval, "CVodeRootInit", 1)) return(1);

``


## 16. Advance solution in time
For each point at which output is desired, call ier = CVode(cvode 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 y(t). See §4.5.7 for details.
每步计算的输出通过引用传递参数yout返回。
itask参数表明了返回值的模式

```c++
    retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
    PrintOutput(t, Ith(y,1), Ith(y,2), Ith(y,3));

    if (retval == CV_ROOT_RETURN) {
      retvalr = CVodeGetRootInfo(cvode_mem, rootsfound);
      if (check_retval(&retvalr, "CVodeGetRootInfo", 1)) return(1);
      PrintRootInfo(rootsfound[0],rootsfound[1]);
    }

    if (check_retval(&retval, "CVode", 1)) break;
    if (retval == CV_SUCCESS) {
      iout++;
      tout *= TMULT;
    }

17. Get optional outputs

Call CVGet functions to obtain optional output. See §4.5.10 for details.

18. Deallocate memory for solution vector

Upon completion of the integration, deallocate memory for the vector y (or yout) by calling the
appropriate destructor function defined by the nvector implementation:
N VDestroy(y);

19. Free solver memory

Call CVodeFree(&cvode mem) to free the memory allocated by cvode.

20. Free nonlinear solver memory (optional )

##21. Free linear solver and matrix memory
Call SUNLinSolFree and SUNMatDestroy to free any memory allocated for the linear solver and
matrix objects created above.

##22. Finalize MPI, if used
Call MPI Finalize() to terminate MPI.