SHARE
TWEET

Untitled

a guest May 25th, 2019 69 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /* -----------------------------------------------------------------
  2.  * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh, and
  3.  *                Radu Serban @ LLNL
  4.  * -----------------------------------------------------------------
  5.  * LLNS Copyright Start
  6.  * Copyright (c) 2017, Lawrence Livermore National Security
  7.  * This work was performed under the auspices of the U.S. Department
  8.  * of Energy by Lawrence Livermore National Laboratory in part under
  9.  * Contract W-7405-Eng-48 and in part under Contract DE-AC52-07NA27344.
  10.  * Produced at the Lawrence Livermore National Laboratory.
  11.  * All rights reserved.
  12.  * For details, see the LICENSE file.
  13.  * LLNS Copyright End
  14.  * -----------------------------------------------------------------
  15.  * Example problem:
  16.  *
  17.  * The following is a simple example problem, with the coding
  18.  * needed for its solution by CVODES for Forward Sensitivity
  19.  * Analysis. The problem is from chemical kinetics, and consists
  20.  * of the following three rate equations:
  21.  *    dy1/dt = -p1*y1 + p2*y2*y3
  22.  *    dy2/dt =  p1*y1 - p2*y2*y3 - p3*(y2)^2
  23.  *    dy3/dt =  p3*(y2)^2
  24.  * on the interval from t = 0.0 to t = 4.e10, with initial
  25.  * conditions y1 = 1.0, y2 = y3 = 0. The reaction rates are: p1=0.04,
  26.  * p2=1e4, and p3=3e7. The problem is stiff.
  27.  * This program solves the problem with the BDF method, Newton
  28.  * iteration with the DENSE linear solver, and a
  29.  * user-supplied Jacobian routine.
  30.  * It uses a scalar relative tolerance and a vector absolute
  31.  * tolerance.
  32.  * Output is printed in decades from t = .4 to t = 4.e10.
  33.  * Run statistics (optional outputs) are printed at the end.
  34.  *
  35.  * Optionally, CVODES can compute sensitivities with respect to the
  36.  * problem parameters p1, p2, and p3.
  37.  * The sensitivity right hand side is given analytically through the
  38.  * user routine fS (of type SensRhs1Fn).
  39.  * Any of three sensitivity methods (SIMULTANEOUS, STAGGERED, and
  40.  * STAGGERED1) can be used and sensitivities may be included in the
  41.  * error test or not (error control set on SUNTRUE or SUNFALSE,
  42.  * respectively).
  43.  *
  44.  * Execution:
  45.  *
  46.  * If no sensitivities are desired:
  47.  *    % cvsRoberts_FSA_dns -nosensi
  48.  * If sensitivities are to be computed:
  49.  *    % cvsRoberts_FSA_dns -sensi sensi_meth err_con
  50.  * where sensi_meth is one of {sim, stg, stg1} and err_con is one of
  51.  * {t, f}.
  52.  * -----------------------------------------------------------------*/
  53.  
  54. #include <stdio.h>
  55. #include <stdlib.h>
  56. #include <string.h>
  57.  
  58. #include <cvodes/cvodes.h>             /* prototypes for CVODE fcts., consts.  */
  59. #include <nvector/nvector_serial.h>    /* access to serial N_Vector            */
  60. #include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix            */
  61. #include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver      */
  62. #include <sundials/sundials_types.h>   /* defs. of realtype, sunindextype      */
  63. #include <sundials/sundials_math.h>    /* definition of ABS */
  64.  
  65.  
  66. #include <iostream>
  67. #include <vector>
  68.  
  69. /* Accessor macros */
  70.  
  71. #define Ith(v,i)    NV_Ith_S(v,i-1)         /* i-th vector component i=1..NEQ */
  72. #define IJth(A,i,j) SM_ELEMENT_D(A,i-1,j-1) /* (i,j)-th matrix component i,j=1..NEQ */
  73.  
  74. /* Problem Constants */
  75.  
  76. #define NEQ   3             /* number of equations  */
  77. #define Y1    RCONST(1.0)   /* initial y components */
  78. #define Y2    RCONST(0.0)
  79. #define Y3    RCONST(0.0)
  80. #define RTOL  RCONST(1e-4)  /* scalar relative tolerance */
  81. #define ATOL1 RCONST(1e-8)  /* vector absolute tolerance components */
  82. #define ATOL2 RCONST(1e-14)
  83. #define ATOL3 RCONST(1e-6)
  84. #define T0    RCONST(0.0)   /* initial time */
  85. #define T1    RCONST(0.4)   /* first output time */
  86. #define TMULT RCONST(10.0)  /* output time factor */
  87. #define NOUT  12            /* number of output times */
  88.  
  89. #define NP    3             /* number of problem parameters */
  90. #define NS    3             /* number of sensitivities computed */
  91.  
  92. #define ZERO  RCONST(0.0)
  93.  
  94. /* Type : UserData */
  95.  
  96. typedef struct {
  97.   // realtype p[3];           /* problem parameters */
  98.   std::vector<double> p;
  99. } *UserData;
  100.  
  101. /* Prototypes of functions by CVODES */
  102.  
  103. static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
  104.  
  105. static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
  106.                void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
  107.  
  108. static int fS(int Ns, realtype t, N_Vector y, N_Vector ydot,
  109.               int iS, N_Vector yS, N_Vector ySdot,
  110.               void *user_data, N_Vector tmp1, N_Vector tmp2);
  111.  
  112. static int ewt(N_Vector y, N_Vector w, void *user_data);
  113.  
  114. /* Prototypes of private functions */
  115.  
  116. static void ProcessArgs(int argc, char *argv[],
  117.                         booleantype *sensi, int *sensi_meth,
  118.                         booleantype *err_con);
  119. static void WrongArgs(char *name);
  120. static void PrintOutput(void *cvode_mem, realtype t, N_Vector u);
  121. static void PrintOutputS(N_Vector *uS);
  122. static void PrintFinalStats(void *cvode_mem, booleantype sensi);
  123. static int check_retval(void *returnvalue, const char *funcname, int opt);
  124.  
  125. /*
  126.  *--------------------------------------------------------------------
  127.  * MAIN PROGRAM
  128.  *--------------------------------------------------------------------
  129.  */
  130.  
  131. int main(int argc, char *argv[])
  132. {
  133.   SUNMatrix A;
  134.   SUNLinearSolver LS;
  135.   void *cvode_mem;
  136.   UserData data;
  137.   realtype t, tout;
  138.   N_Vector y;
  139.   int iout, retval;
  140.  
  141.   realtype pbar[NS];
  142.   int is;
  143.   N_Vector *yS;
  144.   booleantype sensi, err_con;
  145.   int sensi_meth;
  146.  
  147.   cvode_mem = NULL;
  148.   data      = NULL;
  149.   y         = NULL;
  150.   yS        = NULL;
  151.   A         = NULL;
  152.   LS        = NULL;
  153.  
  154.   /* Process arguments */
  155.   ProcessArgs(argc, argv, &sensi, &sensi_meth, &err_con);
  156.  
  157.   /* User data structure */
  158.   data = (UserData) malloc(sizeof *data);
  159.   if (check_retval((void *)data, "malloc", 2)) return(1);
  160.   // data->p[0] = RCONST(0.04);
  161.   // data->p[1] = RCONST(1.0e4);
  162.   // data->p[2] = RCONST(3.0e7);
  163.  
  164.   data->p.push_back(RCONST(0.04));
  165.   data->p.push_back(RCONST(1.0e4));
  166.   data->p.push_back(RCONST(3.0e7));
  167.  
  168.   /* Initial conditions */
  169.   y = N_VNew_Serial(NEQ);
  170.   if (check_retval((void *)y, "N_VNew_Serial", 0)) return(1);
  171.  
  172.   Ith(y,1) = Y1;
  173.   Ith(y,2) = Y2;
  174.   Ith(y,3) = Y3;
  175.  
  176.   /* Create CVODES object */
  177.   cvode_mem = CVodeCreate(CV_BDF);
  178.   if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) return(1);
  179.  
  180.   /* Allocate space for CVODES */
  181.   retval = CVodeInit(cvode_mem, f, T0, y);
  182.   if (check_retval(&retval, "CVodeInit", 1)) return(1);
  183.  
  184.   /* Use private function to compute error weights */
  185.   retval = CVodeWFtolerances(cvode_mem, ewt);
  186.   if (check_retval(&retval, "CVodeSetEwtFn", 1)) return(1);
  187.  
  188.   /* Attach user data */
  189.   retval = CVodeSetUserData(cvode_mem, data);
  190.   if (check_retval(&retval, "CVodeSetUserData", 1)) return(1);
  191.  
  192.   /* Create dense SUNMatrix */
  193.   A = SUNDenseMatrix(NEQ, NEQ);
  194.   if (check_retval((void *)A, "SUNDenseMatrix", 0)) return(1);
  195.  
  196.   /* Create dense SUNLinearSolver */
  197.   LS = SUNLinSol_Dense(y, A);
  198.   if (check_retval((void *)LS, "SUNLinSol_Dense", 0)) return(1);
  199.  
  200.   /* Attach the matrix and linear solver */
  201.   retval = CVodeSetLinearSolver(cvode_mem, LS, A);
  202.   if (check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);
  203.  
  204.   /* Set the user-supplied Jacobian routine Jac */
  205.   /*retval = CVodeSetJacFn(cvode_mem, Jac);
  206.   if (check_retval(&retval, "CVodeSetJacFn", 1)) return(1);
  207.   */
  208.   printf("\n3-species chemical kinetics problem\n");
  209.  
  210.   /* Sensitivity-related settings */
  211.   if (sensi) {
  212.  
  213.     /* Set parameter scaling factor */
  214.     pbar[0] = data->p.at(0);
  215.     pbar[1] = data->p.at(1);
  216.     pbar[2] = data->p.at(2);
  217.  
  218.     /* Set sensitivity initial conditions */
  219.     yS = N_VCloneVectorArray(NS, y);
  220.     if (check_retval((void *)yS, "N_VCloneVectorArray", 0)) return(1);
  221.     for (is=0;is<NS;is++) N_VConst(ZERO, yS[is]);
  222.  
  223.     /* Call CVodeSensInit1 to activate forward sensitivity computations
  224.        and allocate internal memory for COVEDS related to sensitivity
  225.        calculations. Computes the right-hand sides of the sensitivity
  226.        ODE, one at a time */
  227.     retval = CVodeSensInit1(cvode_mem, NS, sensi_meth, NULL, yS);   // fS is third argument
  228.     if(check_retval(&retval, "CVodeSensInit", 1)) return(1);
  229.    
  230.  
  231.     /* Call CVodeSensEEtolerances to estimate tolerances for sensitivity
  232.        variables based on the rolerances supplied for states variables and
  233.        the scaling factor pbar */
  234.     retval = CVodeSensEEtolerances(cvode_mem);
  235.     if(check_retval(&retval, "CVodeSensEEtolerances", 1)) return(1);
  236.    
  237.     /* Set sensitivity analysis optional inputs */
  238.     /* Call CVodeSetSensErrCon to specify the error control strategy for
  239.        sensitivity variables */
  240.     retval = CVodeSetSensErrCon(cvode_mem, err_con);
  241.     if (check_retval(&retval, "CVodeSetSensErrCon", 1)) return(1);
  242.    
  243.     /* Call CVodeSetSensParams to specify problem parameter information for
  244.        sensitivity calculations */
  245.     /*
  246.     int parameter_len = 3;
  247.     realtype *param_array;
  248.     param_array = (realtype*)malloc(parameter_len * sizeof(realtype));
  249.     param_array[0] = RCONST(0.04);
  250.     param_array[1] = RCONST(1.0e4);
  251.     param_array[2] = RCONST(3.0e7);
  252.     */
  253.     retval = CVodeSetSensParams(cvode_mem, data->p[0], pbar, NULL);              // data->p
  254.     if (check_retval(&retval, "CVodeSetSensParams", 1)) return(1);
  255.    
  256.     printf("Sensitivity: YES ");
  257.     if(sensi_meth == CV_SIMULTANEOUS)  
  258.       printf("( SIMULTANEOUS +");
  259.     else
  260.       if(sensi_meth == CV_STAGGERED) printf("( STAGGERED +");
  261.       else                           printf("( STAGGERED1 +");  
  262.     if(err_con) printf(" FULL ERROR CONTROL )");
  263.     else        printf(" PARTIAL ERROR CONTROL )");
  264.  
  265.   } else {
  266.  
  267.     printf("Sensitivity: NO ");
  268.  
  269.   }
  270.  
  271.   /* In loop over output points, call CVode, print results, test for error */
  272.  
  273.   printf("\n\n");
  274.   printf("===========================================");
  275.   printf("============================\n");
  276.   printf("     T     Q       H      NST           y1");
  277.   printf("           y2           y3    \n");
  278.   printf("===========================================");
  279.   printf("============================\n");
  280.  
  281.   for (iout=1, tout=T1; iout <= NOUT; iout++, tout *= TMULT) {
  282.  
  283.     retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
  284.     if (check_retval(&retval, "CVode", 1)) break;
  285.  
  286.     PrintOutput(cvode_mem, t, y);
  287.  
  288.     /* Call CVodeGetSens to get the sensitivity solution vector after a
  289.        successful return from CVode */
  290.     if (sensi) {
  291.       retval = CVodeGetSens(cvode_mem, &t, yS);
  292.       if (check_retval(&retval, "CVodeGetSens", 1)) break;
  293.       PrintOutputS(yS);
  294.     }
  295.     printf("-----------------------------------------");
  296.     printf("------------------------------\n");
  297.  
  298.   }
  299.  
  300.   /* Print final statistics */
  301.   PrintFinalStats(cvode_mem, sensi);
  302.  
  303.   /* Free memory */
  304.  
  305.   N_VDestroy(y);                    /* Free y vector */
  306.   if (sensi) {
  307.     N_VDestroyVectorArray(yS, NS);  /* Free yS vector */
  308.   }
  309.   free(data);                              /* Free user data */
  310.   CVodeFree(&cvode_mem);                   /* Free CVODES memory */
  311.   SUNLinSolFree(LS);                       /* Free the linear solver memory */
  312.   SUNMatDestroy(A);                        /* Free the matrix memory */
  313.  
  314.   return(0);
  315. }
  316.  
  317. /*
  318.  *--------------------------------------------------------------------
  319.  * FUNCTIONS CALLED BY CVODES
  320.  *--------------------------------------------------------------------
  321.  */
  322.  
  323. /*
  324.  * f routine. Compute f(t,y).
  325.  */
  326.  
  327. static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
  328. {
  329.   realtype y1, y2, y3, yd1, yd3;
  330.   UserData data;
  331.   realtype p1, p2, p3;
  332.  
  333.   y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
  334.   data = (UserData) user_data;
  335.   p1 = data->p[0]; p2 = data->p[1]; p3 = data->p[2];
  336.  
  337.   yd1 = Ith(ydot,1) = -p1*y1 + p2*y2*y3;
  338.   yd3 = Ith(ydot,3) = p3*y2*y2;
  339.         Ith(ydot,2) = -yd1 - yd3;
  340.  
  341.   return(0);
  342. }
  343.  
  344.  
  345. /*
  346.  * Jacobian routine. Compute J(t,y).
  347.  */
  348.  
  349. static int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J,
  350.                void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
  351. {
  352.   realtype y2, y3;
  353.   UserData data;
  354.   realtype p1, p2, p3;
  355.  
  356.   y2 = Ith(y,2); y3 = Ith(y,3);
  357.   data = (UserData) user_data;
  358.   p1 = data->p[0]; p2 = data->p[1]; p3 = data->p[2];
  359.  
  360.   IJth(J,1,1) = -p1;  IJth(J,1,2) = p2*y3;          IJth(J,1,3) = p2*y2;
  361.   IJth(J,2,1) =  p1;  IJth(J,2,2) = -p2*y3-2*p3*y2; IJth(J,2,3) = -p2*y2;
  362.                       IJth(J,3,2) = 2*p3*y2;
  363.  
  364.   return(0);
  365. }
  366.  
  367. /*
  368.  * fS routine. Compute sensitivity r.h.s.
  369.  */
  370.  
  371. static int fS(int Ns, realtype t, N_Vector y, N_Vector ydot,
  372.               int iS, N_Vector yS, N_Vector ySdot,
  373.               void *user_data, N_Vector tmp1, N_Vector tmp2)
  374. {
  375.   UserData data;
  376.   realtype p1, p2, p3;
  377.   realtype y1, y2, y3;
  378.   realtype s1, s2, s3;
  379.   realtype sd1, sd2, sd3;
  380.  
  381.   data = (UserData) user_data;
  382.   p1 = data->p[0]; p2 = data->p[1]; p3 = data->p[2];
  383.  
  384.   y1 = Ith(y,1);  y2 = Ith(y,2);  y3 = Ith(y,3);
  385.   s1 = Ith(yS,1); s2 = Ith(yS,2); s3 = Ith(yS,3);
  386.  
  387.   sd1 = -p1*s1 + p2*y3*s2 + p2*y2*s3;
  388.   sd3 = 2*p3*y2*s2;
  389.   sd2 = -sd1-sd3;
  390.  
  391.   switch (iS) {
  392.   case 0:
  393.     sd1 += -y1;
  394.     sd2 +=  y1;
  395.     break;
  396.   case 1:
  397.     sd1 +=  y2*y3;
  398.     sd2 += -y2*y3;
  399.     break;
  400.   case 2:
  401.     sd2 += -y2*y2;
  402.     sd3 +=  y2*y2;
  403.     break;
  404.   }
  405.  
  406.   Ith(ySdot,1) = sd1;
  407.   Ith(ySdot,2) = sd2;
  408.   Ith(ySdot,3) = sd3;
  409.  
  410.   return(0);
  411. }
  412.  
  413. /*
  414.  * EwtSet function. Computes the error weights at the current solution.
  415.  */
  416.  
  417. static int ewt(N_Vector y, N_Vector w, void *user_data)
  418. {
  419.   int i;
  420.   realtype yy, ww, rtol, atol[3];
  421.  
  422.   rtol    = RTOL;
  423.   atol[0] = ATOL1;
  424.   atol[1] = ATOL2;
  425.   atol[2] = ATOL3;
  426.  
  427.   for (i=1; i<=3; i++) {
  428.     yy = Ith(y,i);
  429.     ww = rtol * SUNRabs(yy) + atol[i-1];
  430.     if (ww <= 0.0) return (-1);
  431.     Ith(w,i) = 1.0/ww;
  432.   }
  433.  
  434.   return(0);
  435. }
  436.  
  437. /*
  438.  *--------------------------------------------------------------------
  439.  * PRIVATE FUNCTIONS
  440.  *--------------------------------------------------------------------
  441.  */
  442.  
  443. /*
  444.  * Process and verify arguments to cvsfwddenx.
  445.  */
  446.  
  447. static void ProcessArgs(int argc, char *argv[],
  448.                         booleantype *sensi, int *sensi_meth, booleantype *err_con)
  449. {
  450.   *sensi = SUNFALSE;
  451.   *sensi_meth = -1;
  452.   *err_con = SUNFALSE;
  453.  
  454.   if (argc < 2) WrongArgs(argv[0]);
  455.  
  456.   if (strcmp(argv[1],"-nosensi") == 0)
  457.     *sensi = SUNFALSE;
  458.   else if (strcmp(argv[1],"-sensi") == 0)
  459.     *sensi = SUNTRUE;
  460.   else
  461.     WrongArgs(argv[0]);
  462.  
  463.   if (*sensi) {
  464.  
  465.     if (argc != 4)
  466.       WrongArgs(argv[0]);
  467.  
  468.     if (strcmp(argv[2],"sim") == 0)
  469.       *sensi_meth = CV_SIMULTANEOUS;
  470.     else if (strcmp(argv[2],"stg") == 0)
  471.       *sensi_meth = CV_STAGGERED;
  472.     else if (strcmp(argv[2],"stg1") == 0)
  473.       *sensi_meth = CV_STAGGERED1;
  474.     else
  475.       WrongArgs(argv[0]);
  476.  
  477.     if (strcmp(argv[3],"t") == 0)
  478.       *err_con = SUNTRUE;
  479.     else if (strcmp(argv[3],"f") == 0)
  480.       *err_con = SUNFALSE;
  481.     else
  482.       WrongArgs(argv[0]);
  483.   }
  484.  
  485. }
  486.  
  487. static void WrongArgs(char *name)
  488. {
  489.     printf("\nUsage: %s [-nosensi] [-sensi sensi_meth err_con]\n",name);
  490.     printf("         sensi_meth = sim, stg, or stg1\n");
  491.     printf("         err_con    = t or f\n");
  492.    
  493.     exit(0);
  494. }
  495.  
  496. /*
  497.  * Print current t, step count, order, stepsize, and solution.
  498.  */
  499.  
  500. static void PrintOutput(void *cvode_mem, realtype t, N_Vector u)
  501. {
  502.   long int nst;
  503.   int qu, retval;
  504.   realtype hu, *udata;
  505.  
  506.   udata = N_VGetArrayPointer(u);
  507.  
  508.   retval = CVodeGetNumSteps(cvode_mem, &nst);
  509.   check_retval(&retval, "CVodeGetNumSteps", 1);
  510.   retval = CVodeGetLastOrder(cvode_mem, &qu);
  511.   check_retval(&retval, "CVodeGetLastOrder", 1);
  512.   retval = CVodeGetLastStep(cvode_mem, &hu);
  513.   check_retval(&retval, "CVodeGetLastStep", 1);
  514.  
  515. #if defined(SUNDIALS_EXTENDED_PRECISION)
  516.   printf("%8.3Le %2d  %8.3Le %5ld\n", t, qu, hu, nst);
  517. #elif defined(SUNDIALS_DOUBLE_PRECISION)
  518.   printf("%8.3e %2d  %8.3e %5ld\n", t, qu, hu, nst);
  519. #else
  520.   printf("%8.3e %2d  %8.3e %5ld\n", t, qu, hu, nst);
  521. #endif
  522.  
  523.   printf("                  Solution       ");
  524.  
  525. #if defined(SUNDIALS_EXTENDED_PRECISION)
  526.   printf("%12.4Le %12.4Le %12.4Le \n", udata[0], udata[1], udata[2]);
  527. #elif defined(SUNDIALS_DOUBLE_PRECISION)
  528.   printf("%12.4e %12.4e %12.4e \n", udata[0], udata[1], udata[2]);
  529. #else
  530.   printf("%12.4e %12.4e %12.4e \n", udata[0], udata[1], udata[2]);
  531. #endif
  532.  
  533. }
  534.  
  535. /*
  536.  * Print sensitivities.
  537. */
  538.  
  539. static void PrintOutputS(N_Vector *uS)
  540. {
  541.   realtype *sdata;
  542.  
  543.   sdata = N_VGetArrayPointer(uS[0]);
  544.   printf("                  Sensitivity 1  ");
  545.  
  546. #if defined(SUNDIALS_EXTENDED_PRECISION)
  547.   printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
  548. #elif defined(SUNDIALS_DOUBLE_PRECISION)
  549.   printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  550. #else
  551.   printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  552. #endif
  553.  
  554.   sdata = N_VGetArrayPointer(uS[1]);
  555.   printf("                  Sensitivity 2  ");
  556.  
  557. #if defined(SUNDIALS_EXTENDED_PRECISION)
  558.   printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
  559. #elif defined(SUNDIALS_DOUBLE_PRECISION)
  560.   printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  561. #else
  562.   printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  563. #endif
  564.  
  565.   sdata = N_VGetArrayPointer(uS[2]);
  566.   printf("                  Sensitivity 3  ");
  567.  
  568. #if defined(SUNDIALS_EXTENDED_PRECISION)
  569.   printf("%12.4Le %12.4Le %12.4Le \n", sdata[0], sdata[1], sdata[2]);
  570. #elif defined(SUNDIALS_DOUBLE_PRECISION)
  571.   printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  572. #else
  573.   printf("%12.4e %12.4e %12.4e \n", sdata[0], sdata[1], sdata[2]);
  574. #endif
  575. }
  576.  
  577. /*
  578.  * Print some final statistics from the CVODES memory.
  579.  */
  580.  
  581. static void PrintFinalStats(void *cvode_mem, booleantype sensi)
  582. {
  583.   long int nst;
  584.   long int nfe, nsetups, nni, ncfn, netf;
  585.   long int nfSe, nfeS, nsetupsS, nniS, ncfnS, netfS;
  586.   long int nje, nfeLS;
  587.   int retval;
  588.  
  589.   retval = CVodeGetNumSteps(cvode_mem, &nst);
  590.   check_retval(&retval, "CVodeGetNumSteps", 1);
  591.   retval = CVodeGetNumRhsEvals(cvode_mem, &nfe);
  592.   check_retval(&retval, "CVodeGetNumRhsEvals", 1);
  593.   retval = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
  594.   check_retval(&retval, "CVodeGetNumLinSolvSetups", 1);
  595.   retval = CVodeGetNumErrTestFails(cvode_mem, &netf);
  596.   check_retval(&retval, "CVodeGetNumErrTestFails", 1);
  597.   retval = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
  598.   check_retval(&retval, "CVodeGetNumNonlinSolvIters", 1);
  599.   retval = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
  600.   check_retval(&retval, "CVodeGetNumNonlinSolvConvFails", 1);
  601.  
  602.   if (sensi) {
  603.     retval = CVodeGetSensNumRhsEvals(cvode_mem, &nfSe);
  604.     check_retval(&retval, "CVodeGetSensNumRhsEvals", 1);
  605.     retval = CVodeGetNumRhsEvalsSens(cvode_mem, &nfeS);
  606.     check_retval(&retval, "CVodeGetNumRhsEvalsSens", 1);
  607.     retval = CVodeGetSensNumLinSolvSetups(cvode_mem, &nsetupsS);
  608.     check_retval(&retval, "CVodeGetSensNumLinSolvSetups", 1);
  609.     retval = CVodeGetSensNumErrTestFails(cvode_mem, &netfS);
  610.     check_retval(&retval, "CVodeGetSensNumErrTestFails", 1);
  611.     retval = CVodeGetSensNumNonlinSolvIters(cvode_mem, &nniS);
  612.     check_retval(&retval, "CVodeGetSensNumNonlinSolvIters", 1);
  613.     retval = CVodeGetSensNumNonlinSolvConvFails(cvode_mem, &ncfnS);
  614.     check_retval(&retval, "CVodeGetSensNumNonlinSolvConvFails", 1);
  615.   }
  616.  
  617.   retval = CVodeGetNumJacEvals(cvode_mem, &nje);
  618.   check_retval(&retval, "CVodeGetNumJacEvals", 1);
  619.   retval = CVodeGetNumLinRhsEvals(cvode_mem, &nfeLS);
  620.   check_retval(&retval, "CVodeGetNumLinRhsEvals", 1);
  621.  
  622.   printf("\nFinal Statistics\n\n");
  623.   printf("nst     = %5ld\n\n", nst);
  624.   printf("nfe     = %5ld\n",   nfe);
  625.   printf("netf    = %5ld    nsetups  = %5ld\n", netf, nsetups);
  626.   printf("nni     = %5ld    ncfn     = %5ld\n", nni, ncfn);
  627.  
  628.   if(sensi) {
  629.     printf("\n");
  630.     printf("nfSe    = %5ld    nfeS     = %5ld\n", nfSe, nfeS);
  631.     printf("netfs   = %5ld    nsetupsS = %5ld\n", netfS, nsetupsS);
  632.     printf("nniS    = %5ld    ncfnS    = %5ld\n", nniS, ncfnS);
  633.   }
  634.  
  635.   printf("\n");
  636.   printf("nje    = %5ld    nfeLS     = %5ld\n", nje, nfeLS);
  637.  
  638. }
  639.  
  640. /*
  641.  * Check function return value.
  642.  *    opt == 0 means SUNDIALS function allocates memory so check if
  643.  *             returned NULL pointer
  644.  *    opt == 1 means SUNDIALS function returns an integer value so check if
  645.  *             retval < 0
  646.  *    opt == 2 means function allocates memory so check if returned
  647.  *             NULL pointer
  648.  */
  649.  
  650. static int check_retval(void *returnvalue, const char *funcname, int opt)
  651. {
  652.   int *retval;
  653.  
  654.   /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
  655.   if (opt == 0 && returnvalue == NULL) {
  656.     fprintf(stderr,
  657.             "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
  658.         funcname);
  659.     return(1); }
  660.  
  661.   /* Check if retval < 0 */
  662.   else if (opt == 1) {
  663.     retval = (int *) returnvalue;
  664.     if (*retval < 0) {
  665.       fprintf(stderr,
  666.               "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
  667.           funcname, *retval);
  668.       return(1); }}
  669.  
  670.   /* Check if function returned NULL pointer - no memory allocated */
  671.   else if (opt == 2 && returnvalue == NULL) {
  672.     fprintf(stderr,
  673.             "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
  674.         funcname);
  675.     return(1); }
  676.  
  677.   return(0);
  678. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top