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: 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 *);

     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 $ \dot{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: }