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