PartMC  2.4.0
condense_solver.c
Go to the documentation of this file.
1 /* Copyright (C) 2009-2012 Matthew West
2  * Licensed under the GNU General Public License version 2 or (at your
3  * option) any later version. See the file COPYING for details.
4  */
5 
6 /** \file
7  * \brief Interface to SUNDIALS ODE solver library for condensation.
8  */
9 
10 #include <stdlib.h>
11 #include <stdio.h>
12 #include <cvode/cvode.h>
13 #include <nvector/nvector_serial.h>
14 #include <sundials/sundials_types.h>
15 #include <cvode/cvode_impl.h>
16 
17 /** \brief Result code indicating successful completion.
18  */
19 #define PMC_CONDENSE_SOLVER_SUCCESS 0
20 /** \brief Result code indicating failure to allocate \c y vector.
21  */
22 #define PMC_CONDENSE_SOLVER_INIT_Y 1
23 /** \brief Result code indicating failure to allocate \c abstol vector.
24  */
25 #define PMC_CONDENSE_SOLVER_INIT_ABSTOL 2
26 /** \brief Result code indicating failure to create the solver.
27  */
28 #define PMC_CONDENSE_SOLVER_INIT_CVODE_MEM 3
29 /** \brief Result code indicating failure to initialize the solver.
30  */
31 #define PMC_CONDENSE_SOLVER_INIT_CVODE 4
32 /** \brief Result code indicating failure to set tolerances.
33  */
34 #define PMC_CONDENSE_SOLVER_SVTOL 5
35 /** \brief Result code indicating failure to set maximum steps.
36  */
37 #define PMC_CONDENSE_SOLVER_SET_MAX_STEPS 6
38 /** \brief Result code indicating failure of the solver.
39  */
40 #define PMC_CONDENSE_SOLVER_FAIL 7
41 
42 static int condense_vf(realtype t, N_Vector y, N_Vector ydot, void *user_data);
43 
44 static int condense_check_flag(void *flagvalue, char *funcname, int opt);
45 
46 /*******************************************************/
47 // solver block
48 static int condense_solver_Init(CVodeMem cv_mem);
49 
50 static int condense_solver_Setup(CVodeMem cv_mem, int convfail, N_Vector ypred,
51  N_Vector fpred, booleantype *jcurPtr,
52  N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3);
53 
54 static int condense_solver_Solve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
55  N_Vector ycur, N_Vector fcur);
56 
57 static void condense_solver_Free(CVodeMem cv_mem);
58 /*******************************************************/
59 
60 /** \brief Call the ODE solver.
61  *
62  * \param neq The number of equations.
63  * \param x_f A pointer to a vector of \c neq variables, giving the
64  * initial state vector on entry and the final state vector on exit.
65  * \param abstol_f A pointer to a vector of \c neq variables, giving
66  * the absolute error tolerance for the corresponding state vector
67  * component.
68  * \param reltol_f The scalar relative tolerance.
69  * \param t_initial_f The initial time (s).
70  * \param t_final_f The final time (s).
71  * \return A result code (0 is success).
72  */
73 int condense_solver(int neq, double *x_f, double *abstol_f, double reltol_f,
74  double t_initial_f, double t_final_f)
75 {
76  realtype reltol, t_initial, t_final, t, tout;
77  N_Vector y, abstol;
78  void *cvode_mem;
79  CVodeMem cv_mem;
80  int flag, i, pretype, maxl;
81  realtype *y_data, *abstol_data;
82 
83  y = abstol = NULL;
84  cvode_mem = NULL;
85 
86  y = N_VNew_Serial(neq);
87  if (condense_check_flag((void *)y, "N_VNew_Serial", 0))
89 
90  abstol = N_VNew_Serial(neq);
91  if (condense_check_flag((void *)abstol, "N_VNew_Serial", 0))
93 
94  y_data = NV_DATA_S(y);
95  abstol_data = NV_DATA_S(abstol);
96  for (i = 0; i < neq; i++) {
97  y_data[i] = x_f[i];
98  abstol_data[i] = abstol_f[i];
99  }
100 
101  reltol = reltol_f;
102  t_initial = t_initial_f;
103  t_final = t_final_f;
104 
105  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
106  if (condense_check_flag((void *)cvode_mem, "CVodeCreate", 0))
108 
109  flag = CVodeInit(cvode_mem, condense_vf, t_initial, y);
110  if (condense_check_flag(&flag, "CVodeInit", 1))
112 
113  flag = CVodeSVtolerances(cvode_mem, reltol, abstol);
114  if (condense_check_flag(&flag, "CVodeSVtolerances", 1))
116 
117  flag = CVodeSetMaxNumSteps(cvode_mem, 100000);
118  if (condense_check_flag(&flag, "CVodeSetMaxNumSteps", 1))
120 
121  /*******************************************************/
122  // dense solver
123  //flag = CVDense(cvode_mem, neq);
124  //if (condense_check_flag(&flag, "CVDense", 1)) return(1);
125  /*******************************************************/
126 
127  /*******************************************************/
128  // iterative solver
129  //pretype = PREC_LEFT;
130  //maxl = 0;
131  //flag = CVSptfqmr(cvode_mem, pretype, maxl);
132  //if (condense_check_flag(&flag, "CVSptfqmr", 1)) return(1);
133 
134  //flag = CVSpilsSetJacTimesVecFn(cvode_mem, condense_jtimes);
135  //if (condense_check_flag(&flag, "CVSpilsSetJacTimesVecFn", 1)) return(1);
136 
137  //flag = CVSpilsSetPreconditioner(cvode_mem, NULL, condense_prec);
138  //if (condense_check_flag(&flag, "CVSpilsSetPreconditioner", 1)) return(1);
139  /*******************************************************/
140 
141  /*******************************************************/
142  // explicit solver
143  cv_mem = (CVodeMem)cvode_mem;
144  cv_mem->cv_linit = condense_solver_Init;
145  cv_mem->cv_lsetup = condense_solver_Setup;
146  cv_mem->cv_lsolve = condense_solver_Solve;
147  cv_mem->cv_lfree = condense_solver_Free;
148  /*******************************************************/
149 
150  t = t_initial;
151  flag = CVode(cvode_mem, t_final, y, &t, CV_NORMAL);
152  if (condense_check_flag(&flag, "CVode", 1))
154 
155  for (i = 0; i < neq; i++) {
156  x_f[i] = y_data[i];
157  }
158 
159  N_VDestroy_Serial(y);
160  N_VDestroy_Serial(abstol);
161  CVodeFree(&cvode_mem);
163 }
164 
165 /** \brief The ODE vector field to integrate.
166  *
167  * \param t The current time (s).
168  * \param y The state vector.
169  * \param ydot The rate of change of the state vector.
170  * \param user_data A pointer to user-provided data.
171  * \return A result code (0 is success).
172  */
173 static int condense_vf(realtype t, N_Vector y, N_Vector ydot, void *user_data)
174 {
175  realtype *y_data, *ydot_data;
176  int i, neq;
177  double *y_f, *ydot_f;
178 
179  neq = NV_LENGTH_S(y);
180  y_data = NV_DATA_S(y);
181  ydot_data = NV_DATA_S(ydot);
182 
183  y_f = malloc(neq * sizeof(double));
184  ydot_f = malloc(neq * sizeof(double));
185 
186  for (i = 0; i < neq; i++) {
187  y_f[i] = y_data[i];
188  }
189  condense_vf_f(neq, t, y_f, ydot_f);
190  for (i = 0; i < neq; i++) {
191  ydot_data[i] = ydot_f[i];
192  }
193 
194  free(y_f);
195  free(ydot_f);
196  return(0);
197 }
198 
199 /** \brief Check the return value from a SUNDIALS call.
200  *
201  * - <code>opt == 0</code> means SUNDIALS function allocates memory
202  * so check if \c flagvalue is not a \c NULL pointer.
203 
204  * - <code>opt == 1</code> means SUNDIALS function returns a flag so
205  * check if the \c int pointed to by \c flagvalue has
206  * <code>flag >= 0</code>.
207 
208  * - <code>opt == 2</code> means function allocates memory so check
209  * if \c flagvalue is not a \c NULL pointer.
210  *
211  * \param flagvalue A pointer to check (either for \c NULL, or as an
212  * \c int pointer giving the flag value).
213  * \param funcname A string giving the function name returning this
214  * result code.
215  * \param opt A flag indicating the type of check to perform.
216  * \return A result code (0 is success).
217  */
218 static int condense_check_flag(void *flagvalue, char *funcname, int opt)
219 {
220  int *errflag;
221 
222  /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
223  if (opt == 0 && flagvalue == NULL) {
224  fprintf(stderr,
225  "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
226  funcname);
227  return(1); }
228 
229  /* Check if flag < 0 */
230  else if (opt == 1) {
231  errflag = (int *) flagvalue;
232  if (*errflag < 0) {
233  fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
234  funcname, *errflag);
235  return(1); }}
236 
237  /* Check if function returned NULL pointer - no memory allocated */
238  else if (opt == 2 && flagvalue == NULL) {
239  fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
240  funcname);
241  return(1); }
242 
243  return(0);
244 }
245 
246 /** \brief Initialization routine for the ODE linear equation solver.
247  *
248  * \param cv_mem The \c CVODE solver parameter structure.
249  * \return A result code (0 is success).
250  */
251 static int condense_solver_Init(CVodeMem cv_mem)
252 {
253  return(0);
254 }
255 
256 /** \brief Setup routine for the ODE linear equation solver.
257  *
258  * \param cv_mem The \c CVODE solver parameter structure.
259  * \param convfail
260  * \param ypred
261  * \param fpred
262  * \param jcurPtr
263  * \param vtemp1
264  * \param vtemp2
265  * \param vtemp3
266  * \return A result code (0 is success).
267  */
268 static int condense_solver_Setup(CVodeMem cv_mem, int convfail, N_Vector ypred,
269  N_Vector fpred, booleantype *jcurPtr,
270  N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3)
271 {
272  return(0);
273 }
274 
275 /** \brief Linear solver routine for use by the ODE solver.
276  *
277  * Should solve the system \f$(I - \gamma J) x = b\f$, where \f$J\f$
278  * is the current vector field Jacobian, \f$\gamma\f$ is a given
279  * scalar, and \f$b\f$ is a given vector.
280  *
281  * \param cv_mem The \c CVODE solver parameter structure.
282  * \param b The right-hand-side of the linear system.
283  * \param weight
284  * \param ycur The current state vector.
285  * \param fcur The current vector field vector.
286  * \return A result code (0 is success).
287  */
288 static int condense_solver_Solve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
289  N_Vector ycur, N_Vector fcur)
290 {
291  realtype *b_data, *ycur_data, *fcur_data;
292  int i, neq;
293  double *b_f, *ycur_f, *fcur_f;
294  double t, gamma;
295 
296  neq = NV_LENGTH_S(b);
297  b_data = NV_DATA_S(b);
298  ycur_data = NV_DATA_S(ycur);
299  fcur_data = NV_DATA_S(fcur);
300 
301  b_f = malloc(neq * sizeof(double));
302  ycur_f = malloc(neq * sizeof(double));
303  fcur_f = malloc(neq * sizeof(double));
304 
305  t = cv_mem->cv_tn;
306  gamma = cv_mem->cv_gamma;
307 
308  for (i = 0; i < neq; i++) {
309  b_f[i] = b_data[i];
310  ycur_f[i] = ycur_data[i];
311  fcur_f[i] = fcur_data[i];
312  }
313  condense_jac_solve_f(neq, t, ycur_f, fcur_f, b_f, gamma);
314  for (i = 0; i < neq; i++) {
315  b_data[i] = b_f[i];
316  }
317 
318  free(b_f);
319  free(ycur_f);
320  free(fcur_f);
321  return(0);
322 }
323 
324 /** \brief Finalization routine for the ODE linear equation solver.
325  *
326  * \param cv_mem The \c CVODE solver parameter structure.
327  */
328 static void condense_solver_Free(CVodeMem cv_mem)
329 {
330 }
#define PMC_CONDENSE_SOLVER_INIT_Y
Result code indicating failure to allocate y vector.
#define PMC_CONDENSE_SOLVER_SET_MAX_STEPS
Result code indicating failure to set maximum steps.
static int condense_vf(realtype t, N_Vector y, N_Vector ydot, void *user_data)
The ODE vector field to integrate.
#define PMC_CONDENSE_SOLVER_SUCCESS
Result code indicating successful completion.
#define PMC_CONDENSE_SOLVER_SVTOL
Result code indicating failure to set tolerances.
int condense_solver(int neq, double *x_f, double *abstol_f, double reltol_f, double t_initial_f, double t_final_f)
Call the ODE solver.
static int condense_solver_Setup(CVodeMem cv_mem, int convfail, N_Vector ypred, N_Vector fpred, booleantype *jcurPtr, N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3)
Setup routine for the ODE linear equation solver.
static int condense_check_flag(void *flagvalue, char *funcname, int opt)
Check the return value from a SUNDIALS call.
#define PMC_CONDENSE_SOLVER_INIT_CVODE_MEM
Result code indicating failure to create the solver.
#define PMC_CONDENSE_SOLVER_INIT_CVODE
Result code indicating failure to initialize the solver.
static int condense_solver_Solve(CVodeMem cv_mem, N_Vector b, N_Vector weight, N_Vector ycur, N_Vector fcur)
Linear solver routine for use by the ODE solver.
#define PMC_CONDENSE_SOLVER_INIT_ABSTOL
Result code indicating failure to allocate abstol vector.
static void condense_solver_Free(CVodeMem cv_mem)
Finalization routine for the ODE linear equation solver.
static int condense_solver_Init(CVodeMem cv_mem)
Initialization routine for the ODE linear equation solver.
#define PMC_CONDENSE_SOLVER_FAIL
Result code indicating failure of the solver.