1 static char help[] = "Basic problem for multi-rate method.\n";
2
3 /*F
4
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}
9
10 F*/
11
12 #include <petscts.h>
13
14 typedef struct {
15 PetscReal 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] = -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 }
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] = -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 }
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] = 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 }
60
sol_true(PetscReal t,Vec U)61 static PetscErrorCode sol_true(PetscReal t, Vec U)
62 {
63 PetscScalar *u;
64
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 }
72
main(int argc,char ** argv)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;
87
88 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89 Initialize program
90 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91 PetscFunctionBeginUser;
92 PetscCall(PetscInitialize(&argc, &argv, NULL, 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");
95
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));
105
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));
114
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));
122
123 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124 Create timestepping solver context
125 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
127 PetscCall(TSSetType(ts, TSMPRK));
128
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));
134
135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136 Set initial conditions
137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138 PetscCall(TSSetSolution(ts, U));
139
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));
155
156 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157 Solve linear system
158 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159 PetscCall(TSSolve(ts, U));
160 PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
161
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));
169
170 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171 Print norm2 error
172 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error norm: %g\n", (double)error));
174
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 }
188
189 /*TEST
190 build:
191 requires: !complex
192
193 test:
194
195 TEST*/
196