1 static char help[] = "A fast-slow system for testing ARKIMEX.\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 /*
13 This example demonstrates how to use ARKIMEX for solving a fast-slow system. The system is partitioned additively and component-wise at the same time.
14 ys stands for the slow component and yf stands for the fast component. On the RHS for yf, only the term -\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf} is treated implicitly while the rest is treated explicitly.
15 */
16
17 #include <petscts.h>
18
19 typedef struct {
20 PetscReal Tf, dt;
21 } AppCtx;
22
RHSFunctionslow(TS ts,PetscReal t,Vec U,Vec F,AppCtx * ctx)23 static PetscErrorCode RHSFunctionslow(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
24 {
25 const PetscScalar *u;
26 PetscScalar *f;
27
28 PetscFunctionBegin;
29 PetscCall(VecGetArrayRead(U, &u));
30 PetscCall(VecGetArray(F, &f));
31 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]);
32 PetscCall(VecRestoreArrayRead(U, &u));
33 PetscCall(VecRestoreArray(F, &f));
34 PetscFunctionReturn(PETSC_SUCCESS);
35 }
36
RHSFunctionfast(TS ts,PetscReal t,Vec U,Vec F,AppCtx * ctx)37 static PetscErrorCode RHSFunctionfast(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
38 {
39 const PetscScalar *u;
40 PetscScalar *f;
41
42 PetscFunctionBegin;
43 PetscCall(VecGetArrayRead(U, &u));
44 PetscCall(VecGetArray(F, &f));
45 f[0] = 0.05 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) - 5.0 * PetscSinScalar(5.0 * t) / (2.0 * u[1]);
46 PetscCall(VecRestoreArrayRead(U, &u));
47 PetscCall(VecRestoreArray(F, &f));
48 PetscFunctionReturn(PETSC_SUCCESS);
49 }
50
IFunctionfast(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx * ctx)51 static PetscErrorCode IFunctionfast(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
52 {
53 const PetscScalar *u, *udot;
54 PetscScalar *f;
55
56 PetscFunctionBegin;
57 PetscCall(VecGetArrayRead(U, &u));
58 PetscCall(VecGetArrayRead(Udot, &udot));
59 PetscCall(VecGetArray(F, &f));
60 f[0] = udot[0] + (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]);
61 PetscCall(VecRestoreArrayRead(Udot, &udot));
62 PetscCall(VecRestoreArrayRead(U, &u));
63 PetscCall(VecRestoreArray(F, &f));
64 PetscFunctionReturn(PETSC_SUCCESS);
65 }
66
IJacobianfast(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx * ctx)67 static PetscErrorCode IJacobianfast(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
68 {
69 PetscInt rowcol[] = {0};
70 const PetscScalar *u;
71 PetscScalar J[1][1];
72
73 PetscFunctionBeginUser;
74 PetscCall(VecGetArrayRead(U, &u));
75 J[0][0] = a + (2.0 + PetscCosScalar(5.0 * t)) / (2.0 * u[1] * u[1]) + 0.5;
76 PetscCall(MatSetValues(B, 1, rowcol, 1, rowcol, &J[0][0], INSERT_VALUES));
77 PetscCall(VecRestoreArrayRead(U, &u));
78
79 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
80 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
81 if (A != B) {
82 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
83 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
84 }
85 PetscFunctionReturn(PETSC_SUCCESS);
86 }
87
88 /*
89 Define the analytic solution for check method easily
90 */
sol_true(PetscReal t,Vec U)91 static PetscErrorCode sol_true(PetscReal t, Vec U)
92 {
93 PetscScalar *u;
94
95 PetscFunctionBegin;
96 PetscCall(VecGetArray(U, &u));
97 u[0] = PetscSqrtScalar(1.0 + PetscCosScalar(t));
98 u[1] = PetscSqrtScalar(2.0 + PetscCosScalar(5.0 * t));
99 PetscCall(VecRestoreArray(U, &u));
100 PetscFunctionReturn(PETSC_SUCCESS);
101 }
102
main(int argc,char ** argv)103 int main(int argc, char **argv)
104 {
105 TS ts; /* ODE integrator */
106 Vec U; /* solution will be stored here */
107 Vec Utrue;
108 Mat A;
109 PetscMPIInt size;
110 AppCtx ctx;
111 PetscScalar *u;
112 IS iss;
113 IS isf;
114 PetscInt *indicess;
115 PetscInt *indicesf;
116 PetscInt n = 2;
117 PetscScalar error, tt;
118
119 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120 Initialize program
121 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
122 PetscFunctionBeginUser;
123 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
124 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
125 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
126
127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128 Create index for slow part and fast part
129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130 PetscCall(PetscMalloc1(1, &indicess));
131 indicess[0] = 0;
132 PetscCall(PetscMalloc1(1, &indicesf));
133 indicesf[0] = 1;
134 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicess, PETSC_COPY_VALUES, &iss));
135 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicesf, PETSC_COPY_VALUES, &isf));
136
137 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138 Create necessary vector
139 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140 PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
141 PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
142 PetscCall(VecSetFromOptions(U));
143 PetscCall(VecDuplicate(U, &Utrue));
144 PetscCall(VecCopy(U, Utrue));
145
146 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147 Set runtime options
148 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", "");
150 ctx.Tf = 0.3;
151 ctx.dt = 0.01;
152 PetscCall(PetscOptionsScalar("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL));
153 PetscCall(PetscOptionsScalar("-dt", "", "", ctx.dt, &ctx.dt, NULL));
154 PetscOptionsEnd();
155
156 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157 Initialize the solution
158 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159 PetscCall(VecGetArray(U, &u));
160 u[0] = PetscSqrtScalar(2.0);
161 u[1] = PetscSqrtScalar(3.0);
162 PetscCall(VecRestoreArray(U, &u));
163
164 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165 Create timestepping solver context
166 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
168 PetscCall(TSSetType(ts, TSARKIMEX));
169 PetscCall(TSARKIMEXSetFastSlowSplit(ts, PETSC_TRUE));
170
171 PetscCall(TSRHSSplitSetIS(ts, "slow", iss));
172 PetscCall(TSRHSSplitSetIS(ts, "fast", isf));
173 PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunctionFn *)RHSFunctionslow, &ctx));
174 PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunctionFn *)RHSFunctionfast, &ctx));
175 PetscCall(TSRHSSplitSetIFunction(ts, "fast", NULL, (TSIFunctionFn *)IFunctionfast, &ctx));
176 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
177 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
178 PetscCall(MatSetFromOptions(A));
179 PetscCall(MatSetUp(A));
180 PetscCall(TSRHSSplitSetIJacobian(ts, "fast", A, A, (TSIJacobianFn *)IJacobianfast, &ctx));
181
182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183 Set initial conditions
184 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185 PetscCall(TSSetSolution(ts, U));
186
187 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188 Set solver options
189 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190 PetscCall(TSSetMaxTime(ts, ctx.Tf));
191 PetscCall(TSSetTimeStep(ts, ctx.dt));
192 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
193 PetscCall(TSSetFromOptions(ts));
194
195 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196 Solve linear system
197 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198 PetscCall(TSSolve(ts, U));
199 PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
200
201 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202 Check the error of the PETSc solution
203 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204 PetscCall(TSGetTime(ts, &tt));
205 PetscCall(sol_true(tt, Utrue));
206 PetscCall(VecAXPY(Utrue, -1.0, U));
207 PetscCall(VecNorm(Utrue, NORM_2, &error));
208
209 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210 Print norm2 error
211 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
212 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error norm: %g\n", (double)PetscAbsScalar(error)));
213
214 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215 Free work space. All PETSc objects should be destroyed when they are no longer needed.
216 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
217 PetscCall(VecDestroy(&U));
218 PetscCall(TSDestroy(&ts));
219 PetscCall(VecDestroy(&Utrue));
220 PetscCall(ISDestroy(&iss));
221 PetscCall(ISDestroy(&isf));
222 PetscCall(PetscFree(indicess));
223 PetscCall(PetscFree(indicesf));
224 PetscCall(MatDestroy(&A));
225 PetscCall(PetscFinalize());
226 return 0;
227 }
228
229 /*TEST
230 build:
231 requires: !complex
232
233 test:
234
235 test:
236 suffix: 2
237 args: -ts_arkimex_initial_guess_extrapolate 1
238 output_file: output/ex3fastslowsplit_1.out
239
240 TEST*/
241