1212#include <cvode/cvode.h>
1313#include <nvector/nvector_serial.h>
1414#include <sundials/sundials_types.h>
15- #include <cvode/cvode_impl .h>
15+ #include <sunlinsol/sunlinsol_spgmr .h>
1616
1717/** \brief Result code indicating successful completion.
1818 */
3838/** \brief Result code indicating failure of the solver.
3939 */
4040#define PMC_CONDENSE_SOLVER_FAIL 7
41+ /** \brief Result code indicating failure in creation of the linear solver.
42+ */
43+ #define PMC_CONDENSE_SOLVER_LINSOL_CTOR 8
44+ /** \brief Result code indicating failure in setting the linear solver.
45+ */
46+ #define PMC_CONDENSE_SOLVER_LINSOL_SET 9
47+ /** \brief Result code indicating failure in setting the preconditioner.
48+ */
49+ #define PMC_CONDENSE_SOLVER_LINSOL_PREC 10
4150
4251static int condense_vf (realtype t , N_Vector y , N_Vector ydot , void * user_data );
4352
4453static int condense_check_flag (void * flagvalue , char * funcname , int opt );
4554
4655/*******************************************************/
4756// 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 );
5357
54- static int condense_solver_Solve (CVodeMem cv_mem , N_Vector b , N_Vector weight ,
55- N_Vector ycur , N_Vector fcur );
58+ static int condense_solver_Solve (double t , N_Vector ycur , N_Vector fcur ,
59+ N_Vector r , N_Vector b , double gamma , double delta , int lr , void * user_data );
5660
57- static void condense_solver_Free (CVodeMem cv_mem );
5861/*******************************************************/
5962
63+ void condense_vf_f (int neq , realtype t , double * y_f , double * ydot_f );
64+ void condense_jac_solve_f (int neq , double t , double * ycur_f , double * fcur_f , double * b_f , double gamma );
65+
6066/** \brief Call the ODE solver.
6167 *
6268 * \param neq The number of equations.
@@ -73,11 +79,10 @@ static void condense_solver_Free(CVodeMem cv_mem);
7379int condense_solver (int neq , double * x_f , double * abstol_f , double reltol_f ,
7480 double t_initial_f , double t_final_f )
7581{
76- realtype reltol , t_initial , t_final , t , tout ;
82+ realtype reltol , t_initial , t_final , t ;
7783 N_Vector y , abstol ;
7884 void * cvode_mem ;
79- CVodeMem cv_mem ;
80- int flag , i , pretype , maxl ;
85+ int flag , i ;
8186 realtype * y_data , * abstol_data ;
8287
8388 y = abstol = NULL ;
@@ -102,7 +107,7 @@ int condense_solver(int neq, double *x_f, double *abstol_f, double reltol_f,
102107 t_initial = t_initial_f ;
103108 t_final = t_final_f ;
104109
105- cvode_mem = CVodeCreate (CV_BDF , CV_NEWTON );
110+ cvode_mem = CVodeCreate (CV_BDF );
106111 if (condense_check_flag ((void * )cvode_mem , "CVodeCreate" , 0 ))
107112 return PMC_CONDENSE_SOLVER_INIT_CVODE_MEM ;
108113
@@ -118,34 +123,16 @@ int condense_solver(int neq, double *x_f, double *abstol_f, double reltol_f,
118123 if (condense_check_flag (& flag , "CVodeSetMaxNumSteps" , 1 ))
119124 return PMC_CONDENSE_SOLVER_SET_MAX_STEPS ;
120125
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- /*******************************************************/
126+
127+ SUNLinearSolver LS = SUNLinSol_SPGMR (y , PREC_LEFT , 0 );
128+ if (condense_check_flag ((void * )LS , "SUNLinSol_SPGMR" , 0 ))
129+ return PMC_CONDENSE_SOLVER_LINSOL_CTOR ;
130+ flag = CVodeSetLinearSolver (cvode_mem , LS , NULL );
131+ if (condense_check_flag (& flag , "CVodeSetLinearSolver" , 1 ))
132+ return PMC_CONDENSE_SOLVER_LINSOL_SET ;
133+ flag = CVodeSetPreconditioner (cvode_mem , NULL , condense_solver_Solve );
134+ if (condense_check_flag (& flag , "CVodeSetPreconditioner" , 1 ))
135+ return PMC_CONDENSE_SOLVER_LINSOL_PREC ;
149136
150137 t = t_initial ;
151138 flag = CVode (cvode_mem , t_final , y , & t , CV_NORMAL );
@@ -170,6 +157,8 @@ int condense_solver(int neq, double *x_f, double *abstol_f, double reltol_f,
170157 * \param user_data A pointer to user-provided data.
171158 * \return A result code (0 is success).
172159 */
160+ #pragma GCC diagnostic push
161+ #pragma GCC diagnostic ignored "-Wunused-parameter"
173162static int condense_vf (realtype t , N_Vector y , N_Vector ydot , void * user_data )
174163{
175164 realtype * y_data , * ydot_data ;
@@ -195,6 +184,7 @@ static int condense_vf(realtype t, N_Vector y, N_Vector ydot, void *user_data)
195184 free (ydot_f );
196185 return (0 );
197186}
187+ #pragma GCC diagnostic pop
198188
199189/** \brief Check the return value from a SUNDIALS call.
200190 *
@@ -243,88 +233,45 @@ static int condense_check_flag(void *flagvalue, char *funcname, int opt)
243233 return (0 );
244234}
245235
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-
275236/** \brief Linear solver routine for use by the ODE solver.
276237 *
277238 * Should solve the system \f$(I - \gamma J) x = b\f$, where \f$J\f$
278239 * is the current vector field Jacobian, \f$\gamma\f$ is a given
279240 * scalar, and \f$b\f$ is a given vector.
280241 *
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).
287242 */
288- static int condense_solver_Solve (CVodeMem cv_mem , N_Vector b , N_Vector weight ,
289- N_Vector ycur , N_Vector fcur )
243+ #pragma GCC diagnostic push
244+ #pragma GCC diagnostic ignored "-Wunused-parameter"
245+ static int condense_solver_Solve (double t , N_Vector ycur , N_Vector fcur ,
246+ N_Vector b , N_Vector z , double gamma , double delta , int lr , void * user_data )
290247{
291- realtype * b_data , * ycur_data , * fcur_data ;
248+ realtype * b_data , * ycur_data , * fcur_data , * z_data ;
292249 int i , neq ;
293250 double * b_f , * ycur_f , * fcur_f ;
294- double t , gamma ;
295251
296252 neq = NV_LENGTH_S (b );
297253 b_data = NV_DATA_S (b );
254+ z_data = NV_DATA_S (z );
298255 ycur_data = NV_DATA_S (ycur );
299256 fcur_data = NV_DATA_S (fcur );
300257
301258 b_f = malloc (neq * sizeof (double ));
302259 ycur_f = malloc (neq * sizeof (double ));
303260 fcur_f = malloc (neq * sizeof (double ));
304261
305- t = cv_mem -> cv_tn ;
306- gamma = cv_mem -> cv_gamma ;
307-
308262 for (i = 0 ; i < neq ; i ++ ) {
309263 b_f [i ] = b_data [i ];
310264 ycur_f [i ] = ycur_data [i ];
311265 fcur_f [i ] = fcur_data [i ];
312266 }
313267 condense_jac_solve_f (neq , t , ycur_f , fcur_f , b_f , gamma );
314268 for (i = 0 ; i < neq ; i ++ ) {
315- b_data [i ] = b_f [i ];
269+ z_data [i ] = b_f [i ];
316270 }
317271
318272 free (b_f );
319273 free (ycur_f );
320274 free (fcur_f );
321275 return (0 );
322276}
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- }
277+ #pragma GCC diagnostic pop
0 commit comments