ne_sim_cvode.c

/*===========================================================================
 NetEvo Foundation Library
 Copyright (C) 2009, 2010 Thomas E. Gorochowski <tgorochowski@me.com>
 Bristol Centre for Complexity Sciences, University of Bristol, Bristol, UK
 ---------------------------------------------------------------------------- 
 NetEvo is a computing framework designed to allow researchers to investigate 
 evolutionary aspects of dynamical complex networks. By providing tools to 
 easily integrate each of these factors in a coherent way, it is hoped a 
 greater understanding can be gained of key attributes and features displayed 
 by complex systems.
 
 NetEvo is open-source software released under the Open Source Initiative 
 (OSI) approved Non-Profit Open Software License ("Non-Profit OSL") 3.0. 
 Detailed information about this licence can be found in the COPYING file 
 included as part of the source distribution.
 
 This library is distributed in the hope that it will be useful, but
 WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 ============================================================================*/
/*
 * EXPERIMENTAL... DO NOT USE.
 *
 * Simulator of dynamics for ODE models that uses the CVODE (Sundials) suite of
 * solvers. This is suitable for both non-stiff and stiff problems, however,
 * type should be specified in the parameters when allocating the solver.
 */


#include "ne_sim_cvode.h"
#include <string.h>


#define Ith(v,i)    NV_Ith_S(v,i-1)       /* Ith numbers components 1..NEQ */
#define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1) /* IJth numbers rows,cols 1..NEQ */

ne_vec_dyn_t * ne_sim_cvode_run (ne_system_t *system, void *params, double length, double *iy, 
                                 FILE *file)
{
   /* TODO: Handle errors better during initialisation step */
   int dimension, flag, i;
   double *yCopy;
   realtype cvRelTol, cvAbsTol, t, tout;
   N_Vector y;
   void *cvode_mem;
   ne_vec_dyn_t *vecOut;
   ne_sim_cvode_params_t *cvodeParams;
   
   cvodeParams = (ne_sim_cvode_params_t *)params;
   y = NULL;
   cvode_mem = NULL;
   dimension = system->nodes * system->nodeStates;
   vecOut = ne_vec_dyn_alloc(dimension + 1);

   /* Create serial vector of length NEQ for I.C. and abstol */
   y = N_VNew_Serial(dimension);
   if (check_flag((void *)y, "N_VNew_Serial", 0)) return NULL;
   
   /* Set initial values */
   for (i=0; i<dimension; i++) {
      NV_Ith_S(y,i) = iy[i];
   }
      
   /* Call CVodeCreate to create the solver memory */
   if (cvodeParams->stiff == TRUE) {
      cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
   }
   else {
      cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
   }
   if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return NULL;
   
   /* 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. */
   flag = CVodeInit(cvode_mem, ne_sim_cvode_fn, 0.0, y);
   if (check_flag(&flag, "CVodeInit", 1)) return NULL;
   
   /* TODO: Allow for absTol to be specified as an array */
   /* Call CVodeSVtolerances to specify the scalar relative tolerance
    * and vector absolute tolerances */
   cvRelTol = (realtype)cvodeParams->relTol;
   cvAbsTol = (realtype)cvodeParams->absTol;
   flag = CVodeSStolerances(cvode_mem, cvRelTol, cvAbsTol);
   if (check_flag(&flag, "CVodeSStolerances", 1)) return NULL;
   
   /* Call CVDense to specify the CVDENSE dense linear solver */
   flag = CVDense(cvode_mem, dimension);
   if (check_flag(&flag, "CVDense", 1)) return NULL;
   
   tout = 0.0;
   do {
      if (flag == CV_SUCCESS) {
         tout += cvodeParams->step;
         if (tout > length) {
            tout = length;
         }
      }
      
      /* Run the solver to the next time step */
      flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
      
      /* Make a copy of the results and append to the output vector */
      yCopy = (double*)malloc((sizeof(double)*dimension)+1);
      yCopy[0] = t;
      for (i=0; i<dimension; i++) {
         yCopy[i+1] = NV_Ith_S(y,i);
      }
      ne_vec_dyn_append_item(vecOut, yCopy);
      
      if (check_flag(&flag, "CVode", 1)) break;
      
   } while (tout < length);
   
   /* Free y */
   N_VDestroy_Serial(y);
   
   /* Free integrator memory */
   CVodeFree(&cvode_mem);
   
   return vecOut;
}


int ne_sim_cvode_fn (realtype t, N_Vector y, N_Vector ydot, void *params)
{
   /* Function to calculate: ydot_i = F(x_i) + sum_{j \in E} G(x_i, x_j) */
   
   ne_int_t states, i, j, k;
   ne_system_t *system;
   double *ydotNew, *ydotCum;
   
   system = (ne_system_t *)params;
   states = ne_system_node_states(system);
   
   ydotNew = (double *)malloc(sizeof(double) * states);
   ydotCum = (double *)malloc(sizeof(double) * states);
   
   /* Calculate the node dynamics F(x) */
   for (i=0; i<system->nodes; i++) {
      (ne_system_dyn_node(system, i)->fn)(i, system, (double)t, 
                                          (double *)&NV_Ith_S(y, i*states), 
                                          (double *)&NV_Ith_S(ydot, i*states), 
                                          ne_system_dyn_node_params(system, i));
   }
   
   /* Calculate the edge dynamics G(x) */
   if (ne_system_node_param(system, i, NE_NODE_PARAM_EDGE_INT_TYPE) == NE_INT_SUM) {
      
      /* Sum along edges (NE_INT_SUM) */
      for (i=0; i<system->nodes; i++) {
         memset(ydotCum, 0, sizeof(double)*states);
         for (j=0; j<system->nodes; j++) {
            if (ne_system_edge_exists(system, i, j) == TRUE) {
               memset(ydotNew, 0, sizeof(double)*states);
               (ne_system_dyn_edge(system, i, j)->fn)(i, j, system, (double)t, 
                                                      (double *)NV_DATA_S(y), 
                                                      ydotNew, 
                                                      ne_system_dyn_edge_params(system, i, j));
               /* Keep track of the cumulative value */
               for (k=0; k<states; k++) {
                  ydotCum[k] += ydotNew[k];
               }
            }
         }
         /* Update the actual node values */
         for (k=0; k<states; k++) {
            NV_Ith_S(ydot, (i*states)+k) += ydotCum[k];
         }
      }
   }
   else if (ne_system_node_param(system, i, NE_NODE_PARAM_EDGE_INT_TYPE) == NE_INT_MUT) {
      
      /* Sum along edges (NE_INT_MUT) */
      for (i=0; i<system->nodes; i++) {
         memset(ydotCum, 1, sizeof(double)*states);
         for (j=0; j<system->nodes; j++) {
            if (ne_system_edge_exists(system, i, j) == TRUE) {
               memset(ydotNew, 0, sizeof(double)*states);
               (ne_system_dyn_edge(system, i, j)->fn)(i, j, system, (double)t, 
                                                      (double *)NV_DATA_S(y), 
                                                      ydotNew, 
                                                      ne_system_dyn_edge_params(system, i, j));
               /* Keep track of the cumulative value */
               for (k=0; k<states; k++) {
                  ydotCum[k] *= ydotNew[k];
               }
            }
         }
         /* Update the actual node values */
         for (k=0; k<states; k++) {
            NV_Ith_S(ydot, (i*states)+k) += ydotCum[k];
         }
      }
   }
   
   free(ydotNew);
   free(ydotCum);
   
   return 0;
}



/* 
 * Get and print some final statistics
 */
/*
static void PrintFinalStats(void *cvode_mem)
{
   long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
   int flag;
   
   flag = CVodeGetNumSteps(cvode_mem, &nst);
   check_flag(&flag, "CVodeGetNumSteps", 1);
   flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);
   check_flag(&flag, "CVodeGetNumRhsEvals", 1);
   flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
   check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);
   flag = CVodeGetNumErrTestFails(cvode_mem, &netf);
   check_flag(&flag, "CVodeGetNumErrTestFails", 1);
   flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
   check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);
   flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
   check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);
   
   flag = CVDlsGetNumJacEvals(cvode_mem, &nje);
   check_flag(&flag, "CVDlsGetNumJacEvals", 1);
   flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);
   check_flag(&flag, "CVDlsGetNumRhsEvals", 1);
   
   flag = CVodeGetNumGEvals(cvode_mem, &nge);
   check_flag(&flag, "CVodeGetNumGEvals", 1);
   
   printf("\nFinal Statistics:\n");
   printf("nst = %-6ld nfe  = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld\n",
         nst, nfe, nsetups, nfeLS, nje);
   printf("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n \n",
         nni, ncfn, netf, nge);
}
*/


/*
 * 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, 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);
}