Actual source code: ex1.c
1: static char help[] = "Basic problem for multi-rate method.\n";
3: /*F
5: \begin{eqnarray}
6: ys' = \frac{ys}{a}\\
7: yf' = ys*cos(bt)\\
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: PetscFunctionBegin;
24: PetscCall(VecGetArrayRead(U, &u));
25: PetscCall(VecGetArray(F, &f));
26: f[0] = u[0] / ctx->a;
27: f[1] = u[0] * PetscCosScalar(t * ctx->b);
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] = u[0] / ctx->a;
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] = u[0] * PetscCosScalar(t * ctx->b);
56: PetscCall(VecRestoreArrayRead(U, &u));
57: PetscCall(VecRestoreArray(F, &f));
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: /*
62: Define the analytic solution for check method easily
63: */
64: static PetscErrorCode sol_true(PetscReal t, Vec U, AppCtx *ctx)
65: {
66: PetscScalar *u;
68: PetscFunctionBegin;
69: PetscCall(VecGetArray(U, &u));
70: u[0] = PetscExpScalar(t / ctx->a);
71: u[1] = (ctx->a * PetscCosScalar(ctx->b * t) + ctx->a * ctx->a * ctx->b * PetscSinScalar(ctx->b * t)) * PetscExpScalar(t / ctx->a) / (1 + ctx->a * ctx->a * ctx->b * ctx->b);
72: PetscCall(VecRestoreArray(U, &u));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: int main(int argc, char **argv)
77: {
78: TS ts; /* ODE integrator */
79: Vec U; /* solution will be stored here */
80: Vec Utrue;
81: PetscMPIInt size;
82: AppCtx ctx;
83: PetscScalar *u;
84: IS iss;
85: IS isf;
86: PetscInt *indicess;
87: PetscInt *indicesf;
88: PetscInt n = 2;
89: PetscReal error, tt;
91: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: Initialize program
93: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: PetscFunctionBeginUser;
95: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
96: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
97: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Create index for slow part and fast part
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: PetscCall(PetscMalloc1(1, &indicess));
103: indicess[0] = 0;
104: PetscCall(PetscMalloc1(1, &indicesf));
105: indicesf[0] = 1;
106: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicess, PETSC_COPY_VALUES, &iss));
107: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicesf, PETSC_COPY_VALUES, &isf));
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Create necessary vector
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
113: PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
114: PetscCall(VecSetFromOptions(U));
115: PetscCall(VecDuplicate(U, &Utrue));
116: PetscCall(VecCopy(U, Utrue));
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Set runtime options
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", "");
122: {
123: ctx.a = 2.0;
124: ctx.b = 25.0;
125: PetscCall(PetscOptionsReal("-a", "", "", ctx.a, &ctx.a, NULL));
126: PetscCall(PetscOptionsReal("-b", "", "", ctx.b, &ctx.b, NULL));
127: ctx.Tf = 2;
128: ctx.dt = 0.01;
129: PetscCall(PetscOptionsReal("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL));
130: PetscCall(PetscOptionsReal("-dt", "", "", ctx.dt, &ctx.dt, NULL));
131: }
132: PetscOptionsEnd();
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Initialize the solution
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137: PetscCall(VecGetArray(U, &u));
138: u[0] = 1.0;
139: u[1] = ctx.a / (1 + ctx.a * ctx.a * ctx.b * ctx.b);
140: PetscCall(VecRestoreArray(U, &u));
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Create timestepping solver context
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
146: PetscCall(TSSetType(ts, TSMPRK));
148: PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunctionFn *)RHSFunction, &ctx));
149: PetscCall(TSRHSSplitSetIS(ts, "slow", iss));
150: PetscCall(TSRHSSplitSetIS(ts, "fast", isf));
151: PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunctionFn *)RHSFunctionslow, &ctx));
152: PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunctionFn *)RHSFunctionfast, &ctx));
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Set initial conditions
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: PetscCall(TSSetSolution(ts, U));
159: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: Set solver options
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162: PetscCall(TSSetMaxTime(ts, ctx.Tf));
163: PetscCall(TSSetTimeStep(ts, ctx.dt));
164: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
165: PetscCall(TSSetFromOptions(ts));
167: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: Solve linear system
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: PetscCall(TSSolve(ts, U));
171: PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
173: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: Check the error of the Petsc solution
175: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176: PetscCall(TSGetTime(ts, &tt));
177: PetscCall(sol_true(tt, Utrue, &ctx));
178: PetscCall(VecAXPY(Utrue, -1.0, U));
179: PetscCall(VecNorm(Utrue, NORM_2, &error));
181: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182: Print norm2 error
183: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error norm = %g\n", (double)error));
186: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: Free work space. All PETSc objects should be destroyed when they are no longer needed.
188: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: PetscCall(VecDestroy(&U));
190: PetscCall(VecDestroy(&Utrue));
191: PetscCall(TSDestroy(&ts));
192: PetscCall(ISDestroy(&iss));
193: PetscCall(ISDestroy(&isf));
194: PetscCall(PetscFree(indicess));
195: PetscCall(PetscFree(indicesf));
196: PetscCall(PetscFinalize());
197: return 0;
198: }
200: /*TEST
201: build:
202: requires: !complex
204: test:
206: TEST*/