1 /*-----------------------------------------------------------------
2 * Programmer(s): Daniel R. Reynolds @ SMU
3 *---------------------------------------------------------------
4 * SUNDIALS Copyright Start
5 * Copyright (c) 2002-2020, Lawrence Livermore National Security
6 * and Southern Methodist University.
7 * All rights reserved.
8 *
9 * See the top-level LICENSE and NOTICE files for details.
10 *
11 * SPDX-License-Identifier: BSD-3-Clause
12 * SUNDIALS Copyright End
13 *---------------------------------------------------------------
14 * Example problem:
15 *
16 * The following is a simple example problem with analytical
17 * solution,
18 * dy/dt = (t+1)*exp(-y)
19 * for t in the interval [0.0, 10.0], with initial condition: y=0.
20 * This has analytical solution
21 * y(t) = log(0.5*t^2 + t + 1)
22 *
23 * This program solves the problem with the ERK method.
24 * Output is printed every 1.0 units of time (10 total).
25 * Run statistics (optional outputs) are printed at the end.
26 *-----------------------------------------------------------------*/
27
28 /* Header files */
29 #include <stdio.h>
30 #include <math.h>
31 #include <arkode/arkode_erkstep.h> /* prototypes for ERKStep fcts., consts */
32 #include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
33 #include <sundials/sundials_types.h> /* def. of type 'realtype' */
34 #include <sundials/sundials_math.h> /* def. of SUNRsqrt, etc. */
35
36 #if defined(SUNDIALS_EXTENDED_PRECISION)
37 #define GSYM "Lg"
38 #define ESYM "Le"
39 #define FSYM "Lf"
40 #else
41 #define GSYM "g"
42 #define ESYM "e"
43 #define FSYM "f"
44 #endif
45
46 /* User-supplied Functions Called by the Solver */
47 static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data);
48
49 /* Private function to check function return values */
50 static int check_flag(void *flagvalue, const char *funcname, int opt);
51
52 /* Main Program */
main()53 int main()
54 {
55 /* general problem parameters */
56 realtype T0 = RCONST(0.0); /* initial time */
57 realtype Tf = RCONST(10.0); /* final time */
58 realtype dTout = RCONST(1.0); /* time between outputs */
59 sunindextype NEQ = 1; /* number of dependent vars. */
60 realtype reltol = 1.0e-6; /* tolerances */
61 realtype abstol = 1.0e-10;
62
63 /* general problem variables */
64 int flag; /* reusable error-checking flag */
65 N_Vector y = NULL; /* empty vector for storing solution */
66 void *arkode_mem = NULL; /* empty ARKode memory structure */
67 FILE *UFID;
68 realtype t, tout;
69 long int nst, nst_a, nfe, netf;
70
71 /* Initial problem output */
72 printf("\nAnalytical ODE test problem:\n");
73 printf(" reltol = %.1"ESYM"\n", reltol);
74 printf(" abstol = %.1"ESYM"\n\n",abstol);
75
76 /* Initialize data structures */
77 y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
78 if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
79 NV_Ith_S(y,0) = 0.0; /* Specify initial condition */
80
81 /* Call ERKStepCreate to initialize the ERK timestepper module and
82 specify the right-hand side function in y'=f(t,y), the inital time
83 T0, and the initial dependent variable vector y. */
84 arkode_mem = ERKStepCreate(f, T0, y);
85 if (check_flag((void *)arkode_mem, "ERKStepCreate", 0)) return 1;
86
87 /* Specify tolerances */
88 flag = ERKStepSStolerances(arkode_mem, reltol, abstol);
89 if (check_flag(&flag, "ERKStepSStolerances", 1)) return 1;
90
91 /* Open output stream for results, output comment line */
92 UFID = fopen("solution.txt","w");
93 fprintf(UFID,"# t u\n");
94
95 /* output initial condition to disk */
96 fprintf(UFID," %.16"ESYM" %.16"ESYM"\n", T0, NV_Ith_S(y,0));
97
98 /* Main time-stepping loop: calls ERKStepEvolve to perform the integration, then
99 prints results. Stops when the final time has been reached */
100 t = T0;
101 tout = T0+dTout;
102 printf(" t u\n");
103 printf(" ---------------------\n");
104 while (Tf - t > 1.0e-15) {
105
106 flag = ERKStepEvolve(arkode_mem, tout, y, &t, ARK_NORMAL); /* call integrator */
107 if (check_flag(&flag, "ERKStepEvolve", 1)) break;
108 printf(" %10.6"FSYM" %10.6"FSYM"\n", t, NV_Ith_S(y,0)); /* access/print solution */
109 fprintf(UFID," %.16"ESYM" %.16"ESYM"\n", t, NV_Ith_S(y,0));
110 if (flag >= 0) { /* successful solve: update time */
111 tout += dTout;
112 tout = (tout > Tf) ? Tf : tout;
113 } else { /* unsuccessful solve: break */
114 fprintf(stderr,"Solver failure, stopping integration\n");
115 break;
116 }
117 }
118 printf(" ---------------------\n");
119 fclose(UFID);
120
121 /* Print some final statistics */
122 flag = ERKStepGetNumSteps(arkode_mem, &nst);
123 check_flag(&flag, "ERKStepGetNumSteps", 1);
124 flag = ERKStepGetNumStepAttempts(arkode_mem, &nst_a);
125 check_flag(&flag, "ERKStepGetNumStepAttempts", 1);
126 flag = ERKStepGetNumRhsEvals(arkode_mem, &nfe);
127 check_flag(&flag, "ERKStepGetNumRhsEvals", 1);
128 flag = ERKStepGetNumErrTestFails(arkode_mem, &netf);
129 check_flag(&flag, "ERKStepGetNumErrTestFails", 1);
130
131 printf("\nFinal Solver Statistics:\n");
132 printf(" Internal solver steps = %li (attempted = %li)\n", nst, nst_a);
133 printf(" Total RHS evals = %li\n", nfe);
134 printf(" Total number of error test failures = %li\n\n", netf);
135
136 /* Clean up and return with successful completion */
137 N_VDestroy(y); /* Free y vector */
138 ERKStepFree(&arkode_mem); /* Free integrator memory */
139 return 0;
140 }
141
142 /*-------------------------------
143 * Functions called by the solver
144 *-------------------------------*/
145
146 /* f routine to compute the ODE RHS function f(t,y). */
f(realtype t,N_Vector y,N_Vector ydot,void * user_data)147 static int f(realtype t, N_Vector y, N_Vector ydot, void *user_data)
148 {
149 NV_Ith_S(ydot,0) = (t+1.0)*SUNRexp(-NV_Ith_S(y,0));
150 return 0;
151 }
152
153 /*-------------------------------
154 * Private helper functions
155 *-------------------------------*/
156
157 /* Check function return value...
158 opt == 0 means SUNDIALS function allocates memory so check if
159 returned NULL pointer
160 opt == 1 means SUNDIALS function returns a flag so check if
161 flag >= 0
162 opt == 2 means function allocates memory so check if returned
163 NULL pointer
164 */
check_flag(void * flagvalue,const char * funcname,int opt)165 static int check_flag(void *flagvalue, const char *funcname, int opt)
166 {
167 int *errflag;
168
169 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
170 if (opt == 0 && flagvalue == NULL) {
171 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
172 funcname);
173 return 1; }
174
175 /* Check if flag < 0 */
176 else if (opt == 1) {
177 errflag = (int *) flagvalue;
178 if (*errflag < 0) {
179 fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
180 funcname, *errflag);
181 return 1; }}
182
183 /* Check if function returned NULL pointer - no memory allocated */
184 else if (opt == 2 && flagvalue == NULL) {
185 fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
186 funcname);
187 return 1; }
188
189 return 0;
190 }
191
192
193 /*---- end of file ----*/
194