cvsReinit_ASAi_dns.c

cvsReinit_ASAi_dns.c

/* -----------------------------------------------------------------
 * Adjoint sensitivity example problem.
 * The following is a simple example problem, with the coding
 * needed for its solution by CVODES. 
 *    dy/dt = -p*y
 * on the interval from t = 0.0 to t = 4.0, with initial
 * conditions: y = 1.0. The reaction rates are:
 * p=0.05 
 * 
 * -----------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>

#include <cvodes/cvodes.h>             /* prototypes for CVODE fcts., consts.  */
#include <nvector/nvector_serial.h>    /* access to serial N_Vector            */
#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix            */
#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver      */
#include <cvodes/cvodes_direct.h>      /* access to CVDls interface            */
#include <sundials/sundials_types.h>   /* defs. of realtype, sunindextype      */
#include <sundials/sundials_math.h>    /* defs. of SUNRabs, SUNRexp, etc.      */

/* Accessor macros */

#define Ith(v,i)    NV_Ith_S(v,i-1)         /* i-th vector component, i=1..NEQ */
#define IJth(A,i,j) SM_ELEMENT_D(A,i-1,j-1) /* (i,j)-th matrix el., i,j=1..NEQ */

/* Problem Constants */

#define NEQ      1             /* number of equations                  */

#define RTOL     RCONST(1e-6)  /* scalar relative tolerance            */

#define ATOL1    RCONST(1e-8)  /* vector absolute tolerance components */

#define ATOLl    RCONST(1e-8)  /* absolute tolerance for adjoint vars. */
#define ATOLq    RCONST(1e-6)  /* absolute tolerance for quadratures   */

#define T0       RCONST(0.0)   /* initial time                         */
#define THEAVY   RCONST(2.0)   /* heavyside time                       */
#define TOUT     RCONST(4.0)   /* final time                           */

#define STEPS    1000           /* number of steps between check points */

#define NP       1             /* number of problem parameters         */

#define ZERO     RCONST(0.0)


/* Type : UserData */

typedef struct {
  realtype p[1];
} *UserData;

/* Prototypes of user-supplied functions */

static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
               void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
static int ewt(N_Vector y, N_Vector w, void *user_data);

/* Prototypes of private functions */

static void PrintHead(realtype tB0);
static void PrintOutput(realtype tfinal, N_Vector y, N_Vector yB, N_Vector qB);
static void PrintOutput1(realtype time, realtype t, N_Vector y, N_Vector yB);
static int check_flag(void *flagvalue, const char *funcname, int opt);

/*
 *--------------------------------------------------------------------
 * MAIN PROGRAM
 *--------------------------------------------------------------------
 */

int main(int argc, char *argv[])
{
  UserData data;

  SUNMatrix A;
  SUNLinearSolver LS;
  void *cvode_mem;

  N_Vector y, yres;

  int steps;

  int indexB;

  realtype time;
  int flag, ncheck;

  long int nst;

  CVadjCheckPointRec *ckpnt;

  data = NULL;
  A = NULL;
  LS = NULL;
  cvode_mem = NULL;
  ckpnt = NULL;
  y = NULL;
  yres = NULL;

  /* Print problem description */
  printf("\nAdjoint Sensitivity Example with Reinitialization\n");
  printf("-------------------------------------------------\n\n");
  printf("ODE: dy1/dt = -p*y1\n");

  /* User data structure */
  data = (UserData) malloc(sizeof *data);
  if (check_flag((void *)data, "malloc", 2)) return(1);
  data->p[0] = RCONST(0.05);

  /* Initialize y */
  y = N_VNew_Serial(NEQ);
  yres = N_VNew_Serial(NEQ);
  if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);
  Ith(y,1) = RCONST(1.0);

  /* Create and allocate CVODES memory for forward run */
  printf("Create and allocate CVODES memory for forward runs\n");

  /* Call CVodeCreate to create the solver memory and specify the 
     Backward Differentiation Formula and the use of a Newton iteration */
  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
  if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);

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

  /* Call CVodeWFtolerances to specify a user-supplied function ewt that sets
     the multiplicative error weights w_i for use in the weighted RMS norm */
  flag = CVodeWFtolerances(cvode_mem, ewt);
  if (check_flag(&flag, "CVodeWFtolerances", 1)) return(1);

  /* Attach user data */
  flag = CVodeSetUserData(cvode_mem, data);
  if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);

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

  /* Create dense SUNLinearSolver object */
  LS = SUNDenseLinearSolver(y, A);
  if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return(1);

  /* Attach the matrix and linear solver */
  flag = CVDlsSetLinearSolver(cvode_mem, LS, A);
  if (check_flag(&flag, "CVDlsSetLinearSolver", 1)) return(1);

  /* Set the user-supplied Jacobian routine Jac */
  flag = CVDlsSetJacFn(cvode_mem, Jac);
  if (check_flag(&flag, "CVDlsSetJacFn", 1)) return(1);

  /* Allocate global memory */

  /* Call CVodeAdjInit to update CVODES memory block by allocting the internal 
     memory needed for backward integration.*/
  steps = STEPS; /* no. of integration steps between two consecutive ckeckpoints*/
  flag = CVodeAdjInit(cvode_mem, steps, CV_POLYNOMIAL);
  if (check_flag(&flag, "CVodeAdjInit", 1)) return(1);

  /* Perform forward run */
  printf("Forward integration ... ");

  /* Call CVodeF to integrate the forward problem until THEAVY */
  flag = CVodeF(cvode_mem, THEAVY, y, &time, CV_NORMAL, &ncheck);
  if (check_flag(&flag, "CVodeF", 1)) return(1);

  /* Test check point linked list 
     (uncomment next block to print check point information) */
  
  flag = CVodeGetNumSteps(cvode_mem, &nst);
  if (check_flag(&flag, "CVodeGetNumSteps", 1)) return(1);

  {
    int i;
    int order;
    realtype timepoint;
    printf("\nList of AdjDataPointPolynomial (nst = %ld)\n\n", nst);
    
    for (i=0;i<=nst;i++) {
      CVodeGetAdjDataPointPolynomial(cvode_mem, i,  &timepoint, &order, yres);
      printf("t=%le, y=%le\n",timepoint,Ith(yres,1));
    }
    printf("\n");
  }

  printf("Reinitializing Solver ... \n");
  CVodeReInit(cvode_mem, time, y);


  /* Call CVodeF to integrate the forward problem over an interval in time and
     saves checkpointing data */
  flag = CVodeF(cvode_mem, TOUT, y, &time, CV_NORMAL, &ncheck);
  if (check_flag(&flag, "CVodeF", 1)) return(1);

  {
    int i;
    int order;
    realtype timepoint;
    printf("\nList of AdjDataPointPolynomial (nst = %ld)\n\n", nst);
    
    for (i=0;i<=nst;i++) {
      CVodeGetAdjDataPointPolynomial(cvode_mem, i,  &timepoint, &order, yres);
      printf("t=%le, y=%le\n",timepoint,Ith(yres,1));
    }
    printf("\n");
  }

  printf("done ( nst = %ld )\n",nst);
  printf("\nncheck = %d\n\n", ncheck);

  /* Free memory */
  printf("Free memory\n\n");

  CVodeFree(&cvode_mem);
  N_VDestroy(y); 
  SUNLinSolFree(LS);
  SUNMatDestroy(A);

  if (ckpnt != NULL) free(ckpnt);
  free(data);

  return(0);

}

/*
 *--------------------------------------------------------------------
 * FUNCTIONS CALLED BY CVODES
 *--------------------------------------------------------------------
 */

/*
 * f routine. Compute f(t,y). 
*/

static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
{
  realtype y1, yd1;
  UserData data;
  realtype p;

  y1 = Ith(y,1);
  data = (UserData) user_data;
  p = data->p[0];

  yd1 = Ith(ydot,1) = p*y1;

  return(0);
}

/* 
 * Jacobian routine. Compute J(t,y). 
*/

static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
               void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
{
  UserData data;
  realtype p;
 
  data = (UserData) user_data;
  p = data->p[0];
 
  IJth(J,1,1) = -p;

  return(0);
}

/*
 * EwtSet function. Computes the error weights at the current solution.
 */

static int ewt(N_Vector y, N_Vector w, void *user_data)
{
  int i;
  realtype yy, ww, rtol, atol[3];

  rtol    = RTOL;
  atol[0] = ATOL1;

  for (i=1; i<=1; i++) {
    yy = Ith(y,i);
    ww = rtol * SUNRabs(yy) + atol[i-1];
    if (ww <= 0.0) return (-1);
    Ith(w,i) = 1.0/ww;
  }

  return(0);
}

/*
 *--------------------------------------------------------------------
 * PRIVATE FUNCTIONS
 *--------------------------------------------------------------------
 */

/* 
 * Check function return value.
 *    opt == 0 means SUNDIALS function allocates memory so check if
 *             returned NULL pointer
 *    opt == 1 means SUNDIALS function returns a flag so check if
 *             flag >= 0
 *    opt == 2 means function allocates memory so check if returned
 *             NULL pointer 
 */

static int check_flag(void *flagvalue, const char *funcname, int opt)
{
  int *errflag;

  /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
  if (opt == 0 && flagvalue == NULL) {
    fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
	    funcname);
    return(1); }

  /* Check if flag < 0 */
  else if (opt == 1) {
    errflag = (int *) flagvalue;
    if (*errflag < 0) {
      fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
	      funcname, *errflag);
      return(1); }}

  /* Check if function returned NULL pointer - no memory allocated */
  else if (opt == 2 && flagvalue == NULL) {
    fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
	    funcname);
    return(1); }

  return(0);
}

Output:

Adjoint Sensitivity Example with Reinitialization

ODE: dy1/dt = -p*y1
Create and allocate CVODES memory for forward runs
Forward integration ...
List of AdjDataPointPolynomial (nst = 15)

t=0.000000e+00, y=1.000000e+00
t=1.421267e-02, y=1.000711e+00
t=2.842534e-02, y=1.001423e+00
t=6.714856e-02, y=1.003364e+00
t=1.058718e-01, y=1.005309e+00
t=1.445950e-01, y=1.007257e+00
t=1.833182e-01, y=1.009209e+00
t=2.596545e-01, y=1.013069e+00
t=3.875189e-01, y=1.019566e+00
t=5.153832e-01, y=1.026106e+00
t=7.661383e-01, y=1.039052e+00
t=1.016893e+00, y=1.052162e+00
t=1.267648e+00, y=1.065436e+00
t=1.518403e+00, y=1.078879e+00
t=1.988216e+00, y=1.104522e+00
t=3.312658e+00, y=1.180144e+00

Reinitializing Solver ...

List of AdjDataPointPolynomial (nst = 15)

t=0.000000e+00, y=1.000000e+00
t=2.014206e+00, y=1.105959e+00
t=2.028412e+00, y=1.106745e+00
t=2.067117e+00, y=1.108890e+00
t=2.105822e+00, y=1.111038e+00
t=2.144527e+00, y=1.113190e+00
t=2.183232e+00, y=1.115347e+00
t=2.259533e+00, y=1.119610e+00
t=2.387343e+00, y=1.126788e+00
t=2.515152e+00, y=1.134012e+00
t=2.765776e+00, y=1.148312e+00
t=3.016400e+00, y=1.162792e+00
t=3.267024e+00, y=1.177455e+00
t=3.517648e+00, y=1.192303e+00
t=3.987331e+00, y=1.220634e+00
t=5.312719e+00, y=1.304267e+00

done ( nst = 15 )

ncheck = 0