12 #include <cvode/cvode.h>
13 #include <nvector/nvector_serial.h>
14 #include <sundials/sundials_types.h>
15 #include <sunlinsol/sunlinsol_spgmr.h>
19 #define PMC_CONDENSE_SOLVER_SUCCESS 0
22 #define PMC_CONDENSE_SOLVER_INIT_Y 1
25 #define PMC_CONDENSE_SOLVER_INIT_ABSTOL 2
28 #define PMC_CONDENSE_SOLVER_INIT_CVODE_MEM 3
31 #define PMC_CONDENSE_SOLVER_INIT_CVODE 4
34 #define PMC_CONDENSE_SOLVER_SVTOL 5
37 #define PMC_CONDENSE_SOLVER_SET_MAX_STEPS 6
40 #define PMC_CONDENSE_SOLVER_FAIL 7
43 #define PMC_CONDENSE_SOLVER_LINSOL_CTOR 8
46 #define PMC_CONDENSE_SOLVER_LINSOL_SET 9
49 #define PMC_CONDENSE_SOLVER_LINSOL_PREC 10
51 static int condense_vf(realtype t, N_Vector y, N_Vector ydot,
void *user_data);
59 N_Vector r, N_Vector b,
double gamma,
double delta,
int lr,
void *user_data);
80 double t_initial_f,
double t_final_f)
82 realtype reltol, t_initial, t_final, t;
86 realtype *y_data, *abstol_data;
91 y = N_VNew_Serial(neq);
95 abstol = N_VNew_Serial(neq);
99 y_data = NV_DATA_S(y);
100 abstol_data = NV_DATA_S(abstol);
101 for (i = 0; i < neq; i++) {
103 abstol_data[i] = abstol_f[i];
107 t_initial = t_initial_f;
110 cvode_mem = CVodeCreate(CV_BDF);
114 flag = CVodeInit(cvode_mem,
condense_vf, t_initial, y);
118 flag = CVodeSVtolerances(cvode_mem, reltol, abstol);
122 flag = CVodeSetMaxNumSteps(cvode_mem, 100000);
127 SUNLinearSolver LS = SUNLinSol_SPGMR(y, PREC_LEFT, 0);
130 flag = CVodeSetLinearSolver(cvode_mem, LS, NULL);
138 flag = CVode(cvode_mem, t_final, y, &t, CV_NORMAL);
142 for (i = 0; i < neq; i++) {
146 N_VDestroy_Serial(y);
147 N_VDestroy_Serial(abstol);
148 CVodeFree(&cvode_mem);
160 #pragma GCC diagnostic push
161 #pragma GCC diagnostic ignored "-Wunused-parameter"
162 static int condense_vf(realtype t, N_Vector y, N_Vector ydot,
void *user_data)
164 realtype *y_data, *ydot_data;
166 double *y_f, *ydot_f;
168 neq = NV_LENGTH_S(y);
169 y_data = NV_DATA_S(y);
170 ydot_data = NV_DATA_S(ydot);
172 y_f = malloc(neq *
sizeof(
double));
173 ydot_f = malloc(neq *
sizeof(
double));
175 for (i = 0; i < neq; i++) {
179 for (i = 0; i < neq; i++) {
180 ydot_data[i] = ydot_f[i];
187 #pragma GCC diagnostic pop
213 if (opt == 0 && flagvalue == NULL) {
215 "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
221 errflag = (
int *) flagvalue;
223 fprintf(stderr,
"\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
228 else if (opt == 2 && flagvalue == NULL) {
229 fprintf(stderr,
"\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
243 #pragma GCC diagnostic push
244 #pragma GCC diagnostic ignored "-Wunused-parameter"
246 N_Vector b, N_Vector z,
double gamma,
double delta,
int lr,
void *user_data)
248 realtype *b_data, *ycur_data, *fcur_data, *z_data;
250 double *b_f, *ycur_f, *fcur_f;
252 neq = NV_LENGTH_S(b);
253 b_data = NV_DATA_S(b);
254 z_data = NV_DATA_S(z);
255 ycur_data = NV_DATA_S(ycur);
256 fcur_data = NV_DATA_S(fcur);
258 b_f = malloc(neq *
sizeof(
double));
259 ycur_f = malloc(neq *
sizeof(
double));
260 fcur_f = malloc(neq *
sizeof(
double));
262 for (i = 0; i < neq; i++) {
264 ycur_f[i] = ycur_data[i];
265 fcur_f[i] = fcur_data[i];
268 for (i = 0; i < neq; i++) {
277 #pragma GCC diagnostic pop