1: /*
  2:     Provides a PETSc interface to RADAU5 solver.

  4: */
  5: #include <petsc/private/tsimpl.h>                /*I   "petscts.h"   I*/

  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;

 20:   VecPlaceArray(cvode->work,Y);CHKERRABORT(PETSC_COMM_SELF,ierr);
 21:   VecPlaceArray(cvode->workf,F);CHKERRABORT(PETSC_COMM_SELF,ierr);

 23:   /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */
 24:   TSGetDM(ts,&dm);CHKERRABORT(PETSC_COMM_SELF,ierr);
 25:   DMGetDMTS(dm,&tsdm);CHKERRABORT(PETSC_COMM_SELF,ierr);
 26:   DMTSGetIFunction(dm,&ifunction,NULL);CHKERRABORT(PETSC_COMM_SELF,ierr);
 27:   if (!ifunction) {
 28:     TSComputeRHSFunction(ts,*X,cvode->work,cvode->workf);CHKERRABORT(PETSC_COMM_SELF,ierr);
 29:   } else {       /* If rhsfunction is also set, this computes both parts and scale them to the right hand side */
 30:     Vec yydot;

 32:     VecDuplicate(cvode->work,&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
 33:     VecZeroEntries(yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
 34:     TSComputeIFunction(ts,*X,cvode->work,yydot,cvode->workf,PETSC_FALSE);CHKERRABORT(PETSC_COMM_SELF,ierr);
 35:     VecScale(cvode->workf,-1.);CHKERRABORT(PETSC_COMM_SELF,ierr);
 36:     VecDestroy(&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
 37:   }

 39:   VecResetArray(cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr);
 40:   VecResetArray(cvode->workf);CHKERRABORT(PETSC_COMM_SELF,ierr);
 41: }

 43: void JVPOL(PetscInt *N,PetscScalar *X,PetscScalar *Y,PetscScalar *DFY,int *LDFY,PetscScalar *RPAR,void *IPAR)
 44: {
 45:   TS             ts = (TS) IPAR;
 46:   TS_Radau5      *cvode = (TS_Radau5*)ts->data;
 47:   Vec            yydot;
 48:   Mat            mat;
 49:   PetscInt       n;

 52:   VecPlaceArray(cvode->work,Y);CHKERRABORT(PETSC_COMM_SELF,ierr);
 53:   VecDuplicate(cvode->work,&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
 54:   VecGetSize(yydot,&n);CHKERRABORT(PETSC_COMM_SELF,ierr);
 55:   MatCreateSeqDense(PETSC_COMM_SELF,n,n,DFY,&mat);CHKERRABORT(PETSC_COMM_SELF,ierr);
 56:   VecZeroEntries(yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
 57:   TSComputeIJacobian(ts,*X,cvode->work,yydot,0,mat,mat,PETSC_FALSE);CHKERRABORT(PETSC_COMM_SELF,ierr);
 58:   MatScale(mat,-1.0);CHKERRABORT(PETSC_COMM_SELF,ierr);
 59:   MatDestroy(&mat);CHKERRABORT(PETSC_COMM_SELF,ierr);
 60:   VecDestroy(&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
 61:   VecResetArray(cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr);
 62: }

 64: void SOLOUT(int *NR,double *XOLD,double *X, double *Y,double *CONT,double *LRC,int *N,double *RPAR,void *IPAR,int *IRTRN)
 65: {
 66:   TS             ts = (TS) IPAR;
 67:   TS_Radau5      *cvode = (TS_Radau5*)ts->data;

 70:   VecPlaceArray(cvode->work,Y);CHKERRABORT(PETSC_COMM_SELF,ierr);
 71:   ts->time_step = *X - *XOLD;
 72:   TSMonitor(ts,*NR-1,*X,cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr);
 73:   VecResetArray(cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr);
 74: }


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

 79: PetscErrorCode TSSolve_Radau5(TS ts)
 80: {
 81:   TS_Radau5      *cvode = (TS_Radau5*)ts->data;
 83:   PetscScalar    *Y,*WORK,X,XEND,RTOL,ATOL,H,RPAR;
 84:   PetscInt       ND,*IWORK,LWORK,LIWORK,MUJAC,MLMAS,MUMAS,IDID,ITOL;
 85:   int            IJAC,MLJAC,IMAS,IOUT;

 88:   VecGetArray(ts->vec_sol,&Y);
 89:   VecGetSize(ts->vec_sol,&ND);
 90:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,ND,NULL,&cvode->work);
 91:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,ND,NULL,&cvode->workf);

 93:   LWORK  = 4*ND*ND+12*ND+20;
 94:   LIWORK = 3*ND+20;

 96:   PetscCalloc1(LWORK,&WORK);
 97:   PetscCalloc1(LIWORK,&IWORK);

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

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

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

124:   PetscFree(WORK);
125:   PetscFree(IWORK);
126:   return(0);
127: }

129: PetscErrorCode TSDestroy_Radau5(TS ts)
130: {
131:   TS_Radau5      *cvode = (TS_Radau5*)ts->data;

135:   VecDestroy(&cvode->work);
136:   VecDestroy(&cvode->workf);
137:   PetscFree(ts->data);
138:   return(0);
139: }

141: /*MC
142:       TSRADAU5 - ODE solver using the RADAU5 package

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

149:     Level: beginner

151: .seealso:  TSCreate(), TS, TSSetType()

153: M*/
154: PETSC_EXTERN PetscErrorCode TSCreate_Radau5(TS ts)
155: {
156:   TS_Radau5      *cvode;

160:   ts->ops->destroy        = TSDestroy_Radau5;
161:   ts->ops->solve          = TSSolve_Radau5;
162:   ts->default_adapt_type  = TSADAPTNONE;

164:   PetscNewLog(ts,&cvode);
165:   ts->data = (void*)cvode;
166:   return(0);
167: }