1 /*
2 Formatted test for TS routines.
3
4 Solves U_t=F(t,u)
5 Where:
6
7 [2*u1+u2 ]
8 F(t,u)= [u1+2*u2+u3]
9 [ u2+2*u3]
10
11 When run in parallel, each process solves the same set of equations separately.
12 */
13
14 static char help[] = "Solves a linear ODE. \n\n";
15
16 #include <petscts.h>
17 #include <petscpc.h>
18
19 extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
20 extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
21 extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
22 extern PetscErrorCode Initial(Vec, void *);
23 extern PetscErrorCode MyMatMult(Mat, Vec, Vec);
24
25 extern PetscReal solx(PetscReal);
26 extern PetscReal soly(PetscReal);
27 extern PetscReal solz(PetscReal);
28
main(int argc,char ** argv)29 int main(int argc, char **argv)
30 {
31 PetscInt time_steps = 100, steps;
32 Vec global;
33 PetscReal dt, ftime;
34 TS ts;
35 Mat A, S;
36 PetscBool nest = PETSC_FALSE;
37
38 PetscFunctionBeginUser;
39 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
40 PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL));
41 PetscCall(PetscOptionsGetBool(NULL, NULL, "-nest", &nest, NULL));
42
43 /* create vector to hold state */
44 if (nest) {
45 Vec g[3];
46
47 PetscCall(VecCreate(PETSC_COMM_WORLD, &g[0]));
48 PetscCall(VecSetSizes(g[0], 1, PETSC_DECIDE));
49 PetscCall(VecSetFromOptions(g[0]));
50 PetscCall(VecDuplicate(g[0], &g[1]));
51 PetscCall(VecDuplicate(g[0], &g[2]));
52 PetscCall(VecCreateNest(PETSC_COMM_WORLD, 3, NULL, g, &global));
53 PetscCall(VecDestroy(&g[0]));
54 PetscCall(VecDestroy(&g[1]));
55 PetscCall(VecDestroy(&g[2]));
56 } else {
57 PetscCall(VecCreate(PETSC_COMM_WORLD, &global));
58 PetscCall(VecSetSizes(global, 3, PETSC_DECIDE));
59 PetscCall(VecSetFromOptions(global));
60 }
61
62 /* set initial conditions */
63 PetscCall(Initial(global, NULL));
64
65 /* make timestep context */
66 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
67 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
68 PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
69 dt = 0.001;
70
71 /*
72 The user provides the RHS and Jacobian
73 */
74 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));
75 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
76 PetscCall(MatSetSizes(A, 3, 3, PETSC_DECIDE, PETSC_DECIDE));
77 PetscCall(MatSetFromOptions(A));
78 PetscCall(MatSetUp(A));
79 PetscCall(RHSJacobian(ts, 0.0, global, A, A, NULL));
80 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
81
82 PetscCall(MatCreateShell(PETSC_COMM_WORLD, 3, 3, PETSC_DECIDE, PETSC_DECIDE, NULL, &S));
83 PetscCall(MatShellSetOperation(S, MATOP_MULT, (PetscErrorCodeFn *)MyMatMult));
84 PetscCall(TSSetRHSJacobian(ts, S, A, RHSJacobian, NULL));
85
86 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
87 PetscCall(TSSetFromOptions(ts));
88
89 PetscCall(TSSetTimeStep(ts, dt));
90 PetscCall(TSSetMaxSteps(ts, time_steps));
91 PetscCall(TSSetMaxTime(ts, 1));
92 PetscCall(TSSetSolution(ts, global));
93
94 PetscCall(TSSolve(ts, global));
95 PetscCall(TSGetSolveTime(ts, &ftime));
96 PetscCall(TSGetStepNumber(ts, &steps));
97
98 /* free the memory */
99 PetscCall(TSDestroy(&ts));
100 PetscCall(VecDestroy(&global));
101 PetscCall(MatDestroy(&A));
102 PetscCall(MatDestroy(&S));
103
104 PetscCall(PetscFinalize());
105 return 0;
106 }
107
MyMatMult(Mat S,Vec x,Vec y)108 PetscErrorCode MyMatMult(Mat S, Vec x, Vec y)
109 {
110 const PetscScalar *inptr;
111 PetscScalar *outptr;
112
113 PetscFunctionBeginUser;
114 PetscCall(VecGetArrayRead(x, &inptr));
115 PetscCall(VecGetArrayWrite(y, &outptr));
116
117 outptr[0] = 2.0 * inptr[0] + inptr[1];
118 outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
119 outptr[2] = inptr[1] + 2.0 * inptr[2];
120
121 PetscCall(VecRestoreArrayRead(x, &inptr));
122 PetscCall(VecRestoreArrayWrite(y, &outptr));
123 PetscFunctionReturn(PETSC_SUCCESS);
124 }
125
Initial(Vec global,PetscCtx ctx)126 PetscErrorCode Initial(Vec global, PetscCtx ctx)
127 {
128 PetscScalar *localptr;
129
130 PetscFunctionBeginUser;
131 PetscCall(VecGetArrayWrite(global, &localptr));
132 localptr[0] = solx(0.0);
133 localptr[1] = soly(0.0);
134 localptr[2] = solz(0.0);
135 PetscCall(VecRestoreArrayWrite(global, &localptr));
136 PetscFunctionReturn(PETSC_SUCCESS);
137 }
138
Monitor(TS ts,PetscInt step,PetscReal time,Vec global,PetscCtx ctx)139 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, PetscCtx ctx)
140 {
141 const PetscScalar *tmp;
142 PetscScalar exact[] = {solx(time), soly(time), solz(time)};
143
144 PetscFunctionBeginUser;
145 PetscCall(VecGetArrayRead(global, &tmp));
146 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "At t =%14.6e u = %14.6e %14.6e %14.6e \n", (double)time, (double)PetscRealPart(tmp[0]), (double)PetscRealPart(tmp[1]), (double)PetscRealPart(tmp[2])));
147 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "At t =%14.6e errors = %14.6e %14.6e %14.6e \n", (double)time, (double)PetscRealPart(tmp[0] - exact[0]), (double)PetscRealPart(tmp[1] - exact[1]), (double)PetscRealPart(tmp[2] - exact[2])));
148 PetscCall(VecRestoreArrayRead(global, &tmp));
149 PetscFunctionReturn(PETSC_SUCCESS);
150 }
151
RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,PetscCtx ctx)152 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, PetscCtx ctx)
153 {
154 PetscScalar *outptr;
155 const PetscScalar *inptr;
156
157 PetscFunctionBeginUser;
158 /*Extract income array */
159 PetscCall(VecGetArrayRead(globalin, &inptr));
160
161 /* Extract outcome array*/
162 PetscCall(VecGetArrayWrite(globalout, &outptr));
163
164 outptr[0] = 2.0 * inptr[0] + inptr[1];
165 outptr[1] = inptr[0] + 2.0 * inptr[1] + inptr[2];
166 outptr[2] = inptr[1] + 2.0 * inptr[2];
167
168 PetscCall(VecRestoreArrayRead(globalin, &inptr));
169 PetscCall(VecRestoreArrayWrite(globalout, &outptr));
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172
RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,PetscCtx ctx)173 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, PetscCtx ctx)
174 {
175 PetscScalar v[3];
176 PetscInt idx[3], rst;
177
178 PetscFunctionBeginUser;
179 PetscCall(VecGetOwnershipRange(x, &rst, NULL));
180 idx[0] = 0 + rst;
181 idx[1] = 1 + rst;
182 idx[2] = 2 + rst;
183
184 v[0] = 2.0;
185 v[1] = 1.0;
186 v[2] = 0.0;
187 PetscCall(MatSetValues(BB, 1, idx, 3, idx, v, INSERT_VALUES));
188
189 v[0] = 1.0;
190 v[1] = 2.0;
191 v[2] = 1.0;
192 PetscCall(MatSetValues(BB, 1, idx + 1, 3, idx, v, INSERT_VALUES));
193
194 v[0] = 0.0;
195 v[1] = 1.0;
196 v[2] = 2.0;
197 PetscCall(MatSetValues(BB, 1, idx + 2, 3, idx, v, INSERT_VALUES));
198
199 PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY));
200 PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY));
201
202 if (A != BB) {
203 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
204 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
205 }
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208
209 /*
210 The exact solutions
211 */
solx(PetscReal t)212 PetscReal solx(PetscReal t)
213 {
214 return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0)) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0));
215 }
216
soly(PetscReal t)217 PetscReal soly(PetscReal t)
218 {
219 return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / PetscSqrtReal(2.0) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / PetscSqrtReal(2.0);
220 }
221
solz(PetscReal t)222 PetscReal solz(PetscReal t)
223 {
224 return PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / 2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0)) + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / 2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0)) * t) / (2.0 * PetscSqrtReal(2.0));
225 }
226
227 /*TEST
228
229 test:
230 suffix: euler
231 args: -ts_type euler -nest {{0 1}}
232 requires: !single
233
234 test:
235 suffix: beuler
236 args: -ts_type beuler -nest {{0 1}}
237 requires: !single
238
239 test:
240 suffix: rk
241 args: -ts_type rk -nest {{0 1}} -ts_adapt_monitor
242 requires: !single
243
244 test:
245 diff_args: -j
246 requires: double !complex
247 output_file: output/ex2_be_adapt.out
248 suffix: bdf_1_adapt
249 args: -ts_type bdf -ts_bdf_order 1 -ts_adapt_type basic -ts_adapt_clip 0,2
250
251 test:
252 diff_args: -j
253 requires: double !complex
254 suffix: be_adapt
255 args: -ts_type beuler -ts_adapt_type basic -ts_adapt_clip 0,2
256
257 TEST*/
258