1 /*
2 * -----------------------------------------------------------------
3 * Programmer(s): Daniel R. Reynolds @ SMU
4 * Alan C. Hindmarsh, Radu Serban and Aaron Collier @ LLNL
5 * -----------------------------------------------------------------
6 * SUNDIALS Copyright Start
7 * Copyright (c) 2002-2020, Lawrence Livermore National Security
8 * and Southern Methodist University.
9 * All rights reserved.
10 *
11 * See the top-level LICENSE and NOTICE files for details.
12 *
13 * SPDX-License-Identifier: BSD-3-Clause
14 * SUNDIALS Copyright End
15 * -----------------------------------------------------------------
16 * This is the implementation file for the Fortran interface to
17 * the CVODE package. See fcvode.h for usage.
18 * NOTE: some routines are necessarily stored elsewhere to avoid
19 * linking problems. Therefore, see the othe C files in this folder
20 * for all the options available.
21 * -----------------------------------------------------------------
22 */
23
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <string.h>
27
28 #include "fcvode.h" /* actual function names, prototypes, global vars.*/
29 #include "cvode_impl.h" /* definition of CVodeMem type */
30 #include <sundials/sundials_matrix.h>
31 #include <cvode/cvode_ls.h>
32 #include <cvode/cvode_diag.h>
33
34
35 /***************************************************************************/
36
37 /* Definitions for global variables shared amongst various routines */
38
39 void *CV_cvodemem;
40 long int *CV_iout;
41 realtype *CV_rout;
42 int CV_nrtfn;
43 int CV_ls;
44
45 /***************************************************************************/
46
47 /* private constant(s) */
48 #define ZERO RCONST(0.0)
49
50 /***************************************************************************/
51
52 /* Prototypes of the Fortran routines */
53
54 #ifdef __cplusplus /* wrapper to enable C++ usage */
55 extern "C" {
56 #endif
57 extern void FCV_FUN(realtype*, /* T */
58 realtype*, /* Y */
59 realtype*, /* YDOT */
60 long int*, /* IPAR */
61 realtype*, /* RPAR */
62 int*); /* IER */
63 #ifdef __cplusplus
64 }
65 #endif
66
67 /**************************************************************************/
68
FCV_MALLOC(realtype * t0,realtype * y0,int * meth,int * iatol,realtype * rtol,realtype * atol,long int * iout,realtype * rout,long int * ipar,realtype * rpar,int * ier)69 void FCV_MALLOC(realtype *t0, realtype *y0,
70 int *meth, int *iatol,
71 realtype *rtol, realtype *atol,
72 long int *iout, realtype *rout,
73 long int *ipar, realtype *rpar,
74 int *ier)
75 {
76 int lmm;
77 N_Vector Vatol;
78 FCVUserData CV_userdata;
79
80 *ier = 0;
81
82 /* Check for required vector operations */
83 if(F2C_CVODE_vec->ops->nvgetarraypointer == NULL ||
84 F2C_CVODE_vec->ops->nvsetarraypointer == NULL) {
85 *ier = -1;
86 fprintf(stderr, "A required vector operation is not implemented.\n\n");
87 return;
88 }
89
90 /* Initialize all pointers to NULL */
91 CV_cvodemem = NULL;
92 Vatol = NULL;
93 FCVNullNonlinSol();
94
95 /* initialize global constants to disable each option */
96 CV_nrtfn = 0;
97 CV_ls = -1;
98
99 /* Create CVODE object */
100 lmm = (*meth == 1) ? CV_ADAMS : CV_BDF;
101 CV_cvodemem = CVodeCreate(lmm);
102 if (CV_cvodemem == NULL) {
103 *ier = -1;
104 return;
105 }
106
107 /* Set and attach user data */
108 CV_userdata = NULL;
109 CV_userdata = (FCVUserData) malloc(sizeof *CV_userdata);
110 if (CV_userdata == NULL) {
111 *ier = -1;
112 return;
113 }
114 CV_userdata->rpar = rpar;
115 CV_userdata->ipar = ipar;
116
117 *ier = CVodeSetUserData(CV_cvodemem, CV_userdata);
118 if(*ier != CV_SUCCESS) {
119 free(CV_userdata); CV_userdata = NULL;
120 *ier = -1;
121 return;
122 }
123
124 /* Set data in F2C_CVODE_vec to y0 */
125 N_VSetArrayPointer(y0, F2C_CVODE_vec);
126
127 /* Call CVodeInit */
128 *ier = CVodeInit(CV_cvodemem, FCVf, *t0, F2C_CVODE_vec);
129
130 /* Reset data pointers */
131 N_VSetArrayPointer(NULL, F2C_CVODE_vec);
132
133 /* On failure, exit */
134 if(*ier != CV_SUCCESS) {
135 free(CV_userdata); CV_userdata = NULL;
136 *ier = -1;
137 return;
138 }
139
140 /* Set tolerances */
141 switch (*iatol) {
142 case 1:
143 *ier = CVodeSStolerances(CV_cvodemem, *rtol, *atol);
144 break;
145 case 2:
146 Vatol = NULL;
147 Vatol = N_VCloneEmpty(F2C_CVODE_vec);
148 if (Vatol == NULL) {
149 free(CV_userdata); CV_userdata = NULL;
150 *ier = -1;
151 return;
152 }
153 N_VSetArrayPointer(atol, Vatol);
154 *ier = CVodeSVtolerances(CV_cvodemem, *rtol, Vatol);
155 N_VDestroy(Vatol);
156 break;
157 }
158
159 /* On failure, exit */
160 if(*ier != CV_SUCCESS) {
161 free(CV_userdata); CV_userdata = NULL;
162 *ier = -1;
163 return;
164 }
165
166 /* Grab optional output arrays and store them in global variables */
167 CV_iout = iout;
168 CV_rout = rout;
169
170 /* Store the unit roundoff in rout for user access */
171 CV_rout[5] = UNIT_ROUNDOFF;
172
173 return;
174 }
175
176 /***************************************************************************/
177
FCV_REINIT(realtype * t0,realtype * y0,int * iatol,realtype * rtol,realtype * atol,int * ier)178 void FCV_REINIT(realtype *t0, realtype *y0,
179 int *iatol, realtype *rtol, realtype *atol,
180 int *ier)
181 {
182 N_Vector Vatol;
183
184 *ier = 0;
185
186 /* Initialize all pointers to NULL */
187 Vatol = NULL;
188
189 /* Set data in F2C_CVODE_vec to y0 */
190 N_VSetArrayPointer(y0, F2C_CVODE_vec);
191
192 /* Call CVReInit */
193 *ier = CVodeReInit(CV_cvodemem, *t0, F2C_CVODE_vec);
194
195 /* Reset data pointers */
196 N_VSetArrayPointer(NULL, F2C_CVODE_vec);
197
198 /* On failure, exit */
199 if (*ier != CV_SUCCESS) {
200 *ier = -1;
201 return;
202 }
203
204 /* Set tolerances */
205 switch (*iatol) {
206 case 1:
207 *ier = CVodeSStolerances(CV_cvodemem, *rtol, *atol);
208 break;
209 case 2:
210 Vatol = NULL;
211 Vatol = N_VCloneEmpty(F2C_CVODE_vec);
212 if (Vatol == NULL) {
213 *ier = -1;
214 return;
215 }
216 N_VSetArrayPointer(atol, Vatol);
217 *ier = CVodeSVtolerances(CV_cvodemem, *rtol, Vatol);
218 N_VDestroy(Vatol);
219 break;
220 }
221
222 /* On failure, exit */
223 if (*ier != CV_SUCCESS) {
224 *ier = -1;
225 return;
226 }
227
228 return;
229 }
230
231 /***************************************************************************/
232
FCV_SETIIN(char key_name[],long int * ival,int * ier)233 void FCV_SETIIN(char key_name[], long int *ival, int *ier)
234 {
235 if (!strncmp(key_name,"MAX_ORD",7))
236 *ier = CVodeSetMaxOrd(CV_cvodemem, (int) *ival);
237 else if (!strncmp(key_name,"MAX_NSTEPS",10))
238 *ier = CVodeSetMaxNumSteps(CV_cvodemem, (long int) *ival);
239 else if (!strncmp(key_name,"MAX_ERRFAIL",11))
240 *ier = CVodeSetMaxErrTestFails(CV_cvodemem, (int) *ival);
241 else if (!strncmp(key_name,"MAX_NITERS",10))
242 *ier = CVodeSetMaxNonlinIters(CV_cvodemem, (int) *ival);
243 else if (!strncmp(key_name,"MAX_CONVFAIL",12))
244 *ier = CVodeSetMaxConvFails(CV_cvodemem, (int) *ival);
245 else if (!strncmp(key_name,"HNIL_WARNS",10))
246 *ier = CVodeSetMaxHnilWarns(CV_cvodemem, (int) *ival);
247 else if (!strncmp(key_name,"STAB_LIM",8))
248 *ier = CVodeSetStabLimDet(CV_cvodemem, (booleantype) *ival);
249 else {
250 *ier = -99;
251 fprintf(stderr, "FCVSETIIN: Unrecognized key.\n\n");
252 }
253
254 }
255
256 /***************************************************************************/
257
FCV_SETRIN(char key_name[],realtype * rval,int * ier)258 void FCV_SETRIN(char key_name[], realtype *rval, int *ier)
259 {
260 if (!strncmp(key_name,"INIT_STEP",9))
261 *ier = CVodeSetInitStep(CV_cvodemem, *rval);
262 else if (!strncmp(key_name,"MAX_STEP",8))
263 *ier = CVodeSetMaxStep(CV_cvodemem, *rval);
264 else if (!strncmp(key_name,"MIN_STEP",8))
265 *ier = CVodeSetMinStep(CV_cvodemem, *rval);
266 else if (!strncmp(key_name,"STOP_TIME",9))
267 *ier = CVodeSetStopTime(CV_cvodemem, *rval);
268 else if (!strncmp(key_name,"NLCONV_COEF",11))
269 *ier = CVodeSetNonlinConvCoef(CV_cvodemem, *rval);
270 else {
271 *ier = -99;
272 fprintf(stderr, "FCVSETRIN: Unrecognized key.\n\n");
273 }
274
275 }
276
277 /***************************************************************************/
278
FCV_SETVIN(char key_name[],realtype * vval,int * ier)279 void FCV_SETVIN(char key_name[], realtype *vval, int *ier)
280 {
281 N_Vector Vec;
282
283 *ier = 0;
284
285 if (!strncmp(key_name,"CONSTR_VEC",10)) {
286 Vec = NULL;
287 Vec = N_VCloneEmpty(F2C_CVODE_vec);
288 if (Vec == NULL) {
289 *ier = -1;
290 return;
291 }
292 N_VSetArrayPointer(vval, Vec);
293 *ier = CVodeSetConstraints(CV_cvodemem, Vec);
294 N_VDestroy(Vec);
295 }
296 else {
297 *ier = -99;
298 fprintf(stderr, "FCVSETVIN: Unrecognized key. \n\n");
299 }
300
301 }
302
303 /***************************************************************************/
304
FCV_LSINIT(int * ier)305 void FCV_LSINIT(int *ier) {
306 if ( (CV_cvodemem == NULL) || (F2C_CVODE_linsol == NULL) ) {
307 *ier = -1;
308 return;
309 }
310 *ier = CVodeSetLinearSolver(CV_cvodemem, F2C_CVODE_linsol,
311 F2C_CVODE_matrix);
312 CV_ls = CV_LS_STD;
313 return;
314 }
315
316 /***************************************************************************/
317
318 /* ---DEPRECATED--- */
FCV_DLSINIT(int * ier)319 void FCV_DLSINIT(int *ier)
320 { FCV_LSINIT(ier); }
321
322 /***************************************************************************/
323
324 /* ---DEPRECATED--- */
FCV_SPILSINIT(int * ier)325 void FCV_SPILSINIT(int *ier)
326 { FCV_LSINIT(ier); }
327
328 /*=============================================================*/
329
330 /* ---DEPRECATED--- */
FCV_SPILSSETEPSLIN(realtype * eplifac,int * ier)331 void FCV_SPILSSETEPSLIN(realtype *eplifac, int *ier)
332 { FCV_LSSETEPSLIN(eplifac, ier); }
333
FCV_LSSETEPSLIN(realtype * eplifac,int * ier)334 void FCV_LSSETEPSLIN(realtype *eplifac, int *ier)
335 { *ier = CVodeSetEpsLin(CV_cvodemem, *eplifac); }
336
337 /***************************************************************************/
338
FCV_DIAG(int * ier)339 void FCV_DIAG(int *ier)
340 {
341 if (CV_cvodemem == NULL) {
342 *ier = -1;
343 return;
344 }
345 *ier = CVDiag(CV_cvodemem);
346 CV_ls = CV_LS_DIAG;
347 return;
348 }
349
350 /***************************************************************************/
351
FCV_CVODE(realtype * tout,realtype * t,realtype * y,int * itask,int * ier)352 void FCV_CVODE(realtype *tout, realtype *t, realtype *y, int *itask, int *ier)
353 {
354 /*
355 tout is the t value where output is desired
356 F2C_CVODE_vec is the N_Vector containing the solution on return
357 t is the returned independent variable value
358 itask is the task indicator (1 = CV_NORMAL, 2 = CV_ONE_STEP,
359 3 = CV_NORMAL_TSTOP, 4 = CV_ONE_STEP_TSTOP)
360 */
361
362 int qu, qcur;
363
364 N_VSetArrayPointer(y, F2C_CVODE_vec);
365
366 *ier = CVode(CV_cvodemem, *tout, F2C_CVODE_vec, t, *itask);
367
368 N_VSetArrayPointer(NULL, F2C_CVODE_vec);
369
370 /* Load optional outputs in iout & rout */
371 CVodeGetWorkSpace(CV_cvodemem,
372 &CV_iout[0], /* LENRW */
373 &CV_iout[1]); /* LENIW */
374 CVodeGetIntegratorStats(CV_cvodemem,
375 &CV_iout[2], /* NST */
376 &CV_iout[3], /* NFE */
377 &CV_iout[7], /* NSETUPS */
378 &CV_iout[4], /* NETF */
379 &qu, /* QU */
380 &qcur, /* QCUR */
381 &CV_rout[0], /* H0U */
382 &CV_rout[1], /* HU */
383 &CV_rout[2], /* HCUR */
384 &CV_rout[3]); /* TCUR */
385 CV_iout[8] = (long int) qu;
386 CV_iout[9] = (long int) qcur;
387 CVodeGetTolScaleFactor(CV_cvodemem,
388 &CV_rout[4]); /* TOLSFAC */
389 CVodeGetNonlinSolvStats(CV_cvodemem,
390 &CV_iout[6], /* NNI */
391 &CV_iout[5]); /* NCFN */
392 CVodeGetNumStabLimOrderReds(CV_cvodemem, &CV_iout[10]); /* NOR */
393
394 /* Root finding is on */
395 if (CV_nrtfn != 0)
396 CVodeGetNumGEvals(CV_cvodemem, &CV_iout[11]); /* NGE */
397
398 switch(CV_ls) {
399 case CV_LS_STD:
400 CVodeGetLinWorkSpace(CV_cvodemem, &CV_iout[12], &CV_iout[13]); /* LENRWLS,LENIWLS */
401 CVodeGetLastLinFlag(CV_cvodemem, &CV_iout[14]); /* LSTF */
402 CVodeGetNumLinRhsEvals(CV_cvodemem, &CV_iout[15]); /* NFELS */
403 CVodeGetNumJacEvals(CV_cvodemem, &CV_iout[16]); /* NJE */
404 CVodeGetNumJTSetupEvals(CV_cvodemem, &CV_iout[17]); /* NJTS */
405 CVodeGetNumJtimesEvals(CV_cvodemem, &CV_iout[18]); /* NJTV */
406 CVodeGetNumPrecEvals(CV_cvodemem, &CV_iout[19]); /* NPE */
407 CVodeGetNumPrecSolves(CV_cvodemem, &CV_iout[20]); /* NPS */
408 CVodeGetNumLinIters(CV_cvodemem, &CV_iout[21]); /* NLI */
409 CVodeGetNumLinConvFails(CV_cvodemem, &CV_iout[22]); /* NCFL */
410 break;
411 case CV_LS_DIAG:
412 CVDiagGetWorkSpace(CV_cvodemem, &CV_iout[12], &CV_iout[13]); /* LENRWLS,LENIWLS */
413 CVDiagGetLastFlag(CV_cvodemem, &CV_iout[14]); /* LSTF */
414 CVDiagGetNumRhsEvals(CV_cvodemem, &CV_iout[15]); /* NFELS */
415 }
416 }
417
418 /***************************************************************************/
419
FCV_DKY(realtype * t,int * k,realtype * dky,int * ier)420 void FCV_DKY (realtype *t, int *k, realtype *dky, int *ier)
421 {
422 /*
423 t is the t value where output is desired
424 k is the derivative order
425 F2C_CVODE_vec is the N_Vector containing the solution derivative on return
426 */
427
428 realtype *f2c_data = N_VGetArrayPointer(F2C_CVODE_vec);
429 N_VSetArrayPointer(dky, F2C_CVODE_vec);
430
431 *ier = 0;
432 *ier = CVodeGetDky(CV_cvodemem, *t, *k, F2C_CVODE_vec);
433
434 N_VSetArrayPointer(f2c_data, F2C_CVODE_vec);
435
436 }
437
438 /*************************************************/
439
FCV_GETERRWEIGHTS(realtype * eweight,int * ier)440 void FCV_GETERRWEIGHTS(realtype *eweight, int *ier)
441 {
442 /* Attach user data to vector */
443 realtype *f2c_data = N_VGetArrayPointer(F2C_CVODE_vec);
444 N_VSetArrayPointer(eweight, F2C_CVODE_vec);
445
446 *ier = 0;
447 *ier = CVodeGetErrWeights(CV_cvodemem, F2C_CVODE_vec);
448
449 /* Reset data pointers */
450 N_VSetArrayPointer(f2c_data, F2C_CVODE_vec);
451
452 return;
453 }
454
455 /*************************************************/
456
FCV_GETESTLOCALERR(realtype * ele,int * ier)457 void FCV_GETESTLOCALERR(realtype *ele, int *ier)
458 {
459 /* Attach user data to vector */
460 realtype *f2c_data = N_VGetArrayPointer(F2C_CVODE_vec);
461 N_VSetArrayPointer(ele, F2C_CVODE_vec);
462
463 *ier = 0;
464 *ier = CVodeGetEstLocalErrors(CV_cvodemem, F2C_CVODE_vec);
465
466 /* Reset data pointers */
467 N_VSetArrayPointer(f2c_data, F2C_CVODE_vec);
468
469 return;
470 }
471
472 /***************************************************************************/
473
FCV_FREE()474 void FCV_FREE ()
475 {
476 CVodeMem cv_mem;
477
478 cv_mem = (CVodeMem) CV_cvodemem;
479
480 if (cv_mem->cv_lfree)
481 cv_mem->cv_lfree(cv_mem);
482 cv_mem->cv_lmem = NULL;
483
484 free(cv_mem->cv_user_data); cv_mem->cv_user_data = NULL;
485
486 CVodeFree(&CV_cvodemem);
487
488 N_VSetArrayPointer(NULL, F2C_CVODE_vec);
489 N_VDestroy(F2C_CVODE_vec);
490 if (F2C_CVODE_matrix)
491 SUNMatDestroy(F2C_CVODE_matrix);
492 if (F2C_CVODE_linsol)
493 SUNLinSolFree(F2C_CVODE_linsol);
494 if (F2C_CVODE_nonlinsol)
495 SUNNonlinSolFree(F2C_CVODE_nonlinsol);
496 return;
497 }
498
499 /***************************************************************************/
500
501 /*
502 * C function CVf to interface between CVODE and a Fortran subroutine FCVFUN.
503 * Addresses of t, y, and ydot are passed to CVFUN, using the
504 * routine N_VGetArrayPointer from the NVECTOR module.
505 * Auxiliary data is assumed to be communicated by Common.
506 */
507
FCVf(realtype t,N_Vector y,N_Vector ydot,void * user_data)508 int FCVf(realtype t, N_Vector y, N_Vector ydot, void *user_data)
509 {
510 int ier;
511 realtype *ydata, *dydata;
512 FCVUserData CV_userdata;
513
514 ydata = N_VGetArrayPointer(y);
515 dydata = N_VGetArrayPointer(ydot);
516
517 CV_userdata = (FCVUserData) user_data;
518
519 FCV_FUN(&t, ydata, dydata, CV_userdata->ipar, CV_userdata->rpar, &ier);
520
521 return(ier);
522 }
523
524 /* Fortran interface to C routine CVodeSetNonlinearSolver; see
525 fcvode.h for further details */
FCV_NLSINIT(int * ier)526 void FCV_NLSINIT(int *ier) {
527 if ( (CV_cvodemem == NULL) || (F2C_CVODE_nonlinsol == NULL) ) {
528 *ier = -1;
529 return;
530 }
531 if (((CVodeMem) CV_cvodemem)->cv_lsolve == NULL) {
532 FCVNullMatrix();
533 FCVNullLinsol();
534 }
535
536 *ier = CVodeSetNonlinearSolver(CV_cvodemem, F2C_CVODE_nonlinsol);
537 return;
538 }
539