1 static char help[] = "Basic problem for multi-rate method.\n";
2
3 /*F
4
5 \begin{eqnarray}
6 ys' = \frac{ys}{a}\\
7 yf' = ys*cos(bt)\\
8 \end{eqnarray}
9
10 F*/
11
12 #include <petscts.h>
13
14 typedef struct {
15 PetscReal a, b, Tf, dt;
16 } AppCtx;
17
RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx * ctx)18 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
19 {
20 const PetscScalar *u;
21 PetscScalar *f;
22
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 }
32
RHSFunctionslow(TS ts,PetscReal t,Vec U,Vec F,AppCtx * ctx)33 static PetscErrorCode RHSFunctionslow(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
34 {
35 const PetscScalar *u;
36 PetscScalar *f;
37
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 }
46
RHSFunctionfast(TS ts,PetscReal t,Vec U,Vec F,AppCtx * ctx)47 static PetscErrorCode RHSFunctionfast(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
48 {
49 const PetscScalar *u;
50 PetscScalar *f;
51
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 }
60
61 /*
62 Define the analytic solution for check method easily
63 */
sol_true(PetscReal t,Vec U,AppCtx * ctx)64 static PetscErrorCode sol_true(PetscReal t, Vec U, AppCtx *ctx)
65 {
66 PetscScalar *u;
67
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 }
75
main(int argc,char ** argv)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;
90
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");
98
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));
108
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));
117
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();
133
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));
141
142 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143 Create timestepping solver context
144 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
146 PetscCall(TSSetType(ts, TSMPRK));
147
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));
153
154 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155 Set initial conditions
156 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157 PetscCall(TSSetSolution(ts, U));
158
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));
166
167 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168 Solve linear system
169 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170 PetscCall(TSSolve(ts, U));
171 PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
172
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));
180
181 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182 Print norm2 error
183 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error norm = %g\n", (double)error));
185
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 }
199
200 /*TEST
201 build:
202 requires: !complex
203
204 test:
205
206 TEST*/
207