Merge pull request #26 from seamplex/fix-segfault-in-outward-normal-w…

…ith-O3

Support for sundials > 5

configure.ac

@@ -137,6 +137,12 @@ AS_IF([test "x${with_sundials}" != "xno"] , [
       AC_MSG_FAILURE([--with-sundials was given but test for libsundials-dev library failed])
     )
   )  
+  AC_CHECK_LIB([sundials_core], [SUNContext_Create],
+    [],
+    AS_IF([test "x${with_sundials}" != "xcheck"],
+      AC_MSG_FAILURE([--with-sundials was given but test for libsundials-dev library failed])
+    )
+  )  
   AC_CHECK_HEADER([nvector/nvector_serial.h],
     [],
     AS_IF([test "x${with_sundials}" != "xcheck"],

src/feenox.h

@@ -72,12 +72,16 @@
 
 #ifdef HAVE_SUNDIALS
  #include <ida/ida.h>
- #include <ida/ida_direct.h>
+ // #include <ida/ida_direct.h>
  #include <nvector/nvector_serial.h>
  #include <sundials/sundials_types.h>
  #include <sundials/sundials_math.h>
  #include <sunmatrix/sunmatrix_dense.h>
  #include <sunlinsol/sunlinsol_dense.h>
+
+ #ifndef sunrealtype 
+  #define sunrealty realtype
+ #endif
 #endif
 
 // petsc complains if its headers are included from c++ within extern C
@@ -2094,7 +2098,7 @@ extern int feenox_add_dae(const char *lhs, const char *rhs);
 extern int feenox_dae_init(void);
 extern int feenox_dae_ic(void);
 #ifdef HAVE_SUNDIALS
-extern int feenox_ida_dae(realtype t, N_Vector yy, N_Vector yp, N_Vector rr, void *params);
+extern int feenox_ida_dae(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, void *params);
 #else
 extern int feenox_ida_dae(void);
 #endif

src/math/dae.c

@@ -198,7 +198,7 @@ char *feenox_find_first_dot(const char *s) {
 int feenox_dae_init(void) {
 
 #ifdef HAVE_SUNDIALS
-  if (sizeof(realtype) != sizeof(double)) {
+  if (sizeof(sunrealtype) != sizeof(double)) {
     feenox_push_error_message("SUNDIALS's realtype is not double");
     return FEENOX_ERROR;
   }
@@ -480,7 +480,7 @@ int feenox_dae_ic(void) {
 
 
 #ifdef HAVE_SUNDIALS
-int feenox_ida_dae(realtype t, N_Vector yy, N_Vector yp, N_Vector rr, void *params) {
+int feenox_ida_dae(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, void *params) {
 
   // this function is called many times in a single FeenoX time step for intermediate
   // times which are internal to IDA, nevertheless there might be some residuals

src/parser/parser.c

@@ -2123,7 +2123,7 @@ int feenox_parse_phase_space(void) {
     return FEENOX_ERROR;
   }
 #ifdef HAVE_SUNDIALS
-  if (sizeof(realtype) != sizeof(double)) {
+  if (sizeof(sunrealtype) != sizeof(double)) {
     feenox_push_error_message("\nSUNDIALS was compiled using a different word size than FeenoX, please recompile with double precision floating point arithmetic.");
     return FEENOX_ERROR;
   }

tests/lag.sh

@@ -20,7 +20,7 @@ exitifwrong $?
 # check TIME_PATH
 echo -n "time_path_raw.fee ... "
 
-if [ "x$(${feenox} ${dir}/time_path_sundials.fee | wc -l)" != "x5" ]; then
+if [ "x$(${feenox} ${dir}/time_path_raw.fee | wc -l)" != "x5" ]; then
   echo "failed"
   return 1
 fi

tests/time_path_sundials.fee

@@ -1,6 +1,6 @@
 PHASE_SPACE x
 TIME_PATH 1e-2 5e-2 1e-1 5e-1 1
-end_time = 1
+end_time = 1+0.1
 x_0 = 1
 x_dot = -x
 IF in_time_path