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 PetscErrorCode ierr; 19 20 ierr = VecPlaceArray(cvode->work,Y);CHKERRABORT(PETSC_COMM_SELF,ierr); 21 ierr = VecPlaceArray(cvode->workf,F);CHKERRABORT(PETSC_COMM_SELF,ierr); 22 23 /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */ 24 ierr = TSGetDM(ts,&dm);CHKERRABORT(PETSC_COMM_SELF,ierr); 25 ierr = DMGetDMTS(dm,&tsdm);CHKERRABORT(PETSC_COMM_SELF,ierr); 26 ierr = DMTSGetIFunction(dm,&ifunction,NULL);CHKERRABORT(PETSC_COMM_SELF,ierr); 27 if (!ifunction) { 28 ierr = TSComputeRHSFunction(ts,*X,cvode->work,cvode->workf);CHKERRABORT(PETSC_COMM_SELF,ierr); 29 } else { /* If rhsfunction is also set, this computes both parts and scale them to the right hand side */ 30 Vec yydot; 31 32 ierr = VecDuplicate(cvode->work,&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr); 33 ierr = VecZeroEntries(yydot);CHKERRABORT(PETSC_COMM_SELF,ierr); 34 ierr = TSComputeIFunction(ts,*X,cvode->work,yydot,cvode->workf,PETSC_FALSE);CHKERRABORT(PETSC_COMM_SELF,ierr); 35 ierr = VecScale(cvode->workf,-1.);CHKERRABORT(PETSC_COMM_SELF,ierr); 36 ierr = VecDestroy(&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr); 37 } 38 39 ierr = VecResetArray(cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr); 40 ierr = VecResetArray(cvode->workf);CHKERRABORT(PETSC_COMM_SELF,ierr); 41 } 42 43 void JVPOL(PetscInt *N,PetscScalar *X,PetscScalar *Y,PetscScalar *DFY,int *LDFY,PetscScalar *RPAR,void *IPAR) 44 { 45 TS ts = (TS) IPAR; 46 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 47 Vec yydot; 48 Mat mat; 49 PetscInt n; 50 PetscErrorCode ierr; 51 52 ierr = VecPlaceArray(cvode->work,Y);CHKERRABORT(PETSC_COMM_SELF,ierr); 53 ierr = VecDuplicate(cvode->work,&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr); 54 ierr = VecGetSize(yydot,&n);CHKERRABORT(PETSC_COMM_SELF,ierr); 55 ierr = MatCreateSeqDense(PETSC_COMM_SELF,n,n,DFY,&mat);CHKERRABORT(PETSC_COMM_SELF,ierr); 56 ierr = VecZeroEntries(yydot);CHKERRABORT(PETSC_COMM_SELF,ierr); 57 ierr = TSComputeIJacobian(ts,*X,cvode->work,yydot,0,mat,mat,PETSC_FALSE);CHKERRABORT(PETSC_COMM_SELF,ierr); 58 ierr = MatScale(mat,-1.0);CHKERRABORT(PETSC_COMM_SELF,ierr); 59 ierr = MatDestroy(&mat);CHKERRABORT(PETSC_COMM_SELF,ierr); 60 ierr = VecDestroy(&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr); 61 ierr = VecResetArray(cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr); 62 } 63 64 void SOLOUT(int *NR,double *XOLD,double *X, double *Y,double *CONT,double *LRC,int *N,double *RPAR,void *IPAR,int *IRTRN) 65 { 66 TS ts = (TS) IPAR; 67 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 68 PetscErrorCode ierr; 69 70 ierr = VecPlaceArray(cvode->work,Y);CHKERRABORT(PETSC_COMM_SELF,ierr); 71 ts->time_step = *X - *XOLD; 72 ierr = TSMonitor(ts,*NR-1,*X,cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr); 73 ierr = VecResetArray(cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr); 74 } 75 76 77 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*); 78 79 PetscErrorCode TSSolve_Radau5(TS ts) 80 { 81 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 82 PetscErrorCode ierr; 83 PetscScalar *Y,*WORK,X,XEND,RTOL,ATOL,H,RPAR; 84 PetscInt ND,*IWORK,LWORK,LIWORK,MUJAC,MLMAS,MUMAS,IDID,ITOL; 85 int IJAC,MLJAC,IMAS,IOUT; 86 87 PetscFunctionBegin; 88 ierr = VecGetArray(ts->vec_sol,&Y);CHKERRQ(ierr); 89 ierr = VecGetSize(ts->vec_sol,&ND);CHKERRQ(ierr); 90 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,ND,NULL,&cvode->work);CHKERRQ(ierr); 91 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,ND,NULL,&cvode->workf);CHKERRQ(ierr); 92 93 LWORK = 4*ND*ND+12*ND+20; 94 LIWORK = 3*ND+20; 95 96 ierr = PetscCalloc2(LWORK,&WORK,LIWORK,&IWORK);CHKERRQ(ierr); 97 98 /* C --- PARAMETER IN THE DIFFERENTIAL EQUATION */ 99 RPAR=1.0e-6; 100 /* C --- COMPUTE THE JACOBIAN ANALYTICALLY */ 101 IJAC=1; 102 /* C --- JACOBIAN IS A FULL MATRIX */ 103 MLJAC=ND; 104 /* C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM*/ 105 IMAS=0; 106 /* C --- OUTPUT ROUTINE IS USED DURING INTEGRATION*/ 107 IOUT=1; 108 /* C --- INITIAL VALUES*/ 109 X = ts->ptime; 110 /* C --- ENDPOINT OF INTEGRATION */ 111 XEND = ts->max_time; 112 /* C --- REQUIRED TOLERANCE */ 113 RTOL = ts->rtol; 114 ATOL = ts->atol; 115 ITOL=0; 116 /* C --- INITIAL STEP SIZE */ 117 H = ts->time_step; 118 119 /* output MUJAC MLMAS IDID; currently all ignored */ 120 121 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); 122 123 ierr = PetscFree2(WORK,IWORK);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 127 PetscErrorCode TSDestroy_Radau5(TS ts) 128 { 129 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 ierr = VecDestroy(&cvode->work);CHKERRQ(ierr); 134 ierr = VecDestroy(&cvode->workf);CHKERRQ(ierr); 135 ierr = PetscFree(ts->data);CHKERRQ(ierr); 136 PetscFunctionReturn(0); 137 } 138 139 /*MC 140 TSRADAU5 - ODE solver using the RADAU5 package 141 142 Notes: 143 This uses its own nonlinear solver and dense matrix direct solver so PETSc SNES and KSP options do not apply. 144 Uses its own time-step adaptivity (but uses the TS rtol and atol, and initial timestep) 145 Uses its own memory for the dense matrix storage and factorization 146 Can only handle ODEs of the form \cdot{u} = -F(t,u) + G(t,u) 147 148 Level: beginner 149 150 .seealso: TSCreate(), TS, TSSetType() 151 152 M*/ 153 PETSC_EXTERN PetscErrorCode TSCreate_Radau5(TS ts) 154 { 155 TS_Radau5 *cvode; 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 ts->ops->destroy = TSDestroy_Radau5; 160 ts->ops->solve = TSSolve_Radau5; 161 ts->default_adapt_type = TSADAPTNONE; 162 163 ierr = PetscNewLog(ts,&cvode);CHKERRQ(ierr); 164 ts->data = (void*)cvode; 165 PetscFunctionReturn(0); 166 } 167