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