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 #ifdef foo 84 void FVPOL(int *N,double *X,double *Y,double *F,double *RPAR,void *IPAR) 85 { 86 F[0]=Y[1]; 87 F[1]=((1-Y[0]*Y[0])*Y[1]-Y[0])/(*RPAR); 88 } 89 #endif 90 91 void JVPOL(PetscInt *N,PetscScalar *X,PetscScalar *Y,PetscScalar *DFY,int *LDFY,PetscScalar *RPAR,void *IPAR) 92 { 93 DFY[0]=0.0; 94 DFY[1]=(-2.0*Y[0]*Y[1]-1.0)/(*RPAR); 95 DFY[2]=1.0; 96 DFY[3]=(1.0-Y[0]*Y[0])/(*RPAR); 97 } 98 99 void SOLOUT(int *NR,double *XOLD,double *X, double *Y,double *CONT,double *LRC,int *N,double *RPAR,void *IPAR,int *IRTRN) 100 { 101 TS ts = (TS) IPAR; 102 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 103 104 VecPlaceArray(cvode->work,Y); 105 ts->time_step = *X - *XOLD; 106 TSMonitor(ts,*NR-1,*X,cvode->work); 107 VecResetArray(cvode->work); 108 PetscPrintf(PETSC_COMM_SELF,"X = %g Y = %g %g NSTEP = %d dt = %g\n",*X,Y[0],Y[1],*NR-1,*X - *XOLD); 109 } 110 111 112 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*); 113 114 PetscErrorCode TSSolve_Radau5(TS ts) 115 { 116 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 117 PetscErrorCode ierr; 118 PetscScalar *Y,*WORK,X,XEND,RTOL,ATOL,H,RPAR; 119 PetscInt ND,*IWORK,LWORK,LIWORK,MUJAC,MLMAS,MUMAS,IDID,ITOL; 120 int IJAC,MLJAC,IMAS,IOUT; 121 122 PetscFunctionBegin; 123 ierr = VecGetArray(ts->vec_sol,&Y);CHKERRQ(ierr); 124 ierr = VecGetSize(ts->vec_sol,&ND);CHKERRQ(ierr); 125 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,ND,NULL,&cvode->work);CHKERRQ(ierr); 126 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,ND,NULL,&cvode->workf);CHKERRQ(ierr); 127 128 LWORK = 4*ND*ND+12*ND+20; 129 LIWORK = 3*ND+20; 130 131 ierr = PetscCalloc1(LWORK,&WORK);CHKERRQ(ierr); 132 ierr = PetscCalloc1(LIWORK,&IWORK);CHKERRQ(ierr); 133 134 /* C --- PARAMETER IN THE DIFFERENTIAL EQUATION */ 135 RPAR=1.0e-6; 136 /* C --- COMPUTE THE JACOBIAN ANALYTICALLY */ 137 IJAC=1; 138 /* C --- JACOBIAN IS A FULL MATRIX */ 139 MLJAC=ND; 140 /* C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM*/ 141 IMAS=0; 142 /* C --- OUTPUT ROUTINE IS USED DURING INTEGRATION*/ 143 IOUT=1; 144 /* C --- INITIAL VALUES*/ 145 X = ts->ptime; 146 /* C --- ENDPOINT OF INTEGRATION */ 147 XEND = ts->max_time; 148 /* C --- REQUIRED TOLERANCE */ 149 RTOL = ts->rtol; 150 ATOL = ts->atol; 151 ITOL=0; 152 /* C --- INITIAL STEP SIZE */ 153 H = ts->time_step; 154 155 /* output MUJAC MLMAS IDID */ 156 157 CHKMEMQ; 158 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); 159 CHKMEMQ; 160 161 ierr = PetscFree(WORK);CHKERRQ(ierr); 162 ierr = PetscFree(IWORK);CHKERRQ(ierr); 163 PetscFunctionReturn(0); 164 } 165 166 PetscErrorCode TSDestroy_Radau5(TS ts) 167 { 168 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 169 PetscErrorCode ierr; 170 171 PetscFunctionBegin; 172 ierr = VecDestroy(&cvode->work);CHKERRQ(ierr); 173 ierr = VecDestroy(&cvode->workf);CHKERRQ(ierr); 174 ierr = PetscFree(ts->data);CHKERRQ(ierr); 175 PetscFunctionReturn(0); 176 } 177 178 PetscErrorCode TSSetFromOptions_Radau5(PetscOptionItems *PetscOptionsObject,TS ts) 179 { 180 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 181 PetscErrorCode ierr; 182 183 PetscFunctionBegin; 184 ierr = PetscOptionsHead(PetscOptionsObject,"RADAU5 ODE solver options");CHKERRQ(ierr); 185 ierr = PetscOptionsTail();CHKERRQ(ierr); 186 PetscFunctionReturn(0); 187 } 188 189 PetscErrorCode TSView_Radau5(TS ts,PetscViewer viewer) 190 { 191 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 192 PetscErrorCode ierr; 193 194 PetscFunctionBegin; 195 PetscFunctionReturn(0); 196 } 197 198 /*MC 199 TSRADAU5 - ODE solver using the RADAU5 package 200 201 Notes: This uses its own nonlinear solver and Krylov method so PETSc SNES and KSP options do not apply 202 203 Level: beginner 204 205 .seealso: TSCreate(), TS, TSSetType(), TSSetExactFinalTime() 206 207 M*/ 208 PETSC_EXTERN PetscErrorCode TSCreate_Radau5(TS ts) 209 { 210 TS_Radau5 *cvode; 211 PetscErrorCode ierr; 212 213 PetscFunctionBegin; 214 ts->ops->destroy = TSDestroy_Radau5; 215 ts->ops->view = TSView_Radau5; 216 ts->ops->solve = TSSolve_Radau5; 217 ts->ops->setfromoptions = TSSetFromOptions_Radau5; 218 ts->default_adapt_type = TSADAPTNONE; 219 220 ierr = PetscNewLog(ts,&cvode);CHKERRQ(ierr); 221 ts->data = (void*)cvode; 222 PetscFunctionReturn(0); 223 } 224