Add conda recipe for sundials

sundials/build.sh

@@ -0,0 +1,4 @@
+#!/bin/sh
+mkdir build && cd build
+cmake -DCMAKE_INSTALL_PREFIX=$PREFIX -DCMAKE_BUILD_TYPE:STRING=Release -DBUILD_SHARED_LIBS:BOOL="1" -DBUILD_STATIC_LIBS:BOOL="0" -DEXAMPLES_ENABLE:BOOL="0" -DEXAMPLES_INSTALL:BOOL="0" -DLAPACK_ENABLE:BOOL="1" -DOPENMP_ENABLE:BOOL="0" ..
+make install

sundials/cvRoberts_dns.c

@@ -0,0 +1,360 @@
+/*
+ * -----------------------------------------------------------------
+ * $Revision: 4074 $
+ * $Date: 2014-04-23 14:13:52 -0700 (Wed, 23 Apr 2014) $
+ * -----------------------------------------------------------------
+ * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and
+ *                Radu Serban @ LLNL
+ * -----------------------------------------------------------------
+ * Example problem:
+ * 
+ * The following is a simple example problem, with the coding
+ * needed for its solution by CVODE. The problem is from
+ * chemical kinetics, and consists of the following three rate
+ * equations:         
+ *    dy1/dt = -.04*y1 + 1.e4*y2*y3
+ *    dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2
+ *    dy3/dt = 3.e7*(y2)^2
+ * on the interval from t = 0.0 to t = 4.e10, with initial
+ * conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.
+ * While integrating the system, we also use the rootfinding
+ * feature to find the points at which y1 = 1e-4 or at which
+ * y3 = 0.01. This program solves the problem with the BDF method,
+ * Newton iteration with the CVDENSE dense linear solver, and a
+ * user-supplied Jacobian routine.
+ * It uses a scalar relative tolerance and a vector absolute
+ * tolerance. Output is printed in decades from t = .4 to t = 4.e10.
+ * Run statistics (optional outputs) are printed at the end.
+ * -----------------------------------------------------------------
+ */
+
+#include <stdio.h>
+
+/* Header files with a description of contents used */
+
+#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts. */
+#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */
+#include <cvode/cvode_dense.h>       /* prototype for CVDense */
+#include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */
+#include <sundials/sundials_types.h> /* definition of type realtype */
+
+/* User-defined vector and matrix accessor macros: Ith, IJth */
+
+/* These macros are defined in order to write code which exactly matches
+   the mathematical problem description given above.
+
+   Ith(v,i) references the ith component of the vector v, where i is in
+   the range [1..NEQ] and NEQ is defined below. The Ith macro is defined
+   using the N_VIth macro in nvector.h. N_VIth numbers the components of
+   a vector starting from 0.
+
+   IJth(A,i,j) references the (i,j)th element of the dense matrix A, where
+   i and j are in the range [1..NEQ]. The IJth macro is defined using the
+   DENSE_ELEM macro in dense.h. DENSE_ELEM numbers rows and columns of a
+   dense matrix starting from 0. */
+
+#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 */
+
+
+/* Problem Constants */
+
+#define NEQ   3                /* number of equations  */
+#define Y1    RCONST(1.0)      /* initial y components */
+#define Y2    RCONST(0.0)
+#define Y3    RCONST(0.0)
+#define RTOL  RCONST(1.0e-4)   /* scalar relative tolerance            */
+#define ATOL1 RCONST(1.0e-8)   /* vector absolute tolerance components */
+#define ATOL2 RCONST(1.0e-14)
+#define ATOL3 RCONST(1.0e-6)
+#define T0    RCONST(0.0)      /* initial time           */
+#define T1    RCONST(0.4)      /* first output time      */
+#define TMULT RCONST(10.0)     /* output time factor     */
+#define NOUT  12               /* number of output times */
+
+
+/* Functions Called by the Solver */
+
+static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
+
+static int g(realtype t, N_Vector y, realtype *gout, void *user_data);
+
+static int Jac(long int N, realtype t,
+               N_Vector y, N_Vector fy, DlsMat J, void *user_data,
+               N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
+
+/* Private functions to output results */
+
+static void PrintOutput(realtype t, realtype y1, realtype y2, realtype y3);
+static void PrintRootInfo(int root_f1, int root_f2);
+
+/* Private function to print final statistics */
+
+static void PrintFinalStats(void *cvode_mem);
+
+/* Private function to check function return values */
+
+static int check_flag(void *flagvalue, char *funcname, int opt);
+
+
+/*
+ *-------------------------------
+ * Main Program
+ *-------------------------------
+ */
+
+int main()
+{
+  realtype reltol, t, tout;
+  N_Vector y, abstol;
+  void *cvode_mem;
+  int flag, flagr, iout;
+  int rootsfound[2];
+
+  y = abstol = NULL;
+  cvode_mem = NULL;
+
+  /* Create serial vector of length NEQ for I.C. and abstol */
+  y = N_VNew_Serial(NEQ);
+  if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);
+  abstol = N_VNew_Serial(NEQ); 
+  if (check_flag((void *)abstol, "N_VNew_Serial", 0)) return(1);
+
+  /* Initialize y */
+  Ith(y,1) = Y1;
+  Ith(y,2) = Y2;
+  Ith(y,3) = Y3;
+
+  /* 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 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 inital 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 CVodeSVtolerances to specify the scalar relative tolerance
+   * and vector absolute tolerances */
+  flag = CVodeSVtolerances(cvode_mem, reltol, abstol);
+  if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);
+
+  /* Call CVodeRootInit to specify the root function g with 2 components */
+  flag = CVodeRootInit(cvode_mem, 2, g);
+  if (check_flag(&flag, "CVodeRootInit", 1)) return(1);
+
+  /* Call CVDense to specify the CVDENSE dense linear solver */
+  flag = CVDense(cvode_mem, NEQ);
+  if (check_flag(&flag, "CVDense", 1)) return(1);
+
+  /* Set the Jacobian routine to Jac (user-supplied) */
+  flag = CVDlsSetDenseJacFn(cvode_mem, Jac);
+  if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);
+
+  /* In loop, call CVode, print results, and test for error.
+     Break out of loop when NOUT preset output times have been reached.  */
+  printf(" \n3-species kinetics problem\n\n");
+
+  iout = 0;  tout = T1;
+  while(1) {
+    flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
+    PrintOutput(t, Ith(y,1), Ith(y,2), Ith(y,3));
+
+    if (flag == CV_ROOT_RETURN) {
+      flagr = CVodeGetRootInfo(cvode_mem, rootsfound);
+      if (check_flag(&flagr, "CVodeGetRootInfo", 1)) return(1);
+      PrintRootInfo(rootsfound[0],rootsfound[1]);
+    }
+
+    if (check_flag(&flag, "CVode", 1)) break;
+    if (flag == CV_SUCCESS) {
+      iout++;
+      tout *= TMULT;
+    }
+
+    if (iout == NOUT) break;
+  }
+
+  /* Print some final statistics */
+  PrintFinalStats(cvode_mem);
+
+  /* Free y and abstol vectors */
+  N_VDestroy_Serial(y);
+  N_VDestroy_Serial(abstol);
+
+  /* Free integrator memory */
+  CVodeFree(&cvode_mem);
+
+  return(0);
+}
+
+
+/*
+ *-------------------------------
+ * Functions called by the solver
+ *-------------------------------
+ */
+
+/*
+ * f routine. Compute function f(t,y). 
+ */
+
+static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
+{
+  realtype y1, y2, y3, yd1, yd3;
+
+  y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
+
+  yd1 = Ith(ydot,1) = RCONST(-0.04)*y1 + RCONST(1.0e4)*y2*y3;
+  yd3 = Ith(ydot,3) = RCONST(3.0e7)*y2*y2;
+        Ith(ydot,2) = -yd1 - yd3;
+
+  return(0);
+}
+
+/*
+ * g routine. Compute functions g_i(t,y) for i = 0,1. 
+ */
+
+static int g(realtype t, N_Vector y, realtype *gout, void *user_data)
+{
+  realtype y1, y3;
+
+  y1 = Ith(y,1); y3 = Ith(y,3);
+  gout[0] = y1 - RCONST(0.0001);
+  gout[1] = y3 - RCONST(0.01);
+
+  return(0);
+}
+
+/*
+ * Jacobian routine. Compute J(t,y) = df/dy. *
+ */
+
+static int Jac(long int N, realtype t,
+               N_Vector y, N_Vector fy, DlsMat J, void *user_data,
+               N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
+{
+  realtype y1, y2, y3;
+
+  y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
+
+  IJth(J,1,1) = RCONST(-0.04);
+  IJth(J,1,2) = RCONST(1.0e4)*y3;
+  IJth(J,1,3) = RCONST(1.0e4)*y2;
+  IJth(J,2,1) = RCONST(0.04); 
+  IJth(J,2,2) = RCONST(-1.0e4)*y3-RCONST(6.0e7)*y2;
+  IJth(J,2,3) = RCONST(-1.0e4)*y2;
+  IJth(J,3,2) = RCONST(6.0e7)*y2;
+
+  return(0);
+}
+
+/*
+ *-------------------------------
+ * Private helper functions
+ *-------------------------------
+ */
+
+static void PrintOutput(realtype t, realtype y1, realtype y2, realtype y3)
+{
+#if defined(SUNDIALS_EXTENDED_PRECISION)
+  printf("At t = %0.4Le      y =%14.6Le  %14.6Le  %14.6Le\n", t, y1, y2, y3);
+#elif defined(SUNDIALS_DOUBLE_PRECISION)
+  printf("At t = %0.4e      y =%14.6e  %14.6e  %14.6e\n", t, y1, y2, y3);
+#else
+  printf("At t = %0.4e      y =%14.6e  %14.6e  %14.6e\n", t, y1, y2, y3);
+#endif
+
+  return;
+}
+
+static void PrintRootInfo(int root_f1, int root_f2)
+{
+  printf("    rootsfound[] = %3d %3d\n", root_f1, root_f2);
+
+  return;
+}
+
+/* 
+ * 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);
+}

sundials/cvRoberts_dns.out

@@ -0,0 +1,24 @@
+
+3-species kinetics problem
+
+At t = 2.6391e-01      y =  9.899653e-01    3.470564e-05    1.000000e-02
+    rootsfound[] =   0   1
+At t = 4.0000e-01      y =  9.851641e-01    3.386242e-05    1.480205e-02
+At t = 4.0000e+00      y =  9.055097e-01    2.240338e-05    9.446793e-02
+At t = 4.0000e+01      y =  7.158009e-01    9.185098e-06    2.841900e-01
+At t = 4.0000e+02      y =  4.505440e-01    3.223217e-06    5.494528e-01
+At t = 4.0000e+03      y =  1.831964e-01    8.942051e-07    8.168027e-01
+At t = 4.0000e+04      y =  3.898104e-02    1.621656e-07    9.610188e-01
+At t = 4.0000e+05      y =  4.938672e-03    1.985172e-08    9.950613e-01
+At t = 4.0000e+06      y =  5.166093e-04    2.067499e-09    9.994834e-01
+At t = 2.0800e+07      y =  1.000000e-04    4.000395e-10    9.999000e-01
+    rootsfound[] =  -1   0
+At t = 4.0000e+07      y =  5.206409e-05    2.082671e-10    9.999479e-01
+At t = 4.0000e+08      y =  5.211241e-06    2.084507e-11    9.999948e-01
+At t = 4.0000e+09      y =  5.200520e-07    2.080209e-12    9.999995e-01
+At t = 4.0000e+10      y =  5.699485e-08    2.279794e-13    9.999999e-01
+
+Final Statistics:
+nst = 579    nfe  = 817    nsetups = 118    nfeLS = 0      nje = 12
+nni = 813    ncfn = 0      netf = 31     nge = 615
+