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 { 13 TS ts = (TS) IPAR; 14 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 15 DM dm; 16 DMTS tsdm; 17 TSIFunction 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 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 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 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 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(0); 120 } 121 122 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(0); 131 } 132 133 /*MC 134 TSRADAU5 - ODE solver using the RADAU5 package 135 136 Notes: 137 This uses its own nonlinear solver and dense matrix direct solver so PETSc SNES and KSP options do not apply. 138 Uses its own time-step adaptivity (but uses the TS rtol and atol, and initial timestep) 139 Uses its own memory for the dense matrix storage and factorization 140 Can only handle ODEs of the form \cdot{u} = -F(t,u) + G(t,u) 141 142 Level: beginner 143 144 .seealso: TSCreate(), TS, TSSetType() 145 146 M*/ 147 PETSC_EXTERN PetscErrorCode TSCreate_Radau5(TS ts) 148 { 149 TS_Radau5 *cvode; 150 151 PetscFunctionBegin; 152 ts->ops->destroy = TSDestroy_Radau5; 153 ts->ops->solve = TSSolve_Radau5; 154 ts->default_adapt_type = TSADAPTNONE; 155 156 PetscCall(PetscNewLog(ts,&cvode)); 157 ts->data = (void*)cvode; 158 PetscFunctionReturn(0); 159 } 160