Actual source code: ex3.c
petsc-3.12.0 2019-09-29
1: static char help[] = "Basic problem for multi-rate method.\n";
\begin{eqnarray}
ys' = -2.0*\frac{-1.0+ys^2.0-\cos(t)}{2.0*ys}+0.05*\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-\frac{\sin(t)}{2.0*ys}\\
yf' = 0.05*\frac{-1.0+ys^2-\cos(t)}{2.0*ys}-\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-5.0*\frac{\sin(5.0*t)}{2.0*yf}\\
\end{eqnarray}
12: #include <petscts.h>
14: typedef struct {
15: PetscReal Tf,dt;
16: } AppCtx;
18: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
19: {
20: PetscErrorCode ierr;
21: const PetscScalar *u;
22: PetscScalar *f;
25: VecGetArrayRead(U,&u);
26: VecGetArray(F,&f);
27: f[0] = -2.0*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])+0.05*(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-PetscSinScalar(t)/(2.0*u[0]);
28: f[1] = 0.05*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])-(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-5.0*PetscSinScalar(5.0*t)/(2.0*u[1]);
29: VecRestoreArrayRead(U,&u);
30: VecRestoreArray(F,&f);
31: return(0);
32: }
34: static PetscErrorCode RHSFunctionslow(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
35: {
36: PetscErrorCode ierr;
37: const PetscScalar *u;
38: PetscScalar *f;
41: VecGetArrayRead(U,&u);
42: VecGetArray(F,&f);
43: f[0] = -2.0*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])+0.05*(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-PetscSinScalar(t)/(2.0*u[0]);
44: VecRestoreArrayRead(U,&u);
45: VecRestoreArray(F,&f);
46: return(0);
47: }
49: static PetscErrorCode RHSFunctionfast(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
50: {
51: PetscErrorCode ierr;
52: const PetscScalar *u;
53: PetscScalar *f;
56: VecGetArrayRead(U,&u);
57: VecGetArray(F,&f);
58: f[0] = 0.05*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])-(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-5.0*PetscSinScalar(5.0*t)/(2.0*u[1]);
59: VecRestoreArrayRead(U,&u);
60: VecRestoreArray(F,&f);
61: return(0);
62: }
64: static PetscErrorCode sol_true(PetscReal t,Vec U)
65: {
67: PetscScalar *u;
70: VecGetArray(U,&u);
71: u[0] = PetscSqrtScalar(1.0+PetscCosScalar(t));
72: u[1] = PetscSqrtScalar(2.0+PetscCosScalar(5.0*t));
73: VecRestoreArray(U,&u);
74: return(0);
75: }
77: int main(int argc,char **argv)
78: {
79: TS ts; /* ODE integrator */
80: Vec U; /* solution will be stored here */
81: Vec Utrue;
83: PetscMPIInt size;
84: AppCtx ctx;
85: PetscScalar *u;
86: IS iss;
87: IS isf;
88: PetscInt *indicess;
89: PetscInt *indicesf;
90: PetscInt n=2;
91: PetscReal error,tt;
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Initialize program
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
97: MPI_Comm_size(PETSC_COMM_WORLD,&size);
98: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create index for slow part and fast part
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: PetscMalloc1(1,&indicess);
104: indicess[0]=0;
105: PetscMalloc1(1,&indicesf);
106: indicesf[0]=1;
107: ISCreateGeneral(PETSC_COMM_SELF,1,indicess,PETSC_COPY_VALUES,&iss);
108: ISCreateGeneral(PETSC_COMM_SELF,1,indicesf,PETSC_COPY_VALUES,&isf);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Create necesary vector
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: VecCreate(PETSC_COMM_WORLD,&U);
114: VecSetSizes(U,n,PETSC_DETERMINE);
115: VecSetFromOptions(U);
116: VecDuplicate(U,&Utrue);
117: VecCopy(U,Utrue);
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Set initial condition
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: VecGetArray(U,&u);
123: u[0] = PetscSqrtScalar(2.0);
124: u[1] = PetscSqrtScalar(3.0);
125: VecRestoreArray(U,&u);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Create timestepping solver context
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: TSCreate(PETSC_COMM_WORLD,&ts);
131: TSSetType(ts,TSMPRK);
133: TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
134: TSRHSSplitSetIS(ts,"slow",iss);
135: TSRHSSplitSetIS(ts,"fast",isf);
136: TSRHSSplitSetRHSFunction(ts,"slow",NULL,(TSRHSFunctionslow)RHSFunctionslow,&ctx);
137: TSRHSSplitSetRHSFunction(ts,"fast",NULL,(TSRHSFunctionfast)RHSFunctionfast,&ctx);
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Set initial conditions
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142: TSSetSolution(ts,U);
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Set solver options
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ODE options","");
148: {
149: ctx.Tf = 0.3;
150: ctx.dt = 0.01;
151: PetscOptionsScalar("-Tf","","",ctx.Tf,&ctx.Tf,NULL);
152: PetscOptionsScalar("-dt","","",ctx.dt,&ctx.dt,NULL);
153: }
154: PetscOptionsEnd();
155: TSSetMaxTime(ts,ctx.Tf);
156: TSSetTimeStep(ts,ctx.dt);
157: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
158: TSSetFromOptions(ts);
160: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: Solve linear system
162: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163: TSSolve(ts,U);
164: VecView(U,PETSC_VIEWER_STDOUT_WORLD);
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Check the error of the Petsc solution
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: TSGetTime(ts,&tt);
170: sol_true(tt,Utrue);
171: VecAXPY(Utrue,-1.0,U);
172: VecNorm(Utrue,NORM_2,&error);
174: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: Print norm2 error
176: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177: PetscPrintf(PETSC_COMM_WORLD,"l2 error norm: %g\n", error);
179: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180: Free work space. All PETSc objects should be destroyed when they are no longer needed.
181: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182: VecDestroy(&U);
183: TSDestroy(&ts);
184: VecDestroy(&Utrue);
185: ISDestroy(&iss);
186: ISDestroy(&isf);
187: PetscFree(indicess);
188: PetscFree(indicesf);
189: PetscFinalize();
190: return ierr;
191: }