1 #include <stdio.h>
2
3 #include <cpodes/cpodes.h>
4 #include <cpodes/cpodes_lapack.h>
5 #include <nvector/nvector_serial.h>
6
7 #define Ith(v,i) NV_Ith_S(v,i-1)
8 #define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1)
9
10 /* Problem Constants */
11
12 #define ODE CP_EXPL
13
14 #define NEQ 3
15 #define Y1 RCONST(1.0)
16 #define Y2 RCONST(0.0)
17 #define Y3 RCONST(0.0)
18
19 #define RTOL RCONST(1.0e-4) /* 1e-4 */
20 #define ATOL1 RCONST(1.0e-8) /* 1e-8 */
21 #define ATOL2 RCONST(1.0e-14) /* 1e-14 */
22 #define ATOL3 RCONST(1.0e-6) /* 1e-6 */
23
24 #define T0 RCONST(0.0)
25 #define T1 RCONST(0.4)
26 #define TMULT RCONST(10.0)
27 #define NOUT 12
28
29 #define ONE RCONST(1.0)
30 #define ZERO RCONST(0.0)
31
32 /* Functions Called by the Solver */
33 static int res(realtype t, N_Vector y, N_Vector yp, N_Vector res, void *f_data);
34 static int f(realtype t, N_Vector y, N_Vector ydot, void *f_data);
35 static int g(realtype t, N_Vector y, N_Vector yp, realtype *gout, void *g_data);
36 static int jacI(int N, realtype t, realtype gm,
37 N_Vector y, N_Vector yp, N_Vector r,
38 DlsMat J, void *jac_data,
39 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
40 static int jacE(int N, realtype t,
41 N_Vector y, N_Vector yp,
42 DlsMat J, void *jac_data,
43 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
44
45
46 static void PrintFinalStats(void *cpode_mem);
47
48 /*
49 *-------------------------------
50 * Main Program
51 *-------------------------------
52 */
53
main()54 int main()
55 {
56 void *fct, *jac;
57 void *cpode_mem;
58 N_Vector yy, yp, abstol;
59 realtype reltol, t, tout;
60 int flag, flagr, iout;
61 int rootsfound[2];
62
63 /* Create serial vectors of length NEQ for I.C. and abstol */
64 yy = N_VNew_Serial(NEQ);
65 yp = N_VNew_Serial(NEQ);
66 abstol = N_VNew_Serial(NEQ);
67
68 /* Initialize y */
69 Ith(yy,1) = Y1;
70 Ith(yy,2) = Y2;
71 Ith(yy,3) = Y3;
72
73 /* Set tolerances */
74 reltol = RTOL;
75 Ith(abstol,1) = ATOL1;
76 Ith(abstol,2) = ATOL2;
77 Ith(abstol,3) = ATOL3;
78
79
80 if (ODE == CP_EXPL) {
81 fct = (void *)f;
82 jac = (void *)jacE;
83 } else {
84 f(T0, yy, yp, NULL);
85 fct = (void *)res;
86 jac = (void *)jacI;
87 }
88
89 /* Initialize solver */
90 cpode_mem = CPodeCreate(ODE, CP_BDF, CP_NEWTON);
91 flag = CPodeInit(cpode_mem, fct, NULL, T0, yy, yp, CP_SV, reltol, abstol);
92
93 /* Set initial step size */
94 /*
95 {
96 realtype h0;
97 h0 = 8.5e-14;
98 flag = CPodeSetInitStep(cpode_mem, h0);
99 }
100 */
101
102 /* Call CPodeRootInit to specify the root function g with 2 components */
103 flag = CPodeRootInit(cpode_mem, 2, g, NULL);
104
105 /* Call CPLapackDense to specify the CPLAPACK dense linear solver */
106 flag = CPLapackDense(cpode_mem, NEQ);
107
108 /* Set the Jacobian routine to jac (comment out for internal DQ) */
109 flag = CPDlsSetJacFn(cpode_mem, jac, NULL);
110
111 /* In loop, call CPode, print results, and test for error.
112 Break out of loop when NOUT preset output times have been reached. */
113 iout = 0; tout = T1;
114 while(1) {
115 flag = CPode(cpode_mem, tout, &t, yy, yp, CP_NORMAL);
116 printf("At t = %0.4le y =%14.6le %14.6le %14.6le\n",
117 t, Ith(yy,1), Ith(yy,2), Ith(yy,3));
118
119 if (flag < 0) break;
120
121 if (flag == CP_ROOT_RETURN) {
122 flagr = CPodeGetRootInfo(cpode_mem, rootsfound);
123 printf(" rootsfound[] = %3d %3d\n", rootsfound[0],rootsfound[1]);
124 }
125
126 if (flag == CP_SUCCESS) {
127 iout++;
128 tout *= TMULT;
129 }
130
131 if (iout == NOUT) break;
132 }
133
134 /* Print some final statistics */
135 PrintFinalStats(cpode_mem);
136
137 /* Free memory */
138 N_VDestroy_Serial(yy);
139 N_VDestroy_Serial(yp);
140 N_VDestroy_Serial(abstol);
141 CPodeFree(&cpode_mem);
142
143 return(0);
144 }
145
146
147 /*
148 *-------------------------------
149 * Functions called by the solver
150 *-------------------------------
151 */
152
153 /*
154 * residual function: F(y',y) = y' - f(y)
155 */
156
res(realtype t,N_Vector y,N_Vector yp,N_Vector res,void * f_data)157 static int res(realtype t, N_Vector y, N_Vector yp, N_Vector res, void *f_data)
158 {
159 realtype y1, y2, y3, yp1, yp2, yp3;
160
161 y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
162 yp1 = Ith(yp,1); yp2 = Ith(yp,2); yp3 = Ith(yp,3);
163
164 Ith(res,1) = yp1 - (RCONST(-0.04)*y1 + RCONST(1.0e4)*y2*y3);
165 Ith(res,2) = yp2 + (RCONST(-0.04)*y1 + RCONST(1.0e4)*y2*y3 + RCONST(3.0e7)*y2*y2);
166 Ith(res,3) = yp3 - RCONST(3.0e7)*y2*y2;
167
168 return(0);
169 }
170
171 /*
172 * f routine. Compute function f(t,y).
173 */
174
f(realtype t,N_Vector y,N_Vector ydot,void * f_data)175 static int f(realtype t, N_Vector y, N_Vector ydot, void *f_data)
176 {
177 realtype y1, y2, y3, yd1, yd3;
178
179 y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
180
181 yd1 = Ith(ydot,1) = RCONST(-0.04)*y1 + RCONST(1.0e4)*y2*y3;
182 yd3 = Ith(ydot,3) = RCONST(3.0e7)*y2*y2;
183 Ith(ydot,2) = -yd1 - yd3;
184
185 return(0);
186 }
187
188 /*
189 * Jacobian routine. Compute J(t,y) = dF/dy' + gm * dF/dy = I - gm*df/dy.
190 */
jacI(int N,realtype t,realtype gm,N_Vector y,N_Vector yp,N_Vector r,DlsMat J,void * jac_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)191 static int jacI(int N, realtype t, realtype gm,
192 N_Vector y, N_Vector yp, N_Vector r,
193 DlsMat J, void *jac_data,
194 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
195 {
196
197 realtype y1, y2, y3;
198
199 y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
200
201 IJth(J,1,1) = ONE - gm*RCONST(-0.04);
202 IJth(J,1,2) = -gm*RCONST(1.0e4)*y3;
203 IJth(J,1,3) = -gm*RCONST(1.0e4)*y2;
204
205 IJth(J,2,1) = -gm*RCONST(0.04);
206 IJth(J,2,2) = ONE - gm * (RCONST(-1.0e4)*y3 - RCONST(6.0e7)*y2);
207 IJth(J,2,3) = -gm*RCONST(-1.0e4)*y2;
208
209 IJth(J,3,1) = ZERO;
210 IJth(J,3,2) = -gm*RCONST(6.0e7)*y2;
211 IJth(J,3,3) = ONE;
212
213 return(0);
214 }
215
216 /*
217 * Jacobian routine. Compute J(t,y) = df/dy. *
218 */
219
jacE(int N,realtype t,N_Vector y,N_Vector fy,DlsMat J,void * jac_data,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)220 static int jacE(int N, realtype t,
221 N_Vector y, N_Vector fy,
222 DlsMat J, void *jac_data,
223 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
224 {
225 realtype y1, y2, y3;
226
227 y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
228
229 IJth(J,1,1) = RCONST(-0.04);
230 IJth(J,1,2) = RCONST(1.0e4)*y3;
231 IJth(J,1,3) = RCONST(1.0e4)*y2;
232
233 IJth(J,2,1) = RCONST(0.04);
234 IJth(J,2,2) = RCONST(-1.0e4)*y3-RCONST(6.0e7)*y2;
235 IJth(J,2,3) = RCONST(-1.0e4)*y2;
236
237 IJth(J,3,1) = ZERO;
238 IJth(J,3,2) = RCONST(6.0e7)*y2;
239 IJth(J,3,3) = ZERO;
240
241 return(0);
242 }
243
244
245 /*
246 * g routine. Compute functions g_i(t,y) for i = 0,1.
247 */
248
g(realtype t,N_Vector y,N_Vector yp,realtype * gout,void * g_data)249 static int g(realtype t, N_Vector y, N_Vector yp, realtype *gout, void *g_data)
250 {
251 realtype y1, y3;
252
253 y1 = Ith(y,1); y3 = Ith(y,3);
254 gout[0] = y1 - RCONST(0.0001);
255 gout[1] = y3 - RCONST(0.01);
256
257 return(0);
258 }
259
260 /*
261 * Get and print some final statistics
262 */
263
PrintFinalStats(void * cpode_mem)264 static void PrintFinalStats(void *cpode_mem)
265 {
266 realtype h0u;
267 long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
268 int flag;
269
270 flag = CPodeGetActualInitStep(cpode_mem, &h0u);
271 flag = CPodeGetNumSteps(cpode_mem, &nst);
272 flag = CPodeGetNumFctEvals(cpode_mem, &nfe);
273 flag = CPodeGetNumLinSolvSetups(cpode_mem, &nsetups);
274 flag = CPodeGetNumErrTestFails(cpode_mem, &netf);
275 flag = CPodeGetNumNonlinSolvIters(cpode_mem, &nni);
276 flag = CPodeGetNumNonlinSolvConvFails(cpode_mem, &ncfn);
277
278 flag = CPDlsGetNumJacEvals(cpode_mem, &nje);
279 flag = CPDlsGetNumFctEvals(cpode_mem, &nfeLS);
280
281 flag = CPodeGetNumGEvals(cpode_mem, &nge);
282
283 printf("\nFinal Statistics:\n");
284 printf("h0u = %g\n",h0u);
285 printf("nst = %-6ld nfe = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld\n",
286 nst, nfe, nsetups, nfeLS, nje);
287 printf("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n \n",
288 nni, ncfn, netf, nge);
289 }
290
291