1 /*
2 Provides a PETSc interface to RADAU5 solver.
3
4 */
5 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
6
7 typedef struct {
8 Vec work, workf;
9 } TS_Radau5;
10
FVPOL(int * N,double * X,double * Y,double * F,double * RPAR,void * IPAR)11 static void FVPOL(int *N, double *X, double *Y, double *F, double *RPAR, void *IPAR)
12 {
13 TS ts = (TS)IPAR;
14 TS_Radau5 *cvode = (TS_Radau5 *)ts->data;
15 DM dm;
16 DMTS tsdm;
17 TSIFunctionFn *ifunction;
18
19 PetscCallAbort(PETSC_COMM_SELF, VecPlaceArray(cvode->work, Y));
20 PetscCallAbort(PETSC_COMM_SELF, VecPlaceArray(cvode->workf, F));
21
22 /* Now compute the right-hand side function, via IFunction unless only the more efficient RHSFunction is set */
23 PetscCallAbort(PETSC_COMM_SELF, TSGetDM(ts, &dm));
24 PetscCallAbort(PETSC_COMM_SELF, DMGetDMTS(dm, &tsdm));
25 PetscCallAbort(PETSC_COMM_SELF, DMTSGetIFunction(dm, &ifunction, NULL));
26 if (!ifunction) {
27 PetscCallAbort(PETSC_COMM_SELF, TSComputeRHSFunction(ts, *X, cvode->work, cvode->workf));
28 } else { /* If rhsfunction is also set, this computes both parts and scale them to the right-hand side */
29 Vec yydot;
30
31 PetscCallAbort(PETSC_COMM_SELF, VecDuplicate(cvode->work, &yydot));
32 PetscCallAbort(PETSC_COMM_SELF, VecZeroEntries(yydot));
33 PetscCallAbort(PETSC_COMM_SELF, TSComputeIFunction(ts, *X, cvode->work, yydot, cvode->workf, PETSC_FALSE));
34 PetscCallAbort(PETSC_COMM_SELF, VecScale(cvode->workf, -1.));
35 PetscCallAbort(PETSC_COMM_SELF, VecDestroy(&yydot));
36 }
37
38 PetscCallAbort(PETSC_COMM_SELF, VecResetArray(cvode->work));
39 PetscCallAbort(PETSC_COMM_SELF, VecResetArray(cvode->workf));
40 }
41
JVPOL(PetscInt * N,PetscScalar * X,PetscScalar * Y,PetscScalar * DFY,int * LDFY,PetscScalar * RPAR,void * IPAR)42 static void JVPOL(PetscInt *N, PetscScalar *X, PetscScalar *Y, PetscScalar *DFY, int *LDFY, PetscScalar *RPAR, void *IPAR)
43 {
44 TS ts = (TS)IPAR;
45 TS_Radau5 *cvode = (TS_Radau5 *)ts->data;
46 Vec yydot;
47 Mat mat;
48 PetscInt n;
49
50 PetscCallAbort(PETSC_COMM_SELF, VecPlaceArray(cvode->work, Y));
51 PetscCallAbort(PETSC_COMM_SELF, VecDuplicate(cvode->work, &yydot));
52 PetscCallAbort(PETSC_COMM_SELF, VecGetSize(yydot, &n));
53 PetscCallAbort(PETSC_COMM_SELF, MatCreateSeqDense(PETSC_COMM_SELF, n, n, DFY, &mat));
54 PetscCallAbort(PETSC_COMM_SELF, VecZeroEntries(yydot));
55 PetscCallAbort(PETSC_COMM_SELF, TSComputeIJacobian(ts, *X, cvode->work, yydot, 0, mat, mat, PETSC_FALSE));
56 PetscCallAbort(PETSC_COMM_SELF, MatScale(mat, -1.0));
57 PetscCallAbort(PETSC_COMM_SELF, MatDestroy(&mat));
58 PetscCallAbort(PETSC_COMM_SELF, VecDestroy(&yydot));
59 PetscCallAbort(PETSC_COMM_SELF, VecResetArray(cvode->work));
60 }
61
SOLOUT(int * NR,double * XOLD,double * X,double * Y,double * CONT,double * LRC,int * N,double * RPAR,void * IPAR,int * IRTRN)62 static void SOLOUT(int *NR, double *XOLD, double *X, double *Y, double *CONT, double *LRC, int *N, double *RPAR, void *IPAR, int *IRTRN)
63 {
64 TS ts = (TS)IPAR;
65 TS_Radau5 *cvode = (TS_Radau5 *)ts->data;
66
67 PetscCallAbort(PETSC_COMM_SELF, VecPlaceArray(cvode->work, Y));
68 ts->time_step = *X - *XOLD;
69 PetscCallAbort(PETSC_COMM_SELF, TSMonitor(ts, *NR - 1, *X, cvode->work));
70 PetscCallAbort(PETSC_COMM_SELF, VecResetArray(cvode->work));
71 }
72
73 extern void radau5_(int *, void *, double *, double *, double *, double *, double *, double *, int *, void *, int *, int *, int *, void *, int *, int *, int *, void *, int *, double *, int *, int *, int *, double *, void *, int *);
74
TSSolve_Radau5(TS ts)75 static PetscErrorCode TSSolve_Radau5(TS ts)
76 {
77 TS_Radau5 *cvode = (TS_Radau5 *)ts->data;
78 PetscScalar *Y, *WORK, X, XEND, RTOL, ATOL, H, RPAR;
79 PetscInt ND, *IWORK, LWORK, LIWORK, MUJAC, MLMAS, MUMAS, IDID, ITOL;
80 int IJAC, MLJAC, IMAS, IOUT;
81
82 PetscFunctionBegin;
83 PetscCall(VecGetArray(ts->vec_sol, &Y));
84 PetscCall(VecGetSize(ts->vec_sol, &ND));
85 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, ND, NULL, &cvode->work));
86 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, ND, NULL, &cvode->workf));
87
88 LWORK = 4 * ND * ND + 12 * ND + 20;
89 LIWORK = 3 * ND + 20;
90
91 PetscCall(PetscCalloc2(LWORK, &WORK, LIWORK, &IWORK));
92
93 /* C --- PARAMETER IN THE DIFFERENTIAL EQUATION */
94 RPAR = 1.0e-6;
95 /* C --- COMPUTE THE JACOBIAN ANALYTICALLY */
96 IJAC = 1;
97 /* C --- JACOBIAN IS A FULL MATRIX */
98 MLJAC = ND;
99 /* C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM*/
100 IMAS = 0;
101 /* C --- OUTPUT ROUTINE IS USED DURING INTEGRATION*/
102 IOUT = 1;
103 /* C --- INITIAL VALUES*/
104 X = ts->ptime;
105 /* C --- ENDPOINT OF INTEGRATION */
106 XEND = ts->max_time;
107 /* C --- REQUIRED TOLERANCE */
108 RTOL = ts->rtol;
109 ATOL = ts->atol;
110 ITOL = 0;
111 /* C --- INITIAL STEP SIZE */
112 H = ts->time_step;
113
114 /* output MUJAC MLMAS IDID; currently all ignored */
115
116 radau5_(&ND, FVPOL, &X, Y, &XEND, &H, &RTOL, &ATOL, &ITOL, JVPOL, &IJAC, &MLJAC, &MUJAC, FVPOL, &IMAS, &MLMAS, &MUMAS, SOLOUT, &IOUT, WORK, &LWORK, IWORK, &LIWORK, &RPAR, (void *)ts, &IDID);
117
118 PetscCall(PetscFree2(WORK, IWORK));
119 PetscFunctionReturn(PETSC_SUCCESS);
120 }
121
TSDestroy_Radau5(TS ts)122 static PetscErrorCode TSDestroy_Radau5(TS ts)
123 {
124 TS_Radau5 *cvode = (TS_Radau5 *)ts->data;
125
126 PetscFunctionBegin;
127 PetscCall(VecDestroy(&cvode->work));
128 PetscCall(VecDestroy(&cvode->workf));
129 PetscCall(PetscFree(ts->data));
130 PetscFunctionReturn(PETSC_SUCCESS);
131 }
132
133 /*MC
134 TSRADAU5 - ODE solver using the external RADAU5 package, requires ./configure --download-radau5
135
136 Level: beginner
137
138 Notes:
139 This uses its own nonlinear solver and dense matrix direct solver so PETSc `SNES` and `KSP` options do not apply.
140
141 Uses its own time-step adaptivity (but uses the TS rtol and atol, and initial timestep)
142
143 Uses its own memory for the dense matrix storage and factorization
144
145 Can only handle ODEs of the form $ \dot{u} = -F(t,u) + G(t,u) $
146
147 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSType`
148 M*/
TSCreate_Radau5(TS ts)149 PETSC_EXTERN PetscErrorCode TSCreate_Radau5(TS ts)
150 {
151 TS_Radau5 *cvode;
152
153 PetscFunctionBegin;
154 ts->ops->destroy = TSDestroy_Radau5;
155 ts->ops->solve = TSSolve_Radau5;
156 ts->default_adapt_type = TSADAPTNONE;
157
158 PetscCall(PetscNew(&cvode));
159 ts->data = (void *)cvode;
160 PetscFunctionReturn(PETSC_SUCCESS);
161 }
162