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