1 #include <stdio.h>
2 #include <math.h>
3 #include <cpodes/cpodes.h>
4 #include <cpodes/cpodes_dense.h>
5 #include <nvector/nvector_serial.h>
6
7 #define Ith(v,i) NV_Ith_S(v,i-1)
8
9 #define RTOL RCONST(1.0e-8)
10 #define ATOL RCONST(1.0e-8)
11
12 static int f(realtype t, N_Vector y, N_Vector ydot, void *f_data);
13
14 static int g(realtype t, N_Vector y, N_Vector yp, realtype *gout, void *g_data);
15
16 static void PrintFinalStats(void *cpode_mem);
17
18
19 /*
20 *-------------------------------
21 * Main Program
22 *-------------------------------
23 */
24
main()25 int main()
26 {
27 realtype reltol, abstol, t, t0, tout, delt;
28 N_Vector yy, yp;
29 void *cpode_mem;
30 int flag;
31 int rdir[3], iroots[3];
32 realtype gout[3];
33
34 t0 = 0.0;
35 tout = 10.0;
36 delt = 0.01;
37
38 yy = N_VNew_Serial(1);
39 yp = N_VNew_Serial(1);
40 Ith(yy,1) = 0.0;
41 Ith(yp,1) = 0.0;
42
43 reltol = RTOL;
44 abstol = ATOL;
45
46 cpode_mem = CPodeCreate(CP_EXPL, CP_BDF, CP_NEWTON);
47 flag = CPodeInit(cpode_mem, (void *)f, NULL, t0, yy, NULL, CP_SS, reltol, &abstol);
48 flag = CPDense(cpode_mem, 1);
49 flag = CPodeSetStopTime(cpode_mem, tout);
50 flag = CPodeRootInit(cpode_mem, 3, g, NULL);
51
52 rdir[0] = -1;
53 rdir[1] = 0;
54 rdir[2] = 0;
55 flag = CPodeSetRootDirection(cpode_mem, rdir);
56
57 t = t0;
58 g(t, yy, yp, gout, NULL);
59 printf("%le %le %le %le %le 0 0 0\n", t, Ith(yy,1), gout[0], gout[1], gout[2]);
60
61 while(t<tout) {
62
63 flag = CPode(cpode_mem, tout, &t, yy, yp, CP_ONE_STEP_TSTOP);
64 if (flag < 0) break;
65
66 g(t, yy, yp, gout, NULL);
67 printf("%le %le %le %le %le ", t, Ith(yy,1), gout[0], gout[1], gout[2]);
68 if (flag == CP_ROOT_RETURN) {
69 CPodeGetRootInfo(cpode_mem, iroots);
70 printf(" %d %d %d\n", iroots[0], iroots[1], iroots[2]);
71 } else {
72 printf("0 0 0\n");
73 }
74
75 }
76
77 // PrintFinalStats(cpode_mem);
78
79 /* Clean-up */
80
81 N_VDestroy_Serial(yy);
82 N_VDestroy_Serial(yp);
83 CPodeFree(&cpode_mem);
84
85 return(0);
86 }
87
88
89 /*
90 *-------------------------------
91 * Functions called by the solver
92 *-------------------------------
93 */
94
f(realtype t,N_Vector y,N_Vector ydot,void * f_data)95 static int f(realtype t, N_Vector y, N_Vector ydot, void *f_data)
96 {
97 Ith(ydot,1) = t*cos(t) + sin(t);
98 Ith(ydot,1) = cos(t);
99 return(0);
100 }
101
g(realtype t,N_Vector y,N_Vector yp,realtype * gout,void * g_data)102 static int g(realtype t, N_Vector y, N_Vector yp, realtype *gout, void *g_data)
103 {
104 gout[0] = Ith(y,1);
105
106 if (Ith(y,1)+0.41 >= 0) gout[1] = +0.6;
107 else gout[1] = -0.2;
108
109 if (Ith(y,1)-0.4 >= 0) gout[2] = -0.4;
110 else if (Ith(y,1)+0.4 <= 0) gout[2] = +0.4;
111 else gout[2] = 0.0;
112
113 return(0);
114 }
115
116
PrintFinalStats(void * cpode_mem)117 static void PrintFinalStats(void *cpode_mem)
118 {
119 realtype h0u;
120 long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
121 long int nproj, nce, nsetupsP, nprf;
122 int flag;
123
124 flag = CPodeGetActualInitStep(cpode_mem, &h0u);
125 flag = CPodeGetNumSteps(cpode_mem, &nst);
126 flag = CPodeGetNumFctEvals(cpode_mem, &nfe);
127 flag = CPodeGetNumLinSolvSetups(cpode_mem, &nsetups);
128 flag = CPodeGetNumErrTestFails(cpode_mem, &netf);
129 flag = CPodeGetNumNonlinSolvIters(cpode_mem, &nni);
130 flag = CPodeGetNumNonlinSolvConvFails(cpode_mem, &ncfn);
131
132 flag = CPDlsGetNumJacEvals(cpode_mem, &nje);
133 flag = CPDlsGetNumFctEvals(cpode_mem, &nfeLS);
134
135 flag = CPodeGetProjStats(cpode_mem, &nproj, &nce, &nsetupsP, &nprf);
136
137 flag = CPodeGetNumGEvals(cpode_mem, &nge);
138
139 printf("\nFinal Statistics:\n");
140 printf("h0u = %g\n",h0u);
141 printf("nst = %-6ld nfe = %-6ld nsetups = %-6ld\n",
142 nst, nfe, nsetups);
143 printf("nfeLS = %-6ld nje = %ld\n",
144 nfeLS, nje);
145 printf("nni = %-6ld ncfn = %-6ld netf = %-6ld \n",
146 nni, ncfn, netf);
147 printf("nproj = %-6ld nce = %-6ld nsetupsP = %-6ld nprf = %-6ld\n",
148 nproj, nce, nsetupsP, nprf);
149 printf("nge = %ld\n", nge);
150 }
151
152