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 #ifdef foo 12 /* 13 TSFunction_Radau5 - routine that we provide to Radau5 that applies the right hand side. 14 */ 15 int TSFunction_Radau5(realtype t,N_Vector y,N_Vector ydot,void *ctx) 16 { 17 TS ts = (TS) ctx; 18 DM dm; 19 DMTS tsdm; 20 TSIFunction ifunction; 21 MPI_Comm comm; 22 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 23 Vec yy = cvode->w1,yyd = cvode->w2,yydot = cvode->ydot; 24 PetscScalar *y_data,*ydot_data; 25 PetscErrorCode ierr; 26 27 PetscFunctionBegin; 28 ierr = PetscObjectGetComm((PetscObject)ts,&comm);CHKERRQ(ierr); 29 /* Make the PETSc work vectors yy and yyd point to the arrays in the RADAU5 vectors y and ydot respectively*/ 30 y_data = (PetscScalar*) N_VGetArrayPointer(y); 31 ydot_data = (PetscScalar*) N_VGetArrayPointer(ydot); 32 ierr = VecPlaceArray(yy,y_data);CHKERRABORT(comm,ierr); 33 ierr = VecPlaceArray(yyd,ydot_data);CHKERRABORT(comm,ierr); 34 35 /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */ 36 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 37 ierr = DMGetDMTS(dm,&tsdm);CHKERRQ(ierr); 38 ierr = DMTSGetIFunction(dm,&ifunction,NULL);CHKERRQ(ierr); 39 if (!ifunction) { 40 ierr = TSComputeRHSFunction(ts,t,yy,yyd);CHKERRQ(ierr); 41 } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */ 42 ierr = VecZeroEntries(yydot);CHKERRQ(ierr); 43 ierr = TSComputeIFunction(ts,t,yy,yydot,yyd,PETSC_FALSE);CHKERRABORT(comm,ierr); 44 ierr = VecScale(yyd,-1.);CHKERRQ(ierr); 45 } 46 ierr = VecResetArray(yy);CHKERRABORT(comm,ierr); 47 ierr = VecResetArray(yyd);CHKERRABORT(comm,ierr); 48 PetscFunctionReturn(0); 49 } 50 #endif 51 52 void FVPOL(int *N,double *X,double *Y,double *F,double *RPAR,void *IPAR) 53 { 54 TS ts = (TS) IPAR; 55 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 56 DM dm; 57 DMTS tsdm; 58 TSIFunction ifunction; 59 60 VecPlaceArray(cvode->work,Y); 61 VecPlaceArray(cvode->workf,F); 62 63 /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */ 64 TSGetDM(ts,&dm); 65 DMGetDMTS(dm,&tsdm); 66 DMTSGetIFunction(dm,&ifunction,NULL); 67 if (!ifunction) { 68 TSComputeRHSFunction(ts,*X,cvode->work,cvode->workf); 69 } else { /* If rhsfunction is also set, this computes both parts and shifts them to the right */ 70 Vec yydot; 71 72 VecDuplicate(cvode->work,&yydot); 73 VecZeroEntries(yydot); 74 TSComputeIFunction(ts,*X,cvode->work,yydot,cvode->workf,PETSC_FALSE); 75 VecScale(cvode->workf,-1.); 76 VecDestroy(&yydot); 77 } 78 79 VecResetArray(cvode->work); 80 VecResetArray(cvode->workf); 81 } 82 83 void JVPOL(PetscInt *N,PetscScalar *X,PetscScalar *Y,PetscScalar *DFY,int *LDFY,PetscScalar *RPAR,void *IPAR) 84 { 85 TS ts = (TS) IPAR; 86 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 87 Vec yydot; 88 Mat mat; 89 PetscInt n; 90 91 VecPlaceArray(cvode->work,Y); 92 VecDuplicate(cvode->work,&yydot); 93 VecGetSize(yydot,&n); 94 MatCreateSeqDense(PETSC_COMM_SELF,n,n,DFY,&mat); 95 VecZeroEntries(yydot); 96 TSComputeIJacobian(ts,*X,cvode->work,yydot,0,mat,mat,PETSC_FALSE); 97 MatScale(mat,-1.0); 98 MatDestroy(&mat); 99 VecDestroy(&yydot); 100 VecResetArray(cvode->work); 101 } 102 103 void SOLOUT(int *NR,double *XOLD,double *X, double *Y,double *CONT,double *LRC,int *N,double *RPAR,void *IPAR,int *IRTRN) 104 { 105 TS ts = (TS) IPAR; 106 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 107 108 VecPlaceArray(cvode->work,Y); 109 ts->time_step = *X - *XOLD; 110 TSMonitor(ts,*NR-1,*X,cvode->work); 111 VecResetArray(cvode->work); 112 } 113 114 115 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*); 116 117 PetscErrorCode TSSolve_Radau5(TS ts) 118 { 119 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 120 PetscErrorCode ierr; 121 PetscScalar *Y,*WORK,X,XEND,RTOL,ATOL,H,RPAR; 122 PetscInt ND,*IWORK,LWORK,LIWORK,MUJAC,MLMAS,MUMAS,IDID,ITOL; 123 int IJAC,MLJAC,IMAS,IOUT; 124 125 PetscFunctionBegin; 126 ierr = VecGetArray(ts->vec_sol,&Y);CHKERRQ(ierr); 127 ierr = VecGetSize(ts->vec_sol,&ND);CHKERRQ(ierr); 128 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,ND,NULL,&cvode->work);CHKERRQ(ierr); 129 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,ND,NULL,&cvode->workf);CHKERRQ(ierr); 130 131 LWORK = 4*ND*ND+12*ND+20; 132 LIWORK = 3*ND+20; 133 134 ierr = PetscCalloc1(LWORK,&WORK);CHKERRQ(ierr); 135 ierr = PetscCalloc1(LIWORK,&IWORK);CHKERRQ(ierr); 136 137 /* C --- PARAMETER IN THE DIFFERENTIAL EQUATION */ 138 RPAR=1.0e-6; 139 /* C --- COMPUTE THE JACOBIAN ANALYTICALLY */ 140 IJAC=1; 141 /* C --- JACOBIAN IS A FULL MATRIX */ 142 MLJAC=ND; 143 /* C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM*/ 144 IMAS=0; 145 /* C --- OUTPUT ROUTINE IS USED DURING INTEGRATION*/ 146 IOUT=1; 147 /* C --- INITIAL VALUES*/ 148 X = ts->ptime; 149 /* C --- ENDPOINT OF INTEGRATION */ 150 XEND = ts->max_time; 151 /* C --- REQUIRED TOLERANCE */ 152 RTOL = ts->rtol; 153 ATOL = ts->atol; 154 ITOL=0; 155 /* C --- INITIAL STEP SIZE */ 156 H = ts->time_step; 157 158 /* output MUJAC MLMAS IDID */ 159 160 CHKMEMQ; 161 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); 162 CHKMEMQ; 163 164 ierr = PetscFree(WORK);CHKERRQ(ierr); 165 ierr = PetscFree(IWORK);CHKERRQ(ierr); 166 PetscFunctionReturn(0); 167 } 168 169 PetscErrorCode TSDestroy_Radau5(TS ts) 170 { 171 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 172 PetscErrorCode ierr; 173 174 PetscFunctionBegin; 175 ierr = VecDestroy(&cvode->work);CHKERRQ(ierr); 176 ierr = VecDestroy(&cvode->workf);CHKERRQ(ierr); 177 ierr = PetscFree(ts->data);CHKERRQ(ierr); 178 PetscFunctionReturn(0); 179 } 180 181 PetscErrorCode TSSetFromOptions_Radau5(PetscOptionItems *PetscOptionsObject,TS ts) 182 { 183 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 184 PetscErrorCode ierr; 185 186 PetscFunctionBegin; 187 ierr = PetscOptionsHead(PetscOptionsObject,"RADAU5 ODE solver options");CHKERRQ(ierr); 188 ierr = PetscOptionsTail();CHKERRQ(ierr); 189 PetscFunctionReturn(0); 190 } 191 192 PetscErrorCode TSView_Radau5(TS ts,PetscViewer viewer) 193 { 194 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 195 PetscErrorCode ierr; 196 197 PetscFunctionBegin; 198 PetscFunctionReturn(0); 199 } 200 201 /*MC 202 TSRADAU5 - ODE solver using the RADAU5 package 203 204 Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply 205 206 Level: beginner 207 208 .seealso: TSCreate(), TS, TSSetType(), TSSetExactFinalTime() 209 210 M*/ 211 PETSC_EXTERN PetscErrorCode TSCreate_Radau5(TS ts) 212 { 213 TS_Radau5 *cvode; 214 PetscErrorCode ierr; 215 216 PetscFunctionBegin; 217 ts->ops->destroy = TSDestroy_Radau5; 218 ts->ops->view = TSView_Radau5; 219 ts->ops->solve = TSSolve_Radau5; 220 ts->ops->setfromoptions = TSSetFromOptions_Radau5; 221 ts->default_adapt_type = TSADAPTNONE; 222 223 ierr = PetscNewLog(ts,&cvode);CHKERRQ(ierr); 224 ts->data = (void*)cvode; 225 PetscFunctionReturn(0); 226 } 227