Actual source code: radau5.c
1: /*
2: Provides a PETSc interface to RADAU5 solver.
4: */
5: #include <petsc/private/tsimpl.h>
7: typedef struct {
8: Vec work, workf;
9: } TS_Radau5;
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;
19: PetscCallAbort(PETSC_COMM_SELF, VecPlaceArray(cvode->work, Y));
20: PetscCallAbort(PETSC_COMM_SELF, VecPlaceArray(cvode->workf, F));
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;
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: }
38: PetscCallAbort(PETSC_COMM_SELF, VecResetArray(cvode->work));
39: PetscCallAbort(PETSC_COMM_SELF, VecResetArray(cvode->workf));
40: }
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;
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: }
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;
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: }
73: static 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 *);
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;
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));
88: LWORK = 4 * ND * ND + 12 * ND + 20;
89: LIWORK = 3 * ND + 20;
91: PetscCall(PetscCalloc2(LWORK, &WORK, LIWORK, &IWORK));
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;
114: /* output MUJAC MLMAS IDID; currently all ignored */
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);
118: PetscCall(PetscFree2(WORK, IWORK));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode TSDestroy_Radau5(TS ts)
123: {
124: TS_Radau5 *cvode = (TS_Radau5 *)ts->data;
126: PetscFunctionBegin;
127: PetscCall(VecDestroy(&cvode->work));
128: PetscCall(VecDestroy(&cvode->workf));
129: PetscCall(PetscFree(ts->data));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /*MC
134: TSRADAU5 - ODE solver using the external RADAU5 package, requires ./configure --download-radau5
136: Level: beginner
138: Notes:
139: This uses its own nonlinear solver and dense matrix direct solver so PETSc `SNES` and `KSP` options do not apply.
141: Uses its own time-step adaptivity (but uses the TS rtol and atol, and initial timestep)
143: Uses its own memory for the dense matrix storage and factorization
145: Can only handle ODEs of the form \cdot{u} = -F(t,u) + G(t,u)
147: .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSType`
148: M*/
149: PETSC_EXTERN PetscErrorCode TSCreate_Radau5(TS ts)
150: {
151: TS_Radau5 *cvode;
153: PetscFunctionBegin;
154: ts->ops->destroy = TSDestroy_Radau5;
155: ts->ops->solve = TSSolve_Radau5;
156: ts->default_adapt_type = TSADAPTNONE;
158: PetscCall(PetscNew(&cvode));
159: ts->data = (void *)cvode;
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }