Actual source code: radau5.c
petsc-3.10.5 2019-03-28
1: /*
2: Provides a PETSc interface to RADAU5 solver.
4: */
5: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
7: typedef struct {
8: Vec work,workf;
9: } TS_Radau5;
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;
20: VecPlaceArray(cvode->work,Y);CHKERRABORT(PETSC_COMM_SELF,ierr);
21: VecPlaceArray(cvode->workf,F);CHKERRABORT(PETSC_COMM_SELF,ierr);
23: /* Now compute the right hand side function, via IFunction unless only the more efficient RHSFunction is set */
24: TSGetDM(ts,&dm);CHKERRABORT(PETSC_COMM_SELF,ierr);
25: DMGetDMTS(dm,&tsdm);CHKERRABORT(PETSC_COMM_SELF,ierr);
26: DMTSGetIFunction(dm,&ifunction,NULL);CHKERRABORT(PETSC_COMM_SELF,ierr);
27: if (!ifunction) {
28: 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;
32: VecDuplicate(cvode->work,&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
33: VecZeroEntries(yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
34: TSComputeIFunction(ts,*X,cvode->work,yydot,cvode->workf,PETSC_FALSE);CHKERRABORT(PETSC_COMM_SELF,ierr);
35: VecScale(cvode->workf,-1.);CHKERRABORT(PETSC_COMM_SELF,ierr);
36: VecDestroy(&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
37: }
39: VecResetArray(cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr);
40: VecResetArray(cvode->workf);CHKERRABORT(PETSC_COMM_SELF,ierr);
41: }
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;
52: VecPlaceArray(cvode->work,Y);CHKERRABORT(PETSC_COMM_SELF,ierr);
53: VecDuplicate(cvode->work,&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
54: VecGetSize(yydot,&n);CHKERRABORT(PETSC_COMM_SELF,ierr);
55: MatCreateSeqDense(PETSC_COMM_SELF,n,n,DFY,&mat);CHKERRABORT(PETSC_COMM_SELF,ierr);
56: VecZeroEntries(yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
57: TSComputeIJacobian(ts,*X,cvode->work,yydot,0,mat,mat,PETSC_FALSE);CHKERRABORT(PETSC_COMM_SELF,ierr);
58: MatScale(mat,-1.0);CHKERRABORT(PETSC_COMM_SELF,ierr);
59: MatDestroy(&mat);CHKERRABORT(PETSC_COMM_SELF,ierr);
60: VecDestroy(&yydot);CHKERRABORT(PETSC_COMM_SELF,ierr);
61: VecResetArray(cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr);
62: }
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;
70: VecPlaceArray(cvode->work,Y);CHKERRABORT(PETSC_COMM_SELF,ierr);
71: ts->time_step = *X - *XOLD;
72: TSMonitor(ts,*NR-1,*X,cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr);
73: VecResetArray(cvode->work);CHKERRABORT(PETSC_COMM_SELF,ierr);
74: }
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*);
79: PetscErrorCode TSSolve_Radau5(TS ts)
80: {
81: TS_Radau5 *cvode = (TS_Radau5*)ts->data;
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;
88: VecGetArray(ts->vec_sol,&Y);
89: VecGetSize(ts->vec_sol,&ND);
90: VecCreateSeqWithArray(PETSC_COMM_SELF,1,ND,NULL,&cvode->work);
91: VecCreateSeqWithArray(PETSC_COMM_SELF,1,ND,NULL,&cvode->workf);
93: LWORK = 4*ND*ND+12*ND+20;
94: LIWORK = 3*ND+20;
96: PetscCalloc1(LWORK,&WORK);
97: PetscCalloc1(LIWORK,&IWORK);
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;
120: /* output MUJAC MLMAS IDID; currently all ignored */
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);
124: PetscFree(WORK);
125: PetscFree(IWORK);
126: return(0);
127: }
129: PetscErrorCode TSDestroy_Radau5(TS ts)
130: {
131: TS_Radau5 *cvode = (TS_Radau5*)ts->data;
135: VecDestroy(&cvode->work);
136: VecDestroy(&cvode->workf);
137: PetscFree(ts->data);
138: return(0);
139: }
141: /*MC
142: TSRADAU5 - ODE solver using the RADAU5 package
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)
150: Level: beginner
152: .seealso: TSCreate(), TS, TSSetType()
154: M*/
155: PETSC_EXTERN PetscErrorCode TSCreate_Radau5(TS ts)
156: {
157: TS_Radau5 *cvode;
161: ts->ops->destroy = TSDestroy_Radau5;
162: ts->ops->solve = TSSolve_Radau5;
163: ts->default_adapt_type = TSADAPTNONE;
165: PetscNewLog(ts,&cvode);
166: ts->data = (void*)cvode;
167: return(0);
168: }