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