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: 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:   TSIFunction ifunction;

 19:   PETSC_COMM_SELF, VecPlaceArray(cvode->work, Y);
 20:   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:   PETSC_COMM_SELF, TSGetDM(ts, &dm);
 24:   PETSC_COMM_SELF, DMGetDMTS(dm, &tsdm);
 25:   PETSC_COMM_SELF, DMTSGetIFunction(dm, &ifunction, NULL);
 26:   if (!ifunction) {
 27:     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:     PETSC_COMM_SELF, VecDuplicate(cvode->work, &yydot);
 32:     PETSC_COMM_SELF, VecZeroEntries(yydot);
 33:     PETSC_COMM_SELF, TSComputeIFunction(ts, *X, cvode->work, yydot, cvode->workf, PETSC_FALSE);
 34:     PETSC_COMM_SELF, VecScale(cvode->workf, -1.);
 35:     PETSC_COMM_SELF, VecDestroy(&yydot);
 36:   }

 38:   PETSC_COMM_SELF, VecResetArray(cvode->work);
 39:   PETSC_COMM_SELF, VecResetArray(cvode->workf);
 40: }

 42: 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:   PETSC_COMM_SELF, VecPlaceArray(cvode->work, Y);
 51:   PETSC_COMM_SELF, VecDuplicate(cvode->work, &yydot);
 52:   PETSC_COMM_SELF, VecGetSize(yydot, &n);
 53:   PETSC_COMM_SELF, MatCreateSeqDense(PETSC_COMM_SELF, n, n, DFY, &mat);
 54:   PETSC_COMM_SELF, VecZeroEntries(yydot);
 55:   PETSC_COMM_SELF, TSComputeIJacobian(ts, *X, cvode->work, yydot, 0, mat, mat, PETSC_FALSE);
 56:   PETSC_COMM_SELF, MatScale(mat, -1.0);
 57:   PETSC_COMM_SELF, MatDestroy(&mat);
 58:   PETSC_COMM_SELF, VecDestroy(&yydot);
 59:   PETSC_COMM_SELF, VecResetArray(cvode->work);
 60: }

 62: 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:   PETSC_COMM_SELF, VecPlaceArray(cvode->work, Y);
 68:   ts->time_step = *X - *XOLD;
 69:   PETSC_COMM_SELF, TSMonitor(ts, *NR - 1, *X, cvode->work);
 70:   PETSC_COMM_SELF, VecResetArray(cvode->work);
 71: }

 73: 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: 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:   VecGetArray(ts->vec_sol, &Y);
 83:   VecGetSize(ts->vec_sol, &ND);
 84:   VecCreateSeqWithArray(PETSC_COMM_SELF, 1, ND, NULL, &cvode->work);
 85:   VecCreateSeqWithArray(PETSC_COMM_SELF, 1, ND, NULL, &cvode->workf);

 87:   LWORK  = 4 * ND * ND + 12 * ND + 20;
 88:   LIWORK = 3 * ND + 20;

 90:   PetscCalloc2(LWORK, &WORK, LIWORK, &IWORK);

 92:   /* C --- PARAMETER IN THE DIFFERENTIAL EQUATION */
 93:   RPAR = 1.0e-6;
 94:   /* C --- COMPUTE THE JACOBIAN ANALYTICALLY */
 95:   IJAC = 1;
 96:   /* C --- JACOBIAN IS A FULL MATRIX */
 97:   MLJAC = ND;
 98:   /* C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM*/
 99:   IMAS = 0;
100:   /* C --- OUTPUT ROUTINE IS USED DURING INTEGRATION*/
101:   IOUT = 1;
102:   /* C --- INITIAL VALUES*/
103:   X = ts->ptime;
104:   /* C --- ENDPOINT OF INTEGRATION */
105:   XEND = ts->max_time;
106:   /* C --- REQUIRED TOLERANCE */
107:   RTOL = ts->rtol;
108:   ATOL = ts->atol;
109:   ITOL = 0;
110:   /* C --- INITIAL STEP SIZE */
111:   H = ts->time_step;

113:   /* output MUJAC MLMAS IDID; currently all ignored */

115:   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:   PetscFree2(WORK, IWORK);
118:   return 0;
119: }

121: PetscErrorCode TSDestroy_Radau5(TS ts)
122: {
123:   TS_Radau5 *cvode = (TS_Radau5 *)ts->data;

125:   VecDestroy(&cvode->work);
126:   VecDestroy(&cvode->workf);
127:   PetscFree(ts->data);
128:   return 0;
129: }

131: /*MC
132:       TSRADAU5 - ODE solver using the RADAU5 package

134:     Notes:
135:     This uses its own nonlinear solver and dense matrix direct solver so PETSc SNES and KSP options do not apply.
136:            Uses its own time-step adaptivity (but uses the TS rtol and atol, and initial timestep)
137:            Uses its own memory for the dense matrix storage and factorization
138:            Can only handle ODEs of the form \cdot{u} = -F(t,u) + G(t,u)

140:     Level: beginner

142: .seealso: `TSCreate()`, `TS`, `TSSetType()`

144: M*/
145: PETSC_EXTERN PetscErrorCode TSCreate_Radau5(TS ts)
146: {
147:   TS_Radau5 *cvode;

149:   ts->ops->destroy       = TSDestroy_Radau5;
150:   ts->ops->solve         = TSSolve_Radau5;
151:   ts->default_adapt_type = TSADAPTNONE;

153:   PetscNew(&cvode);
154:   ts->data = (void *)cvode;
155:   return 0;
156: }