1 static char help[] = "Testing integrators on the simple harmonic oscillator\n";
2
3 /*
4 In order to check long time behavior, you can give
5
6 -ts_max_steps 10000 -ts_max_time 10.0 -output_step 1000
7
8 To look at the long time behavior for different integrators, we can use the nergy monitor. Below we see that Euler is almost 100 times worse at conserving energy than the symplectic integrator of the same order.
9
10 make -f ./gmakefile test search="dm_impls_swarm_tests-ex4_bsi_1"
11 EXTRA_OPTIONS="-ts_max_steps 10000 -ts_max_time 10.0 -output_step 1000 -ts_type euler" | grep error
12
13 energy/exact energy 398.236 / 396.608 (0.4104399231)
14 energy/exact energy 1579.52 / 1573.06 (0.4104399231)
15 energy/exact energy 397.421 / 396.608 (0.2050098479)
16 energy/exact energy 1576.29 / 1573.06 (0.2050098479)
17 energy/exact energy 397.014 / 396.608 (0.1024524454)
18 energy/exact energy 1574.68 / 1573.06 (0.1024524454)
19
20 make -f ./gmakefile test search="dm_impls_swarm_tests-ex4_bsi_1"
21 EXTRA_OPTIONS="-ts_max_steps 10000 -ts_max_time 10.0 -output_step 1000" | grep error
22
23 energy/exact energy 396.579 / 396.608 (0.0074080434)
24 energy/exact energy 1572.95 / 1573.06 (0.0074080434)
25 energy/exact energy 396.593 / 396.608 (0.0037040885)
26 energy/exact energy 1573.01 / 1573.06 (0.0037040886)
27 energy/exact energy 396.601 / 396.608 (0.0018520613)
28 energy/exact energy 1573.04 / 1573.06 (0.0018520613)
29
30 We can look at third order integrators in the same way, but we need to use more steps.
31
32 make -f ./gmakefile test search="dm_impls_swarm_tests-ex4_bsi_1"
33 EXTRA_OPTIONS="-ts_max_steps 1000000 -ts_max_time 1000.0 -output_step 100000 -ts_type rk -ts_adapt_type none" | grep error
34
35 energy/exact energy 396.608 / 396.608 (0.0000013981)
36 energy/exact energy 1573.06 / 1573.06 (0.0000013981)
37 energy/exact energy 396.608 / 396.608 (0.0000001747)
38 energy/exact energy 1573.06 / 1573.06 (0.0000001748)
39 energy/exact energy 396.608 / 396.608 (0.0000000218)
40 energy/exact energy 1573.06 / 1573.06 (0.0000000218)
41
42 make -f ./gmakefile test search="dm_impls_swarm_tests-ex4_bsi_3"
43 EXTRA_OPTIONS="-ts_max_steps 1000000 -ts_max_time 1000.0 -output_step 100000 -ts_adapt_type none" | grep error
44
45 energy/exact energy 396.608 / 396.608 (0.0000000007)
46 energy/exact energy 1573.06 / 1573.06 (0.0000000007)
47 energy/exact energy 396.608 / 396.608 (0.0000000001)
48 energy/exact energy 1573.06 / 1573.06 (0.0000000001)
49 energy/exact energy 396.608 / 396.608 (0.0000000000)
50 energy/exact energy 1573.06 / 1573.06 (0.0000000000)
51 */
52
53 #include <petscts.h>
54 #include <petscdmplex.h>
55 #include <petscdmswarm.h>
56 #include <petsc/private/dmpleximpl.h> /* For norm and dot */
57
58 typedef struct {
59 PetscReal omega; /* Oscillation frequency omega */
60 PetscBool error; /* Flag for printing the error */
61 PetscInt ostep; /* print the energy at each ostep time steps */
62 } AppCtx;
63
ProcessOptions(MPI_Comm comm,AppCtx * options)64 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
65 {
66 PetscFunctionBeginUser;
67 options->omega = 64.0;
68 options->error = PETSC_FALSE;
69 options->ostep = 100;
70
71 PetscOptionsBegin(comm, "", "Harmonic Oscillator Options", "DMSWARM");
72 PetscCall(PetscOptionsReal("-omega", "Oscillator frequency", "ex4.c", options->omega, &options->omega, NULL));
73 PetscCall(PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL));
74 PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, NULL));
75 PetscOptionsEnd();
76 PetscFunctionReturn(PETSC_SUCCESS);
77 }
78
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)79 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
80 {
81 PetscFunctionBeginUser;
82 PetscCall(DMCreate(comm, dm));
83 PetscCall(DMSetType(*dm, DMPLEX));
84 PetscCall(DMSetFromOptions(*dm));
85 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
86 PetscFunctionReturn(PETSC_SUCCESS);
87 }
88
CreateSwarm(DM dm,AppCtx * user,DM * sw)89 static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
90 {
91 PetscReal v0[1] = {0.}; // Initialize velocities to 0
92 PetscInt dim;
93
94 PetscFunctionBeginUser;
95 PetscCall(DMGetDimension(dm, &dim));
96 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
97 PetscCall(DMSetType(*sw, DMSWARM));
98 PetscCall(DMSetDimension(*sw, dim));
99 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
100 PetscCall(DMSwarmSetCellDM(*sw, dm));
101 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
102 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
103 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
104 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initCoordinates", dim, PETSC_REAL));
105 PetscCall(DMSwarmFinalizeFieldRegister(*sw));
106 PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw));
107 PetscCall(DMSwarmInitializeCoordinates(*sw));
108 PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, v0));
109 PetscCall(DMSetFromOptions(*sw));
110 PetscCall(DMSetApplicationContext(*sw, user));
111 PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
112 PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
113 {
114 Vec gc, gc0;
115
116 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
117 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initCoordinates", &gc0));
118 PetscCall(VecCopy(gc, gc0));
119 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
120 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initCoordinates", &gc0));
121 }
122 {
123 const char *fieldnames[2] = {DMSwarmPICField_coor, "velocity"};
124 PetscCall(DMSwarmVectorDefineFields(*sw, 2, fieldnames));
125 }
126 PetscFunctionReturn(PETSC_SUCCESS);
127 }
128
RHSFunction(TS ts,PetscReal t,Vec U,Vec G,PetscCtx ctx)129 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, PetscCtx ctx)
130 {
131 const PetscReal omega = ((AppCtx *)ctx)->omega;
132 DM sw;
133 const PetscScalar *u;
134 PetscScalar *g;
135 PetscInt dim, d, Np, p;
136
137 PetscFunctionBeginUser;
138 PetscCall(TSGetDM(ts, &sw));
139 PetscCall(DMGetDimension(sw, &dim));
140 PetscCall(VecGetLocalSize(U, &Np));
141 PetscCall(VecGetArray(G, &g));
142 PetscCall(VecGetArrayRead(U, &u));
143 Np /= 2 * dim;
144 for (p = 0; p < Np; ++p) {
145 for (d = 0; d < dim; ++d) {
146 g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
147 g[(p * 2 + 1) * dim + d] = -PetscSqr(omega) * u[(p * 2 + 0) * dim + d];
148 }
149 }
150 PetscCall(VecRestoreArrayRead(U, &u));
151 PetscCall(VecRestoreArray(G, &g));
152 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154
155 /* J_{ij} = dF_i/dx_j
156 J_p = ( 0 1)
157 (-w^2 0)
158 */
RHSJacobian(TS ts,PetscReal t,Vec U,Mat J,Mat P,PetscCtx ctx)159 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, PetscCtx ctx)
160 {
161 PetscScalar vals[4] = {0., 1., -PetscSqr(((AppCtx *)ctx)->omega), 0.};
162 DM sw;
163 PetscInt dim, d, Np, p, rStart;
164
165 PetscFunctionBeginUser;
166 PetscCall(TSGetDM(ts, &sw));
167 PetscCall(DMGetDimension(sw, &dim));
168 PetscCall(VecGetLocalSize(U, &Np));
169 PetscCall(MatGetOwnershipRange(J, &rStart, NULL));
170 Np /= 2 * dim;
171 for (p = 0; p < Np; ++p) {
172 for (d = 0; d < dim; ++d) {
173 const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
174 PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES));
175 }
176 }
177 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
178 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
179 PetscFunctionReturn(PETSC_SUCCESS);
180 }
181
RHSFunctionX(TS ts,PetscReal t,Vec V,Vec Xres,PetscCtx ctx)182 static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, PetscCtx ctx)
183 {
184 const PetscScalar *v;
185 PetscScalar *xres;
186 PetscInt Np, p;
187
188 PetscFunctionBeginUser;
189 PetscCall(VecGetLocalSize(Xres, &Np));
190 PetscCall(VecGetArrayRead(V, &v));
191 PetscCall(VecGetArray(Xres, &xres));
192 for (p = 0; p < Np; ++p) xres[p] = v[p];
193 PetscCall(VecRestoreArrayRead(V, &v));
194 PetscCall(VecRestoreArray(Xres, &xres));
195 PetscFunctionReturn(PETSC_SUCCESS);
196 }
197
RHSFunctionV(TS ts,PetscReal t,Vec X,Vec Vres,PetscCtx ctx)198 static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, PetscCtx ctx)
199 {
200 const PetscReal omega = ((AppCtx *)ctx)->omega;
201 const PetscScalar *x;
202 PetscScalar *vres;
203 PetscInt Np, p;
204
205 PetscFunctionBeginUser;
206 PetscCall(VecGetArray(Vres, &vres));
207 PetscCall(VecGetArrayRead(X, &x));
208 PetscCall(VecGetLocalSize(Vres, &Np));
209 for (p = 0; p < Np; ++p) vres[p] = -PetscSqr(omega) * x[p];
210 PetscCall(VecRestoreArrayRead(X, &x));
211 PetscCall(VecRestoreArray(Vres, &vres));
212 PetscFunctionReturn(PETSC_SUCCESS);
213 }
214
RHSJacobianS(TS ts,PetscReal t,Vec U,Mat S,PetscCtx ctx)215 PetscErrorCode RHSJacobianS(TS ts, PetscReal t, Vec U, Mat S, PetscCtx ctx)
216 {
217 PetscScalar vals[4] = {0., 1., -1., 0.};
218 DM sw;
219 PetscInt dim, d, Np, p, rStart;
220
221 PetscFunctionBeginUser;
222 PetscCall(TSGetDM(ts, &sw));
223 PetscCall(DMGetDimension(sw, &dim));
224 PetscCall(VecGetLocalSize(U, &Np));
225 PetscCall(MatGetOwnershipRange(S, &rStart, NULL));
226 Np /= 2 * dim;
227 for (p = 0; p < Np; ++p) {
228 for (d = 0; d < dim; ++d) {
229 const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
230 PetscCall(MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES));
231 }
232 }
233 PetscCall(MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY));
234 PetscCall(MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY));
235 PetscFunctionReturn(PETSC_SUCCESS);
236 }
237
RHSObjectiveF(TS ts,PetscReal t,Vec U,PetscScalar * F,PetscCtx ctx)238 PetscErrorCode RHSObjectiveF(TS ts, PetscReal t, Vec U, PetscScalar *F, PetscCtx ctx)
239 {
240 const PetscReal omega = ((AppCtx *)ctx)->omega;
241 DM sw;
242 const PetscScalar *u;
243 PetscInt dim, Np, p;
244
245 PetscFunctionBeginUser;
246 PetscCall(TSGetDM(ts, &sw));
247 PetscCall(DMGetDimension(sw, &dim));
248 PetscCall(VecGetArrayRead(U, &u));
249 PetscCall(VecGetLocalSize(U, &Np));
250 Np /= 2 * dim;
251 for (p = 0; p < Np; ++p) {
252 const PetscReal x2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 0) * dim], &u[(p * 2 + 0) * dim]);
253 const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
254 const PetscReal E = 0.5 * (v2 + PetscSqr(omega) * x2);
255
256 *F += E;
257 }
258 PetscCall(VecRestoreArrayRead(U, &u));
259 PetscFunctionReturn(PETSC_SUCCESS);
260 }
261
262 /* dF/dx = omega^2 x dF/dv = v */
RHSFunctionG(TS ts,PetscReal t,Vec U,Vec G,PetscCtx ctx)263 PetscErrorCode RHSFunctionG(TS ts, PetscReal t, Vec U, Vec G, PetscCtx ctx)
264 {
265 const PetscReal omega = ((AppCtx *)ctx)->omega;
266 DM sw;
267 const PetscScalar *u;
268 PetscScalar *g;
269 PetscInt dim, d, Np, p;
270
271 PetscFunctionBeginUser;
272 PetscCall(TSGetDM(ts, &sw));
273 PetscCall(DMGetDimension(sw, &dim));
274 PetscCall(VecGetArray(G, &g));
275 PetscCall(VecGetArrayRead(U, &u));
276 PetscCall(VecGetLocalSize(U, &Np));
277 Np /= 2 * dim;
278 for (p = 0; p < Np; ++p) {
279 for (d = 0; d < dim; ++d) {
280 g[(p * 2 + 0) * dim + d] = PetscSqr(omega) * u[(p * 2 + 0) * dim + d];
281 g[(p * 2 + 1) * dim + d] = u[(p * 2 + 1) * dim + d];
282 }
283 }
284 PetscCall(VecRestoreArrayRead(U, &u));
285 PetscCall(VecRestoreArray(G, &g));
286 PetscFunctionReturn(PETSC_SUCCESS);
287 }
288
CreateSolution(TS ts)289 static PetscErrorCode CreateSolution(TS ts)
290 {
291 DM sw;
292 Vec u;
293 PetscInt dim, Np;
294
295 PetscFunctionBegin;
296 PetscCall(TSGetDM(ts, &sw));
297 PetscCall(DMGetDimension(sw, &dim));
298 PetscCall(DMSwarmGetLocalSize(sw, &Np));
299 PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
300 PetscCall(VecSetBlockSize(u, dim));
301 PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE));
302 PetscCall(VecSetUp(u));
303 PetscCall(TSSetSolution(ts, u));
304 PetscCall(VecDestroy(&u));
305 PetscFunctionReturn(PETSC_SUCCESS);
306 }
307
SetProblem(TS ts)308 static PetscErrorCode SetProblem(TS ts)
309 {
310 AppCtx *user;
311 DM sw;
312
313 PetscFunctionBegin;
314 PetscCall(TSGetDM(ts, &sw));
315 PetscCall(DMGetApplicationContext(sw, &user));
316 // Define unified system for (X, V)
317 {
318 Mat J;
319 PetscInt dim, Np;
320
321 PetscCall(DMGetDimension(sw, &dim));
322 PetscCall(DMSwarmGetLocalSize(sw, &Np));
323 PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
324 PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE));
325 PetscCall(MatSetBlockSize(J, 2 * dim));
326 PetscCall(MatSetFromOptions(J));
327 PetscCall(MatSetUp(J));
328 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user));
329 PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user));
330 PetscCall(MatDestroy(&J));
331 }
332 // Define split system for X and V
333 {
334 IS isx, isv, istmp;
335 const PetscInt *idx;
336 PetscInt dim, Np;
337
338 PetscCall(DMGetDimension(sw, &dim));
339 PetscCall(DMSwarmGetLocalSize(sw, &Np));
340 PetscCall(ISCreateStride(PETSC_COMM_SELF, Np, 0, 2, &istmp));
341 PetscCall(ISGetIndices(istmp, &idx));
342 PetscCall(ISCreateBlock(PETSC_COMM_SELF, dim, Np, idx, PETSC_COPY_VALUES, &isx));
343 PetscCall(ISRestoreIndices(istmp, &idx));
344 PetscCall(ISDestroy(&istmp));
345 PetscCall(ISCreateStride(PETSC_COMM_SELF, Np, 1, 2, &istmp));
346 PetscCall(ISGetIndices(istmp, &idx));
347 PetscCall(ISCreateBlock(PETSC_COMM_SELF, dim, Np, idx, PETSC_COPY_VALUES, &isv));
348 PetscCall(ISRestoreIndices(istmp, &idx));
349 PetscCall(ISDestroy(&istmp));
350 PetscCall(TSRHSSplitSetIS(ts, "position", isx));
351 PetscCall(TSRHSSplitSetIS(ts, "momentum", isv));
352 PetscCall(ISDestroy(&isx));
353 PetscCall(ISDestroy(&isv));
354 PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user));
355 PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user));
356 }
357 // Define symplectic formulation U_t = S . G, where G = grad F
358 {
359 PetscCall(TSDiscGradSetFormulation(ts, RHSJacobianS, RHSObjectiveF, RHSFunctionG, user));
360 }
361 PetscFunctionReturn(PETSC_SUCCESS);
362 }
363
InitializeSolve(TS ts,Vec u)364 static PetscErrorCode InitializeSolve(TS ts, Vec u)
365 {
366 DM sw;
367 Vec gc, gc0;
368 IS isx, isv;
369 AppCtx *user;
370
371 PetscFunctionBeginUser;
372 PetscCall(TSGetDM(ts, &sw));
373 PetscCall(DMGetApplicationContext(sw, &user));
374 {
375 PetscReal v0[1] = {1.};
376
377 PetscCall(DMSwarmInitializeCoordinates(sw));
378 PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0));
379 PetscCall(SetProblem(ts));
380 }
381 PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
382 PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
383 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
384 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0));
385 PetscCall(VecCopy(gc, gc0));
386 PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc));
387 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
388 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0));
389 PetscCall(VecISSet(u, isv, 0.));
390 PetscFunctionReturn(PETSC_SUCCESS);
391 }
392
ComputeError(TS ts,Vec U,Vec E)393 static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
394 {
395 MPI_Comm comm;
396 DM sw;
397 AppCtx *user;
398 const PetscScalar *u;
399 const PetscReal *coords;
400 PetscScalar *e;
401 PetscReal t;
402 PetscInt dim, d, Np, p;
403
404 PetscFunctionBeginUser;
405 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
406 PetscCall(TSGetDM(ts, &sw));
407 PetscCall(DMGetApplicationContext(sw, &user));
408 PetscCall(DMGetDimension(sw, &dim));
409 PetscCall(TSGetSolveTime(ts, &t));
410 PetscCall(VecGetArray(E, &e));
411 PetscCall(VecGetArrayRead(U, &u));
412 PetscCall(VecGetLocalSize(U, &Np));
413 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
414 Np /= 2 * dim;
415 for (p = 0; p < Np; ++p) {
416 const PetscReal omega = user->omega;
417 const PetscReal ct = PetscCosReal(omega * t);
418 const PetscReal st = PetscSinReal(omega * t);
419 const PetscReal x0 = DMPlex_NormD_Internal(dim, &coords[p * dim]);
420 const PetscReal ex = x0 * ct;
421 const PetscReal ev = -x0 * omega * st;
422 const PetscReal x = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 0) * dim], &coords[p * dim]) / x0;
423 const PetscReal v = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &coords[p * dim]) / x0;
424
425 if (user->error) {
426 const PetscReal en = 0.5 * (v * v + PetscSqr(omega) * x * x);
427 const PetscReal exen = 0.5 * PetscSqr(omega * x0);
428 PetscCall(PetscPrintf(comm, "p%" PetscInt_FMT " error [%.2f %.2f] sol [%.6lf %.6lf] exact [%.6lf %.6lf] energy/exact energy %g / %g (%.10lf%%)\n", p, (double)PetscAbsReal(x - ex), (double)PetscAbsReal(v - ev), (double)x, (double)v, (double)ex, (double)ev, (double)en, (double)exen, (double)(PetscAbsReal(exen - en) * 100. / exen)));
429 }
430 for (d = 0; d < dim; ++d) {
431 e[(p * 2 + 0) * dim + d] = u[(p * 2 + 0) * dim + d] - coords[p * dim + d] * ct;
432 e[(p * 2 + 1) * dim + d] = u[(p * 2 + 1) * dim + d] + coords[p * dim + d] * omega * st;
433 }
434 }
435 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
436 PetscCall(VecRestoreArrayRead(U, &u));
437 PetscCall(VecRestoreArray(E, &e));
438 PetscFunctionReturn(PETSC_SUCCESS);
439 }
440
EnergyMonitor(TS ts,PetscInt step,PetscReal t,Vec U,PetscCtx ctx)441 static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, PetscCtx ctx)
442 {
443 const PetscReal omega = ((AppCtx *)ctx)->omega;
444 const PetscInt ostep = ((AppCtx *)ctx)->ostep;
445 DM sw;
446 const PetscScalar *u;
447 PetscReal dt;
448 PetscInt dim, Np, p;
449 MPI_Comm comm;
450
451 PetscFunctionBeginUser;
452 if (step % ostep == 0) {
453 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
454 PetscCall(TSGetDM(ts, &sw));
455 PetscCall(TSGetTimeStep(ts, &dt));
456 PetscCall(DMGetDimension(sw, &dim));
457 PetscCall(VecGetArrayRead(U, &u));
458 PetscCall(VecGetLocalSize(U, &Np));
459 Np /= 2 * dim;
460 if (!step) PetscCall(PetscPrintf(comm, "Time Step Part Energy Mod Energy\n"));
461 for (p = 0; p < Np; ++p) {
462 const PetscReal x2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 0) * dim], &u[(p * 2 + 0) * dim]);
463 const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
464 const PetscReal E = 0.5 * (v2 + PetscSqr(omega) * x2);
465 const PetscReal mE = 0.5 * (v2 + PetscSqr(omega) * x2 - PetscSqr(omega) * dt * PetscSqrtReal(x2 * v2));
466
467 PetscCall(PetscPrintf(comm, "%.6lf %4" PetscInt_FMT " %4" PetscInt_FMT " %10.4lf %10.4lf\n", (double)t, step, p, (double)E, (double)mE));
468 }
469 PetscCall(VecRestoreArrayRead(U, &u));
470 }
471 PetscFunctionReturn(PETSC_SUCCESS);
472 }
473
main(int argc,char ** argv)474 int main(int argc, char **argv)
475 {
476 DM dm, sw;
477 TS ts;
478 Vec u;
479 AppCtx user;
480
481 PetscFunctionBeginUser;
482 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
483 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
484 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
485 PetscCall(CreateSwarm(dm, &user, &sw));
486 PetscCall(DMSetApplicationContext(sw, &user));
487
488 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
489 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
490 PetscCall(TSSetDM(ts, sw));
491 PetscCall(TSSetMaxTime(ts, 0.1));
492 PetscCall(TSSetTimeStep(ts, 0.00001));
493 PetscCall(TSSetMaxSteps(ts, 100));
494 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
495 PetscCall(TSMonitorSet(ts, EnergyMonitor, &user, NULL));
496 PetscCall(TSSetFromOptions(ts));
497
498 PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
499 PetscCall(TSSetComputeExactError(ts, ComputeError));
500
501 PetscCall(CreateSolution(ts));
502 PetscCall(TSGetSolution(ts, &u));
503 PetscCall(TSComputeInitialCondition(ts, u));
504 PetscCall(TSSolve(ts, NULL));
505
506 PetscCall(TSDestroy(&ts));
507 PetscCall(DMDestroy(&sw));
508 PetscCall(DMDestroy(&dm));
509 PetscCall(PetscFinalize());
510 return 0;
511 }
512
513 /*TEST
514
515 build:
516 requires: triangle !single !complex
517
518 testset:
519 args: -dm_plex_dim 1 -dm_plex_box_faces 1 -dm_plex_box_lower -1 -dm_plex_box_upper 1 \
520 -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \
521 -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
522 -dm_view -output_step 50 -error
523 test:
524 suffix: bsi_1
525 args: -ts_basicsymplectic_type 1
526 test:
527 suffix: bsi_2
528 args: -ts_basicsymplectic_type 2
529 test:
530 suffix: bsi_3
531 args: -ts_basicsymplectic_type 3
532 test:
533 suffix: bsi_4
534 args: -ts_basicsymplectic_type 4 -ts_time_step 0.0001
535
536 testset:
537 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1,-1 -dm_plex_box_upper 1,1 \
538 -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \
539 -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
540 -dm_view -output_step 50 -error
541 test:
542 suffix: bsi_2d_1
543 args: -ts_basicsymplectic_type 1
544 test:
545 suffix: bsi_2d_2
546 args: -ts_basicsymplectic_type 2
547 test:
548 suffix: bsi_2d_3
549 args: -ts_basicsymplectic_type 3
550 test:
551 suffix: bsi_2d_4
552 args: -ts_basicsymplectic_type 4 -ts_time_step 0.0001
553
554 testset:
555 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1,-1,-1 -dm_plex_box_upper 1,1,1 \
556 -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \
557 -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
558 -dm_view -output_step 50 -error
559 test:
560 suffix: bsi_3d_1
561 args: -ts_basicsymplectic_type 1
562 test:
563 suffix: bsi_3d_2
564 args: -ts_basicsymplectic_type 2
565 test:
566 suffix: bsi_3d_3
567 args: -ts_basicsymplectic_type 3
568 test:
569 suffix: bsi_3d_4
570 args: -ts_basicsymplectic_type 4 -ts_time_step 0.0001
571
572 testset:
573 args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \
574 -ts_type theta -ts_theta_theta 0.5 -ts_convergence_estimate -convest_num_refine 2 \
575 -mat_type baij -ksp_error_if_not_converged -pc_type lu \
576 -dm_view -output_step 50 -error
577 test:
578 suffix: im_1d
579 args: -dm_plex_dim 1 -dm_plex_box_faces 1 -dm_plex_box_lower -1 -dm_plex_box_upper 1
580 test:
581 suffix: im_2d
582 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1,-1 -dm_plex_box_upper 1,1
583 test:
584 suffix: im_3d
585 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1,-1,-1 -dm_plex_box_upper 1,1,1
586
587 testset:
588 args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \
589 -ts_type discgrad -ts_discgrad_type none -ts_convergence_estimate -convest_num_refine 2 \
590 -mat_type baij -ksp_error_if_not_converged -pc_type lu \
591 -dm_view -output_step 50 -error
592 test:
593 suffix: dg_1d
594 args: -dm_plex_dim 1 -dm_plex_box_faces 1 -dm_plex_box_lower -1 -dm_plex_box_upper 1
595 test:
596 suffix: dg_2d
597 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -1,-1 -dm_plex_box_upper 1,1
598 test:
599 suffix: dg_3d
600 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 -dm_plex_box_lower -1,-1,-1 -dm_plex_box_upper 1,1,1
601
602 testset:
603 args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_density constant \
604 -ts_type discgrad -ts_convergence_estimate -convest_num_refine 2 \
605 -mat_type baij -ksp_error_if_not_converged -pc_type lu \
606 -dm_view -output_step 50 -error
607 test:
608 suffix: dg_gonzalez
609 args: -dm_plex_dim 1 -dm_plex_box_faces 1 -dm_plex_box_lower -1 -dm_plex_box_upper 1 -ts_discgrad_type gonzalez
610 test:
611 suffix: dg_average
612 args: -dm_plex_dim 1 -dm_plex_box_faces 1 -dm_plex_box_lower -1 -dm_plex_box_upper 1 -ts_discgrad_type average
613
614 TEST*/
615