Actual source code: ex2.c
1: static char help[] = "Basic problem for multi-rate method.\n";
\begin{eqnarray}
ys' = -sin(a*t)\\
yf' = bcos(b*t)ys-sin(b*t)sin(a*t)\\
\end{eqnarray}
12: #include <petscts.h>
14: typedef struct {
15: PetscReal a,b,Tf,dt;
16: } AppCtx;
18: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
19: {
20: const PetscScalar *u;
21: PetscScalar *f;
23: VecGetArrayRead(U,&u);
24: VecGetArray(F,&f);
25: f[0] = -PetscSinScalar(ctx->a*t);
26: f[1] = ctx->b*PetscCosScalar(ctx->b*t)*u[0]-PetscSinScalar(ctx->b*t)*PetscSinScalar(ctx->a*t);
27: VecRestoreArrayRead(U,&u);
28: VecRestoreArray(F,&f);
29: return 0;
30: }
32: static PetscErrorCode RHSFunctionslow(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
33: {
34: const PetscScalar *u;
35: PetscScalar *f;
37: VecGetArrayRead(U,&u);
38: VecGetArray(F,&f);
39: f[0] = -PetscSinScalar(ctx->a*t);
40: VecRestoreArrayRead(U,&u);
41: VecRestoreArray(F,&f);
42: return 0;
43: }
45: static PetscErrorCode RHSFunctionfast(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
46: {
47: const PetscScalar *u;
48: PetscScalar *f;
50: VecGetArrayRead(U,&u);
51: VecGetArray(F,&f);
52: f[0] = ctx->b*PetscCosScalar(ctx->b*t)*u[0]-PetscSinScalar(ctx->b*t)*PetscSinScalar(ctx->a*t);
53: VecRestoreArrayRead(U,&u);
54: VecRestoreArray(F,&f);
55: return 0;
56: }
58: /*
59: Define the analytic solution for check method easily
60: */
61: static PetscErrorCode sol_true(PetscReal t,Vec U,AppCtx *ctx)
62: {
63: PetscScalar *u;
65: VecGetArray(U,&u);
66: u[0] = PetscCosScalar(ctx->a*t)/ctx->a;
67: u[1] = PetscSinScalar(ctx->b*t)*PetscCosScalar(ctx->a*t)/ctx->a;
68: VecRestoreArray(U,&u);
69: return 0;
70: }
72: int main(int argc,char **argv)
73: {
74: TS ts; /* ODE integrator */
75: Vec U; /* solution will be stored here */
76: Vec Utrue;
78: PetscMPIInt size;
79: AppCtx ctx;
80: PetscScalar *u;
81: IS iss;
82: IS isf;
83: PetscInt *indicess;
84: PetscInt *indicesf;
85: PetscInt n=2;
86: PetscScalar error,tt;
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Initialize program
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: PetscInitialize(&argc,&argv,(char*)0,help);
92: MPI_Comm_size(PETSC_COMM_WORLD,&size);
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Create index for slow part and fast part
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: PetscMalloc1(1,&indicess);
99: indicess[0]=0;
100: PetscMalloc1(1,&indicesf);
101: indicesf[0]=1;
102: ISCreateGeneral(PETSC_COMM_SELF,1,indicess,PETSC_COPY_VALUES,&iss);
103: ISCreateGeneral(PETSC_COMM_SELF,1,indicesf,PETSC_COPY_VALUES,&isf);
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Create necessary vector
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: VecCreate(PETSC_COMM_WORLD,&U);
109: VecSetSizes(U,n,PETSC_DETERMINE);
110: VecSetFromOptions(U);
111: VecDuplicate(U,&Utrue);
112: VecCopy(U,Utrue);
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Set runtime options
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ODE options","");
118: {
119: ctx.a = 1.0;
120: ctx.b = 25.0;
121: PetscOptionsScalar("-a","","",ctx.a,&ctx.a,NULL);
122: PetscOptionsScalar("-b","","",ctx.b,&ctx.b,NULL);
123: ctx.Tf = 5.0;
124: ctx.dt = 0.01;
125: PetscOptionsScalar("-Tf","","",ctx.Tf,&ctx.Tf,NULL);
126: PetscOptionsScalar("-dt","","",ctx.dt,&ctx.dt,NULL);
127: }
128: PetscOptionsEnd();
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Initialize the solution
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: VecGetArray(U,&u);
134: u[0] = 1.0/ctx.a;
135: u[1] = 0.0;
136: VecRestoreArray(U,&u);
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Create timestepping solver context
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: TSCreate(PETSC_COMM_WORLD,&ts);
142: TSSetType(ts,TSMPRK);
144: TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
145: TSRHSSplitSetIS(ts,"slow",iss);
146: TSRHSSplitSetIS(ts,"fast",isf);
147: TSRHSSplitSetRHSFunction(ts,"slow",NULL,(TSRHSFunction)RHSFunctionslow,&ctx);
148: TSRHSSplitSetRHSFunction(ts,"fast",NULL,(TSRHSFunction)RHSFunctionfast,&ctx);
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Set initial conditions
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: TSSetSolution(ts,U);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: Set solver options
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: TSSetMaxTime(ts,ctx.Tf);
159: TSSetTimeStep(ts,ctx.dt);
160: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
161: TSSetFromOptions(ts);
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Solve linear system
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: TSSolve(ts,U);
167: VecView(U,PETSC_VIEWER_STDOUT_WORLD);
169: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170: Check the error of the Petsc solution
171: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
172: TSGetTime(ts,&tt);
173: sol_true(tt,Utrue,&ctx);
174: VecAXPY(Utrue,-1.0,U);
175: VecNorm(Utrue,NORM_2,&error);
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Print norm2 error
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: PetscPrintf(PETSC_COMM_WORLD,"l2 error norm: %g\n", error);
182: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: Free work space. All PETSc objects should be destroyed when they are no longer needed.
184: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185: VecDestroy(&U);
186: TSDestroy(&ts);
187: VecDestroy(&Utrue);
188: ISDestroy(&iss);
189: ISDestroy(&isf);
190: PetscFree(indicess);
191: PetscFree(indicesf);
192: PetscFinalize();
193: return 0;
194: }
196: /*TEST
197: build:
198: requires: !complex
200: test:
202: TEST*/