summaryrefslogtreecommitdiff
path: root/packages/sundials/src/Numeric
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2018-03-29 08:29:56 +0100
committerDominic Steinitz <dominic@steinitz.org>2018-03-29 08:29:56 +0100
commit72d9adbfb0ee5593281aa9e443c6e703ef10f36d (patch)
treea82731f89fb85244552cca0386221a97dddfe841 /packages/sundials/src/Numeric
parent6148374bbf11e443ed3c64eb4add3f20d612f362 (diff)
Make tolerances parameters
Diffstat (limited to 'packages/sundials/src/Numeric')
-rw-r--r--packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs23
1 files changed, 10 insertions, 13 deletions
diff --git a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
index d0c58dc..2ede918 100644
--- a/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
+++ b/packages/sundials/src/Numeric/Sundials/ARKode/ODE.hs
@@ -151,7 +151,7 @@ odeSolve :: ODEMethod
151 -> Vector Double -- ^ desired solution times 151 -> Vector Double -- ^ desired solution times
152 -> Matrix Double -- ^ solution 152 -> Matrix Double -- ^ solution
153odeSolve method f y0 ts = 153odeSolve method f y0 ts =
154 case solveOde method g (V.fromList y0) (V.fromList $ toList ts) of 154 case solveOde method 1.0e-6 1.0e-10 g (V.fromList y0) (V.fromList $ toList ts) of
155 Left c -> error $ show c -- FIXME 155 Left c -> error $ show c -- FIXME
156 Right (v, _) -> (nR >< nC) (V.toList v) 156 Right (v, _) -> (nR >< nC) (V.toList v)
157 where 157 where
@@ -162,22 +162,27 @@ odeSolve method f y0 ts =
162 162
163solveOde :: 163solveOde ::
164 ODEMethod 164 ODEMethod
165 -> Double
166 -> Double
165 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 167 -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
166 -> V.Vector Double -- ^ Initial conditions 168 -> V.Vector Double -- ^ Initial conditions
167 -> V.Vector Double -- ^ Desired solution times 169 -> V.Vector Double -- ^ Desired solution times
168 -> Either Int ((V.Vector Double), SundialsDiagnostics) -- ^ Error code or solution 170 -> Either Int ((V.Vector Double), SundialsDiagnostics) -- ^ Error code or solution
169solveOde method f y0 tt = 171solveOde method relTol absTol f y0 tt =
170 case solveOdeC (fromIntegral $ fromEnum method) (coerce f) (coerce y0) (coerce tt) of 172 case solveOdeC (fromIntegral $ fromEnum method) (CDouble relTol) (CDouble absTol)
173 (coerce f) (coerce y0) (coerce tt) of
171 Left c -> Left $ fromIntegral c 174 Left c -> Left $ fromIntegral c
172 Right (v, d) -> Right (coerce v, d) 175 Right (v, d) -> Right (coerce v, d)
173 176
174solveOdeC :: 177solveOdeC ::
175 CInt -> 178 CInt ->
179 CDouble ->
180 CDouble ->
176 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\) 181 (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
177 -> V.Vector CDouble -- ^ Initial conditions 182 -> V.Vector CDouble -- ^ Initial conditions
178 -> V.Vector CDouble -- ^ Desired solution times 183 -> V.Vector CDouble -- ^ Desired solution times
179 -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution 184 -> Either CInt ((V.Vector CDouble), SundialsDiagnostics) -- ^ Error code or solution
180solveOdeC method fun f0 ts = unsafePerformIO $ do 185solveOdeC method relTol absTol fun f0 ts = unsafePerformIO $ do
181 let dim = V.length f0 186 let dim = V.length f0
182 nEq :: CLong 187 nEq :: CLong
183 nEq = fromIntegral dim 188 nEq = fromIntegral dim
@@ -216,15 +221,7 @@ solveOdeC method fun f0 ts = unsafePerformIO $ do
216 221
217 /* general problem parameters */ 222 /* general problem parameters */
218 realtype T0 = RCONST(($vec-ptr:(double *tMut))[0]); /* initial time */ 223 realtype T0 = RCONST(($vec-ptr:(double *tMut))[0]); /* initial time */
219
220 sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars. */ 224 sunindextype NEQ = $(sunindextype nEq); /* number of dependent vars. */
221 realtype reltol = 1.0e-6; /* tolerances */
222 realtype abstol = 1.0e-10;
223
224 /* Initial diagnostics output */
225 printf("\nAnalytical ODE test problem:\n");
226 printf(" reltol = %.1"ESYM"\n", reltol);
227 printf(" abstol = %.1"ESYM"\n\n",abstol);
228 225
229 /* Initialize data structures */ 226 /* Initialize data structures */
230 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */ 227 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
@@ -247,7 +244,7 @@ solveOdeC method fun f0 ts = unsafePerformIO $ do
247 if (check_flag(&flag, "ARKodeInit", 1)) return 1; 244 if (check_flag(&flag, "ARKodeInit", 1)) return 1;
248 245
249 /* Set routines */ 246 /* Set routines */
250 flag = ARKodeSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */ 247 flag = ARKodeSStolerances(arkode_mem, $(double relTol), $(double absTol));
251 if (check_flag(&flag, "ARKodeSStolerances", 1)) return 1; 248 if (check_flag(&flag, "ARKodeSStolerances", 1)) return 1;
252 249
253 /* Initialize dense matrix data structure and solver */ 250 /* Initialize dense matrix data structure and solver */