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 = PetscCalloc1(LWORK,&WORK);CHKERRQ(ierr); 97 ierr = PetscCalloc1(LIWORK,&IWORK);CHKERRQ(ierr); 98 99 /* C --- PARAMETER IN THE DIFFERENTIAL EQUATION */ 100 RPAR=1.0e-6; 101 /* C --- COMPUTE THE JACOBIAN ANALYTICALLY */ 102 IJAC=1; 103 /* C --- JACOBIAN IS A FULL MATRIX */ 104 MLJAC=ND; 105 /* C --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM*/ 106 IMAS=0; 107 /* C --- OUTPUT ROUTINE IS USED DURING INTEGRATION*/ 108 IOUT=1; 109 /* C --- INITIAL VALUES*/ 110 X = ts->ptime; 111 /* C --- ENDPOINT OF INTEGRATION */ 112 XEND = ts->max_time; 113 /* C --- REQUIRED TOLERANCE */ 114 RTOL = ts->rtol; 115 ATOL = ts->atol; 116 ITOL=0; 117 /* C --- INITIAL STEP SIZE */ 118 H = ts->time_step; 119 120 /* output MUJAC MLMAS IDID; currently all ignored */ 121 122 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); 123 124 ierr = PetscFree(WORK);CHKERRQ(ierr); 125 ierr = PetscFree(IWORK);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 PetscErrorCode TSDestroy_Radau5(TS ts) 130 { 131 TS_Radau5 *cvode = (TS_Radau5*)ts->data; 132 PetscErrorCode ierr; 133 134 PetscFunctionBegin; 135 ierr = VecDestroy(&cvode->work);CHKERRQ(ierr); 136 ierr = VecDestroy(&cvode->workf);CHKERRQ(ierr); 137 ierr = PetscFree(ts->data);CHKERRQ(ierr); 138 PetscFunctionReturn(0); 139 } 140 141 /*MC 142 TSRADAU5 - ODE solver using the RADAU5 package 143 144 Notes: 145 This uses its own nonlinear solver and dense matrix direct solver so PETSc SNES and KSP options do not apply. 146 Uses its own time-step adaptivity (but uses the TS rtol and atol, and initial timestep) 147 Uses its own memory for the dense matrix storage and factorization 148 Can only handle ODEs of the form \cdot{u} = -F(t,u) + G(t,u) 149 150 Level: beginner 151 152 .seealso: TSCreate(), TS, TSSetType() 153 154 M*/ 155 PETSC_EXTERN PetscErrorCode TSCreate_Radau5(TS ts) 156 { 157 TS_Radau5 *cvode; 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 ts->ops->destroy = TSDestroy_Radau5; 162 ts->ops->solve = TSSolve_Radau5; 163 ts->default_adapt_type = TSADAPTNONE; 164 165 ierr = PetscNewLog(ts,&cvode);CHKERRQ(ierr); 166 ts->data = (void*)cvode; 167 PetscFunctionReturn(0); 168 } 169