Actual source code: ex2.c
1: static char help[] = "Basic problem for multi-rate method.\n";
3: /*F
5: \begin{eqnarray}
6: ys' = -sin(a*t)\\
7: yf' = bcos(b*t)ys-sin(b*t)sin(a*t)\\
8: \end{eqnarray}
10: F*/
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;
77: PetscMPIInt size;
78: AppCtx ctx;
79: PetscScalar *u;
80: IS iss;
81: IS isf;
82: PetscInt *indicess;
83: PetscInt *indicesf;
84: PetscInt n = 2;
85: PetscScalar error, tt;
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Initialize program
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
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", (double)PetscAbsScalar(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*/