Merge pull request #2 from idontgetoutmuch/revert-1-sundials

Revert "Cleanups to Sundials PR"

2 files changed, 34 insertions, 31 deletions

diff --git a/packages/sundials/src/Main.hs b/packages/sundials/src/Main.hs
index d1f6755..5e51372 100644
--- a/packages/sundials/src/Main.hs
+++ b/packages/sundials/src/Main.hs
@@ -1,5 +1,6 @@
-module Main where
+{-# OPTIONS_GHC -Wall #-}
 
+import qualified Data.Vector.Storable as V
 import           Numeric.Sundials.Arkode.ODE
 import           Numeric.LinearAlgebra
 
@@ -8,10 +9,11 @@ import qualified Diagrams.Prelude as D
 import           Diagrams.Backend.Rasterific
 
 import           Control.Lens
+import           Data.List (zip4)
 
-brusselator :: a -> [Double] -> [Double]
-brusselator _t x = [ a - (w + 1) * u + v * u^(2::Int)
-                   , w * u - v * u^(2::Int)
+
+brusselator _t x = [ a - (w + 1) * u + v * u^2
+                   , w * u - v * u^2
                    , (b - w) / eps - w * u
                    ]
   where
@@ -22,7 +24,6 @@ brusselator _t x = [ a - (w + 1) * u + v * u^(2::Int)
     v = x !! 1
     w = x !! 2
 
-stiffish :: Double -> [Double] -> [Double]
 stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
   where
     lamda = -100.0
@@ -44,14 +45,15 @@ kSaxis xs = P.r2Axis &~ do
 
 main :: IO ()
 main = do
-  do let res = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
-     putStrLn $ show res
-     renderRasterific "diagrams/brusselator.png"
-                      (D.dims2D 500.0 500.0)
-                      (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res))
-
-  do let res = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
-     putStrLn $ show res
-     renderRasterific "diagrams/stiffish.png"
-                      (D.dims2D 500.0 500.0)
-                      (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res))
+  let res = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
+  putStrLn $ show res
+  renderRasterific "diagrams/brusselator.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res))
+
+  let res = odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
+  putStrLn $ show res
+  renderRasterific "diagrams/stiffish.png"
+                   (D.dims2D 500.0 500.0)
+                   (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res))
+  
diff --git a/packages/sundials/src/Numeric/Sundials/Arkode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index 6d9a1b2..f432951 100644
--- a/packages/sundials/src/Numeric/Sundials/Arkode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -1,14 +1,14 @@
+{-# OPTIONS_GHC -Wall #-}
+
+{-# LANGUAGE QuasiQuotes #-}
+{-# LANGUAGE TemplateHaskell #-}
 {-# LANGUAGE MultiWayIf #-}
 {-# LANGUAGE OverloadedStrings #-}
-{-# LANGUAGE QuasiQuotes #-}
 {-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TemplateHaskell #-}
 
-module Numeric.Sundials.Arkode.ODE
-  ( SundialsDiagnostics(..)
-  , solveOde
-  , odeSolve
-  ) where
+module Numeric.Sundials.Arkode.ODE ( solveOde
+                                   , odeSolve
+                                   ) where
 
 import qualified Language.C.Inline as C
 import qualified Language.C.Inline.Unsafe as CU
@@ -40,10 +40,10 @@ C.include "<stdio.h>"
 C.include "<math.h>"
 C.include "<arkode/arkode.h>"                 -- prototypes for ARKODE fcts., consts.
 C.include "<nvector/nvector_serial.h>"        -- serial N_Vector types, fcts., macros
-C.include "<sunmatrix/sunmatrix_dense.h>"     -- access to dense SUNMatrix
-C.include "<sunlinsol/sunlinsol_dense.h>"     -- access to dense SUNLinearSolver
-C.include "<arkode/arkode_direct.h>"          -- access to ARKDls interface
-C.include "<sundials/sundials_types.h>"       -- definition of type realtype
+C.include "<sunmatrix/sunmatrix_dense.h>"     -- access to dense SUNMatrix           
+C.include "<sunlinsol/sunlinsol_dense.h>"     -- access to dense SUNLinearSolver     
+C.include "<arkode/arkode_direct.h>"          -- access to ARKDls interface          
+C.include "<sundials/sundials_types.h>"       -- definition of type realtype         
 C.include "<sundials/sundials_math.h>"
 C.include "../../../helpers.h"
 
@@ -108,7 +108,7 @@ odeSolve f y0 ts = case solveOde g (V.fromList y0) (V.fromList $ toList ts) of
     nC = length y0
     g t x0 = V.fromList $ f t (V.toList x0)
 
-solveOde ::
+solveOde :: 
   (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
           -> V.Vector Double                     -- ^ Initial conditions
           -> V.Vector Double                     -- ^ Desired solution times
@@ -116,7 +116,7 @@ solveOde ::
 solveOde f y0 tt = case solveOdeC (coerce f) (coerce y0) (coerce tt) of
                      Left c -> Left $ fromIntegral c
                      Right (v, d) -> Right (coerce v, d)
-
+                        
 solveOdeC ::
   (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
           -> V.Vector CDouble -- ^ Initial conditions
@@ -285,13 +285,13 @@ solveOdeC fun f0 ts = unsafePerformIO $ do
                          flag = ARKDlsGetNumRhsEvals(arkode_mem, &nfeLS);
                          check_flag(&flag, "ARKDlsGetNumRhsEvals", 1);
                          ($vec-ptr:(long int *diagMut))[9] = ncfn;
-
+                         
                          /* Clean up and return */
                          N_VDestroy(y);            /* Free y vector */
                          ARKodeFree(&arkode_mem);  /* Free integrator memory */
                          SUNLinSolFree(LS);        /* Free linear solver */
                          SUNMatDestroy(A);         /* Free A matrix */
-
+ 
                          return flag;
                        } |]
   if res == 0
@@ -311,3 +311,4 @@ solveOdeC fun f0 ts = unsafePerformIO $ do
       return $ Right (m, d)
     else do
       return $ Left res
+