Actual source code: ex3.c
1: static char help[] = "Basic problem for multi-rate method.\n";
3: /*F
5: \begin{eqnarray}
6: 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}\\
7: 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}\\
8: \end{eqnarray}
10: F*/
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: const PetscScalar *u;
21: PetscScalar *f;
23: VecGetArrayRead(U, &u);
24: VecGetArray(F, &f);
25: 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]);
26: 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]);
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] = -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]);
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] = 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]);
53: VecRestoreArrayRead(U, &u);
54: VecRestoreArray(F, &f);
55: return 0;
56: }
58: static PetscErrorCode sol_true(PetscReal t, Vec U)
59: {
60: PetscScalar *u;
62: VecGetArray(U, &u);
63: u[0] = PetscSqrtScalar(1.0 + PetscCosScalar(t));
64: u[1] = PetscSqrtScalar(2.0 + PetscCosScalar(5.0 * t));
65: VecRestoreArray(U, &u);
66: return 0;
67: }
69: int main(int argc, char **argv)
70: {
71: TS ts; /* ODE integrator */
72: Vec U; /* solution will be stored here */
73: Vec Utrue;
74: PetscMPIInt size;
75: AppCtx ctx;
76: PetscScalar *u;
77: IS iss;
78: IS isf;
79: PetscInt *indicess;
80: PetscInt *indicesf;
81: PetscInt n = 2;
82: PetscReal error, tt;
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Initialize program
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: PetscInitialize(&argc, &argv, (char *)0, help);
89: MPI_Comm_size(PETSC_COMM_WORLD, &size);
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Create index for slow part and fast part
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: PetscMalloc1(1, &indicess);
96: indicess[0] = 0;
97: PetscMalloc1(1, &indicesf);
98: indicesf[0] = 1;
99: ISCreateGeneral(PETSC_COMM_SELF, 1, indicess, PETSC_COPY_VALUES, &iss);
100: ISCreateGeneral(PETSC_COMM_SELF, 1, indicesf, PETSC_COPY_VALUES, &isf);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Create necessary vector
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: VecCreate(PETSC_COMM_WORLD, &U);
106: VecSetSizes(U, n, PETSC_DETERMINE);
107: VecSetFromOptions(U);
108: VecDuplicate(U, &Utrue);
109: VecCopy(U, Utrue);
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Set initial condition
113: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: VecGetArray(U, &u);
115: u[0] = PetscSqrtScalar(2.0);
116: u[1] = PetscSqrtScalar(3.0);
117: VecRestoreArray(U, &u);
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Create timestepping solver context
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122: TSCreate(PETSC_COMM_WORLD, &ts);
123: TSSetType(ts, TSMPRK);
125: TSSetRHSFunction(ts, NULL, (TSRHSFunction)RHSFunction, &ctx);
126: TSRHSSplitSetIS(ts, "slow", iss);
127: TSRHSSplitSetIS(ts, "fast", isf);
128: TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunction)RHSFunctionslow, &ctx);
129: TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunction)RHSFunctionfast, &ctx);
131: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: Set initial conditions
133: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: TSSetSolution(ts, U);
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Set solver options
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", "");
140: {
141: ctx.Tf = 0.3;
142: ctx.dt = 0.01;
143: PetscOptionsScalar("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL);
144: PetscOptionsScalar("-dt", "", "", ctx.dt, &ctx.dt, NULL);
145: }
146: PetscOptionsEnd();
147: TSSetMaxTime(ts, ctx.Tf);
148: TSSetTimeStep(ts, ctx.dt);
149: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
150: TSSetFromOptions(ts);
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Solve linear system
154: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155: TSSolve(ts, U);
156: VecView(U, PETSC_VIEWER_STDOUT_WORLD);
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Check the error of the Petsc solution
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: TSGetTime(ts, &tt);
162: sol_true(tt, Utrue);
163: VecAXPY(Utrue, -1.0, U);
164: VecNorm(Utrue, NORM_2, &error);
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Print norm2 error
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: PetscPrintf(PETSC_COMM_WORLD, "l2 error norm: %g\n", (double)error);
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Free work space. All PETSc objects should be destroyed when they are no longer needed.
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174: VecDestroy(&U);
175: TSDestroy(&ts);
176: VecDestroy(&Utrue);
177: ISDestroy(&iss);
178: ISDestroy(&isf);
179: PetscFree(indicess);
180: PetscFree(indicesf);
181: PetscFinalize();
182: return 0;
183: }
185: /*TEST
186: build:
187: requires: !complex
189: test:
191: TEST*/