diff options
Diffstat (limited to 'packages/sundials/src/helpers.c')
-rw-r--r-- | packages/sundials/src/helpers.c | 36 |
1 files changed, 35 insertions, 1 deletions
diff --git a/packages/sundials/src/helpers.c b/packages/sundials/src/helpers.c index 3498119..3e1a9a8 100644 --- a/packages/sundials/src/helpers.c +++ b/packages/sundials/src/helpers.c | |||
@@ -1,5 +1,12 @@ | |||
1 | #include <stdio.h> | 1 | #include <stdio.h> |
2 | #include "helpers.h" | 2 | #include <math.h> |
3 | #include <arkode/arkode.h> /* prototypes for ARKODE fcts., consts. */ | ||
4 | #include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */ | ||
5 | #include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */ | ||
6 | #include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */ | ||
7 | #include <arkode/arkode_direct.h> /* access to ARKDls interface */ | ||
8 | #include <sundials/sundials_types.h> /* definition of type realtype */ | ||
9 | #include <sundials/sundials_math.h> | ||
3 | 10 | ||
4 | /* Check function return value... | 11 | /* Check function return value... |
5 | opt == 0 means SUNDIALS function allocates memory so check if | 12 | opt == 0 means SUNDIALS function allocates memory so check if |
@@ -35,3 +42,30 @@ int check_flag(void *flagvalue, const char *funcname, int opt) | |||
35 | 42 | ||
36 | return 0; | 43 | return 0; |
37 | } | 44 | } |
45 | |||
46 | /* f routine to compute the ODE RHS function f(t,y). */ | ||
47 | int f(realtype t, N_Vector y, N_Vector ydot, void *user_data) | ||
48 | { | ||
49 | realtype *rdata = (realtype *) user_data; /* cast user_data to realtype */ | ||
50 | realtype lamda = rdata[0]; /* set shortcut for stiffness parameter */ | ||
51 | realtype u = NV_Ith_S(y,0); /* access current solution value */ | ||
52 | |||
53 | /* fill in the RHS function: "NV_Ith_S" accesses the 0th entry of ydot */ | ||
54 | NV_Ith_S(ydot,0) = lamda*u + 1.0/(1.0+t*t) - lamda*atan(t); | ||
55 | |||
56 | return 0; /* return with success */ | ||
57 | } | ||
58 | |||
59 | /* Jacobian routine to compute J(t,y) = df/dy. */ | ||
60 | int Jac(realtype t, N_Vector y, N_Vector fy, SUNMatrix J, | ||
61 | void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) | ||
62 | { | ||
63 | realtype *rdata = (realtype *) user_data; /* cast user_data to realtype */ | ||
64 | realtype lamda = rdata[0]; /* set shortcut for stiffness parameter */ | ||
65 | realtype *Jdata = SUNDenseMatrix_Data(J); | ||
66 | |||
67 | /* Fill in Jacobian of f: set the first entry of the data array to set the (0,0) entry */ | ||
68 | Jdata[0] = lamda; | ||
69 | |||
70 | return 0; /* return with success */ | ||
71 | } | ||