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 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 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 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 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 *); 74 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 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 \cdot{u} = -F(t,u) + G(t,u) 146 147 .seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSType` 148 M*/ 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