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: PetscFunctionBegin;
24: PetscCall(VecGetArrayRead(U, &u));
25: PetscCall(VecGetArray(F, &f));
26: 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]);
27: 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]);
28: PetscCall(VecRestoreArrayRead(U, &u));
29: PetscCall(VecRestoreArray(F, &f));
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: static PetscErrorCode RHSFunctionslow(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
34: {
35: const PetscScalar *u;
36: PetscScalar *f;
38: PetscFunctionBegin;
39: PetscCall(VecGetArrayRead(U, &u));
40: PetscCall(VecGetArray(F, &f));
41: 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]);
42: PetscCall(VecRestoreArrayRead(U, &u));
43: PetscCall(VecRestoreArray(F, &f));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: static PetscErrorCode RHSFunctionfast(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
48: {
49: const PetscScalar *u;
50: PetscScalar *f;
52: PetscFunctionBegin;
53: PetscCall(VecGetArrayRead(U, &u));
54: PetscCall(VecGetArray(F, &f));
55: 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]);
56: PetscCall(VecRestoreArrayRead(U, &u));
57: PetscCall(VecRestoreArray(F, &f));
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode sol_true(PetscReal t, Vec U)
62: {
63: PetscScalar *u;
65: PetscFunctionBegin;
66: PetscCall(VecGetArray(U, &u));
67: u[0] = PetscSqrtScalar(1.0 + PetscCosScalar(t));
68: u[1] = PetscSqrtScalar(2.0 + PetscCosScalar(5.0 * t));
69: PetscCall(VecRestoreArray(U, &u));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: int main(int argc, char **argv)
74: {
75: TS ts; /* ODE integrator */
76: Vec U; /* solution will be stored here */
77: 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: PetscReal error, tt;
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Initialize program
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: PetscFunctionBeginUser;
92: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
93: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
94: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Create index for slow part and fast part
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: PetscCall(PetscMalloc1(1, &indicess));
100: indicess[0] = 0;
101: PetscCall(PetscMalloc1(1, &indicesf));
102: indicesf[0] = 1;
103: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicess, PETSC_COPY_VALUES, &iss));
104: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicesf, PETSC_COPY_VALUES, &isf));
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Create necessary vector
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
110: PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
111: PetscCall(VecSetFromOptions(U));
112: PetscCall(VecDuplicate(U, &Utrue));
113: PetscCall(VecCopy(U, Utrue));
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Set initial condition
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118: PetscCall(VecGetArray(U, &u));
119: u[0] = PetscSqrtScalar(2.0);
120: u[1] = PetscSqrtScalar(3.0);
121: PetscCall(VecRestoreArray(U, &u));
123: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: Create timestepping solver context
125: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
127: PetscCall(TSSetType(ts, TSMPRK));
129: PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunctionFn *)RHSFunction, &ctx));
130: PetscCall(TSRHSSplitSetIS(ts, "slow", iss));
131: PetscCall(TSRHSSplitSetIS(ts, "fast", isf));
132: PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunctionFn *)RHSFunctionslow, &ctx));
133: PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunctionFn *)RHSFunctionfast, &ctx));
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Set initial conditions
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: PetscCall(TSSetSolution(ts, U));
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Set solver options
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", "");
144: {
145: ctx.Tf = 0.3;
146: ctx.dt = 0.01;
147: PetscCall(PetscOptionsScalar("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL));
148: PetscCall(PetscOptionsScalar("-dt", "", "", ctx.dt, &ctx.dt, NULL));
149: }
150: PetscOptionsEnd();
151: PetscCall(TSSetMaxTime(ts, ctx.Tf));
152: PetscCall(TSSetTimeStep(ts, ctx.dt));
153: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
154: PetscCall(TSSetFromOptions(ts));
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Solve linear system
158: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159: PetscCall(TSSolve(ts, U));
160: PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Check the error of the Petsc solution
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: PetscCall(TSGetTime(ts, &tt));
166: PetscCall(sol_true(tt, Utrue));
167: PetscCall(VecAXPY(Utrue, -1.0, U));
168: PetscCall(VecNorm(Utrue, NORM_2, &error));
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Print norm2 error
172: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error norm: %g\n", (double)error));
175: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: Free work space. All PETSc objects should be destroyed when they are no longer needed.
177: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
178: PetscCall(VecDestroy(&U));
179: PetscCall(TSDestroy(&ts));
180: PetscCall(VecDestroy(&Utrue));
181: PetscCall(ISDestroy(&iss));
182: PetscCall(ISDestroy(&isf));
183: PetscCall(PetscFree(indicess));
184: PetscCall(PetscFree(indicesf));
185: PetscCall(PetscFinalize());
186: return 0;
187: }
189: /*TEST
190: build:
191: requires: !complex
193: test:
195: TEST*/