1 static char help[] = "Vlasov example of central orbits\n";
2
3 /*
4 To visualize the orbit, we can used
5
6 -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain -1 -ts_monitor_sp_swarm_phase 0 -draw_size 500,500
7
8 and we probably want it to run fast and not check convergence
9
10 -convest_num_refine 0 -ts_time_step 0.01 -ts_max_steps 100 -ts_max_time 100 -output_step 10
11 */
12
13 #include <petscts.h>
14 #include <petscdmplex.h>
15 #include <petscdmswarm.h>
16 #include <petsc/private/dmpleximpl.h> /* For norm and dot */
17
18 PETSC_EXTERN PetscErrorCode circleSingleX(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
19 PETSC_EXTERN PetscErrorCode circleSingleV(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
20 PETSC_EXTERN PetscErrorCode circleMultipleX(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
21 PETSC_EXTERN PetscErrorCode circleMultipleV(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
22
23 typedef struct {
24 PetscBool error; /* Flag for printing the error */
25 PetscInt ostep; /* print the energy at each ostep time steps */
26 } AppCtx;
27
ProcessOptions(MPI_Comm comm,AppCtx * options)28 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
29 {
30 PetscFunctionBeginUser;
31 options->error = PETSC_FALSE;
32 options->ostep = 100;
33
34 PetscOptionsBegin(comm, "", "Central Orbit Options", "DMSWARM");
35 PetscCall(PetscOptionsBool("-error", "Flag to print the error", "ex5.c", options->error, &options->error, NULL));
36 PetscCall(PetscOptionsInt("-output_step", "Number of time steps between output", "ex5.c", options->ostep, &options->ostep, NULL));
37 PetscOptionsEnd();
38 PetscFunctionReturn(PETSC_SUCCESS);
39 }
40
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)41 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
42 {
43 PetscFunctionBeginUser;
44 PetscCall(DMCreate(comm, dm));
45 PetscCall(DMSetType(*dm, DMPLEX));
46 PetscCall(DMSetFromOptions(*dm));
47 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
48 PetscFunctionReturn(PETSC_SUCCESS);
49 }
50
CreateSwarm(DM dm,AppCtx * user,DM * sw)51 static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
52 {
53 PetscReal v0[1] = {1.};
54 PetscInt dim;
55
56 PetscFunctionBeginUser;
57 PetscCall(DMGetDimension(dm, &dim));
58 PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), sw));
59 PetscCall(DMSetType(*sw, DMSWARM));
60 PetscCall(DMSetDimension(*sw, dim));
61 PetscCall(DMSwarmSetType(*sw, DMSWARM_PIC));
62 PetscCall(DMSwarmSetCellDM(*sw, dm));
63 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR));
64 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL));
65 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT));
66 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initCoordinates", dim, PETSC_REAL));
67 PetscCall(DMSwarmRegisterPetscDatatypeField(*sw, "initVelocity", dim, PETSC_REAL));
68 PetscCall(DMSwarmFinalizeFieldRegister(*sw));
69 PetscCall(DMSwarmComputeLocalSizeFromOptions(*sw));
70 PetscCall(DMSwarmInitializeCoordinates(*sw));
71 PetscCall(DMSwarmInitializeVelocitiesFromOptions(*sw, v0));
72 PetscCall(DMSetFromOptions(*sw));
73 PetscCall(DMSetApplicationContext(*sw, user));
74 PetscCall(PetscObjectSetName((PetscObject)*sw, "Particles"));
75 PetscCall(DMViewFromOptions(*sw, NULL, "-sw_view"));
76 {
77 Vec gc, gc0, gv, gv0;
78
79 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
80 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initCoordinates", &gc0));
81 PetscCall(VecCopy(gc, gc0));
82 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc));
83 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initCoordinates", &gc0));
84 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "velocity", &gv));
85 PetscCall(DMSwarmCreateGlobalVectorFromField(*sw, "initVelocity", &gv0));
86 PetscCall(VecCopy(gv, gv0));
87 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "velocity", &gv));
88 PetscCall(DMSwarmDestroyGlobalVectorFromField(*sw, "initVelocity", &gv0));
89 }
90 PetscFunctionReturn(PETSC_SUCCESS);
91 }
92
RHSFunction(TS ts,PetscReal t,Vec U,Vec G,PetscCtx ctx)93 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, PetscCtx ctx)
94 {
95 DM sw;
96 const PetscReal *coords, *vel;
97 const PetscScalar *u;
98 PetscScalar *g;
99 PetscInt dim, d, Np, p;
100
101 PetscFunctionBeginUser;
102 PetscCall(TSGetDM(ts, &sw));
103 PetscCall(DMGetDimension(sw, &dim));
104 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
105 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
106 PetscCall(VecGetLocalSize(U, &Np));
107 PetscCall(VecGetArrayRead(U, &u));
108 PetscCall(VecGetArray(G, &g));
109 Np /= 2 * dim;
110 for (p = 0; p < Np; ++p) {
111 const PetscReal x0 = coords[p * dim + 0];
112 const PetscReal vy0 = vel[p * dim + 1];
113 const PetscReal omega = vy0 / x0;
114
115 for (d = 0; d < dim; ++d) {
116 g[(p * 2 + 0) * dim + d] = u[(p * 2 + 1) * dim + d];
117 g[(p * 2 + 1) * dim + d] = -PetscSqr(omega) * u[(p * 2 + 0) * dim + d];
118 }
119 }
120 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
121 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
122 PetscCall(VecRestoreArrayRead(U, &u));
123 PetscCall(VecRestoreArray(G, &g));
124 PetscFunctionReturn(PETSC_SUCCESS);
125 }
126
127 /* J_{ij} = dF_i/dx_j
128 J_p = ( 0 1)
129 (-w^2 0)
130 */
RHSJacobian(TS ts,PetscReal t,Vec U,Mat J,Mat P,PetscCtx ctx)131 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, PetscCtx ctx)
132 {
133 DM sw;
134 const PetscReal *coords, *vel;
135 PetscInt dim, d, Np, p, rStart;
136
137 PetscFunctionBeginUser;
138 PetscCall(TSGetDM(ts, &sw));
139 PetscCall(DMGetDimension(sw, &dim));
140 PetscCall(VecGetLocalSize(U, &Np));
141 PetscCall(MatGetOwnershipRange(J, &rStart, NULL));
142 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
143 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
144 Np /= 2 * dim;
145 for (p = 0; p < Np; ++p) {
146 const PetscReal x0 = coords[p * dim + 0];
147 const PetscReal vy0 = vel[p * dim + 1];
148 const PetscReal omega = vy0 / x0;
149 PetscScalar vals[4] = {0., 1., -PetscSqr(omega), 0.};
150
151 for (d = 0; d < dim; ++d) {
152 const PetscInt rows[2] = {(p * 2 + 0) * dim + d + rStart, (p * 2 + 1) * dim + d + rStart};
153 PetscCall(MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES));
154 }
155 }
156 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
157 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
158 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
159 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
160 PetscFunctionReturn(PETSC_SUCCESS);
161 }
162
RHSFunctionX(TS ts,PetscReal t,Vec V,Vec Xres,PetscCtx ctx)163 static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, PetscCtx ctx)
164 {
165 DM sw;
166 const PetscScalar *v;
167 PetscScalar *xres;
168 PetscInt Np, p, dim, d;
169
170 PetscFunctionBeginUser;
171 PetscCall(TSGetDM(ts, &sw));
172 PetscCall(DMGetDimension(sw, &dim));
173 PetscCall(VecGetLocalSize(Xres, &Np));
174 PetscCall(VecGetArrayRead(V, &v));
175 PetscCall(VecGetArray(Xres, &xres));
176 Np /= dim;
177 for (p = 0; p < Np; ++p)
178 for (d = 0; d < dim; ++d) xres[p * dim + d] = v[p * dim + d];
179 PetscCall(VecRestoreArrayRead(V, &v));
180 PetscCall(VecRestoreArray(Xres, &xres));
181 PetscFunctionReturn(PETSC_SUCCESS);
182 }
183
RHSFunctionV(TS ts,PetscReal t,Vec X,Vec Vres,PetscCtx ctx)184 static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, PetscCtx ctx)
185 {
186 DM sw;
187 const PetscScalar *x;
188 const PetscReal *coords, *vel;
189 PetscScalar *vres;
190 PetscInt Np, p, dim, d;
191
192 PetscFunctionBeginUser;
193 PetscCall(TSGetDM(ts, &sw));
194 PetscCall(DMGetDimension(sw, &dim));
195 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
196 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
197 PetscCall(VecGetLocalSize(Vres, &Np));
198 PetscCall(VecGetArrayRead(X, &x));
199 PetscCall(VecGetArray(Vres, &vres));
200 PetscCheck(dim == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension must be 2");
201 Np /= dim;
202 for (p = 0; p < Np; ++p) {
203 const PetscReal x0 = coords[p * dim + 0];
204 const PetscReal vy0 = vel[p * dim + 1];
205 const PetscReal omega = vy0 / x0;
206
207 for (d = 0; d < dim; ++d) vres[p * dim + d] = -PetscSqr(omega) * x[p * dim + d];
208 }
209 PetscCall(VecRestoreArrayRead(X, &x));
210 PetscCall(VecRestoreArray(Vres, &vres));
211 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
212 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
213 PetscFunctionReturn(PETSC_SUCCESS);
214 }
215
CreateSolution(TS ts)216 static PetscErrorCode CreateSolution(TS ts)
217 {
218 DM sw;
219 Vec u;
220 PetscInt dim, Np;
221
222 PetscFunctionBegin;
223 PetscCall(TSGetDM(ts, &sw));
224 PetscCall(DMGetDimension(sw, &dim));
225 PetscCall(DMSwarmGetLocalSize(sw, &Np));
226 PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
227 PetscCall(VecSetBlockSize(u, dim));
228 PetscCall(VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE));
229 PetscCall(VecSetUp(u));
230 PetscCall(TSSetSolution(ts, u));
231 PetscCall(VecDestroy(&u));
232 PetscFunctionReturn(PETSC_SUCCESS);
233 }
234
SetProblem(TS ts)235 static PetscErrorCode SetProblem(TS ts)
236 {
237 AppCtx *user;
238 DM sw;
239
240 PetscFunctionBegin;
241 PetscCall(TSGetDM(ts, &sw));
242 PetscCall(DMGetApplicationContext(sw, &user));
243 // Define unified system for (X, V)
244 {
245 Mat J;
246 PetscInt dim, Np;
247
248 PetscCall(DMGetDimension(sw, &dim));
249 PetscCall(DMSwarmGetLocalSize(sw, &Np));
250 PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
251 PetscCall(MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE));
252 PetscCall(MatSetBlockSize(J, 2 * dim));
253 PetscCall(MatSetFromOptions(J));
254 PetscCall(MatSetUp(J));
255 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, user));
256 PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, user));
257 PetscCall(MatDestroy(&J));
258 }
259 // Define split system for X and V
260 {
261 Vec u;
262 IS isx, isv, istmp;
263 const PetscInt *idx;
264 PetscInt dim, Np, rstart;
265
266 PetscCall(TSGetSolution(ts, &u));
267 PetscCall(DMGetDimension(sw, &dim));
268 PetscCall(DMSwarmGetLocalSize(sw, &Np));
269 PetscCall(VecGetOwnershipRange(u, &rstart, NULL));
270 PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp));
271 PetscCall(ISGetIndices(istmp, &idx));
272 PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx));
273 PetscCall(ISRestoreIndices(istmp, &idx));
274 PetscCall(ISDestroy(&istmp));
275 PetscCall(ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp));
276 PetscCall(ISGetIndices(istmp, &idx));
277 PetscCall(ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv));
278 PetscCall(ISRestoreIndices(istmp, &idx));
279 PetscCall(ISDestroy(&istmp));
280 PetscCall(TSRHSSplitSetIS(ts, "position", isx));
281 PetscCall(TSRHSSplitSetIS(ts, "momentum", isv));
282 PetscCall(ISDestroy(&isx));
283 PetscCall(ISDestroy(&isv));
284 PetscCall(TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, user));
285 PetscCall(TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, user));
286 }
287 // Define symplectic formulation U_t = S . G, where G = grad F
288 {
289 //PetscCall(TSDiscGradSetFormulation(ts, RHSJacobianS, RHSObjectiveF, RHSFunctionG, user));
290 }
291 PetscFunctionReturn(PETSC_SUCCESS);
292 }
293
DMSwarmTSRedistribute(TS ts)294 static PetscErrorCode DMSwarmTSRedistribute(TS ts)
295 {
296 DM sw;
297 Vec u;
298 PetscReal t, maxt, dt;
299 PetscInt n, maxn;
300
301 PetscFunctionBegin;
302 PetscCall(TSGetDM(ts, &sw));
303 PetscCall(TSGetTime(ts, &t));
304 PetscCall(TSGetMaxTime(ts, &maxt));
305 PetscCall(TSGetTimeStep(ts, &dt));
306 PetscCall(TSGetStepNumber(ts, &n));
307 PetscCall(TSGetMaxSteps(ts, &maxn));
308
309 PetscCall(TSReset(ts));
310 PetscCall(TSSetDM(ts, sw));
311 /* TODO Check whether TS was set from options */
312 PetscCall(TSSetFromOptions(ts));
313 PetscCall(TSSetTime(ts, t));
314 PetscCall(TSSetMaxTime(ts, maxt));
315 PetscCall(TSSetTimeStep(ts, dt));
316 PetscCall(TSSetStepNumber(ts, n));
317 PetscCall(TSSetMaxSteps(ts, maxn));
318
319 PetscCall(CreateSolution(ts));
320 PetscCall(SetProblem(ts));
321 PetscCall(TSGetSolution(ts, &u));
322 PetscFunctionReturn(PETSC_SUCCESS);
323 }
324
circleSingleX(PetscInt dim,PetscReal time,const PetscReal unused[],PetscInt p,PetscScalar x[],PetscCtx ctx)325 PetscErrorCode circleSingleX(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar x[], PetscCtx ctx)
326 {
327 x[0] = p + 1.;
328 x[1] = 0.;
329 return PETSC_SUCCESS;
330 }
331
circleSingleV(PetscInt dim,PetscReal time,const PetscReal unused[],PetscInt p,PetscScalar v[],PetscCtx ctx)332 PetscErrorCode circleSingleV(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar v[], PetscCtx ctx)
333 {
334 v[0] = 0.;
335 v[1] = PetscSqrtReal(1000. / (p + 1.));
336 return PETSC_SUCCESS;
337 }
338
339 /* Put 5 particles into each circle */
circleMultipleX(PetscInt dim,PetscReal time,const PetscReal unused[],PetscInt p,PetscScalar x[],PetscCtx ctx)340 PetscErrorCode circleMultipleX(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar x[], PetscCtx ctx)
341 {
342 const PetscInt n = 5;
343 const PetscReal r0 = (p / n) + 1.;
344 const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
345
346 x[0] = r0 * PetscCosReal(th0);
347 x[1] = r0 * PetscSinReal(th0);
348 return PETSC_SUCCESS;
349 }
350
351 /* Put 5 particles into each circle */
circleMultipleV(PetscInt dim,PetscReal time,const PetscReal unused[],PetscInt p,PetscScalar v[],PetscCtx ctx)352 PetscErrorCode circleMultipleV(PetscInt dim, PetscReal time, const PetscReal unused[], PetscInt p, PetscScalar v[], PetscCtx ctx)
353 {
354 const PetscInt n = 5;
355 const PetscReal r0 = (p / n) + 1.;
356 const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
357 const PetscReal omega = PetscSqrtReal(1000. / r0) / r0;
358
359 v[0] = -r0 * omega * PetscSinReal(th0);
360 v[1] = r0 * omega * PetscCosReal(th0);
361 return PETSC_SUCCESS;
362 }
363
364 /*
365 InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
366
367 Input Parameters:
368 + ts - The TS
369 - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem
370
371 Output Parameter:
372 . u - The initialized solution vector
373
374 Level: advanced
375
376 .seealso: InitializeSolve()
377 */
InitializeSolveAndSwarm(TS ts,PetscBool useInitial)378 static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
379 {
380 DM sw;
381 Vec u, gc, gv, gc0, gv0;
382 IS isx, isv;
383 AppCtx *user;
384
385 PetscFunctionBeginUser;
386 PetscCall(TSGetDM(ts, &sw));
387 PetscCall(DMGetApplicationContext(sw, &user));
388 if (useInitial) {
389 PetscReal v0[1] = {1.};
390
391 PetscCall(DMSwarmInitializeCoordinates(sw));
392 PetscCall(DMSwarmInitializeVelocitiesFromOptions(sw, v0));
393 PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
394 PetscCall(DMSwarmTSRedistribute(ts));
395 }
396 PetscCall(TSGetSolution(ts, &u));
397 PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
398 PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
399 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
400 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0));
401 if (useInitial) PetscCall(VecCopy(gc, gc0));
402 PetscCall(VecISCopy(u, isx, SCATTER_FORWARD, gc));
403 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
404 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0));
405 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
406 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "initVelocity", &gv0));
407 if (useInitial) PetscCall(VecCopy(gv, gv0));
408 PetscCall(VecISCopy(u, isv, SCATTER_FORWARD, gv));
409 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
410 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "initVelocity", &gv0));
411 PetscFunctionReturn(PETSC_SUCCESS);
412 }
413
InitializeSolve(TS ts,Vec u)414 static PetscErrorCode InitializeSolve(TS ts, Vec u)
415 {
416 PetscFunctionBegin;
417 PetscCall(TSSetSolution(ts, u));
418 PetscCall(InitializeSolveAndSwarm(ts, PETSC_TRUE));
419 PetscFunctionReturn(PETSC_SUCCESS);
420 }
421
ComputeError(TS ts,Vec U,Vec E)422 static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
423 {
424 MPI_Comm comm;
425 DM sw;
426 AppCtx *user;
427 const PetscScalar *u;
428 const PetscReal *coords, *vel;
429 PetscScalar *e;
430 PetscReal t;
431 PetscInt dim, Np, p;
432
433 PetscFunctionBeginUser;
434 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
435 PetscCall(TSGetDM(ts, &sw));
436 PetscCall(DMGetApplicationContext(sw, &user));
437 PetscCall(DMGetDimension(sw, &dim));
438 PetscCall(TSGetSolveTime(ts, &t));
439 PetscCall(VecGetArray(E, &e));
440 PetscCall(VecGetArrayRead(U, &u));
441 PetscCall(VecGetLocalSize(U, &Np));
442 PetscCall(DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
443 PetscCall(DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel));
444 Np /= 2 * dim;
445 for (p = 0; p < Np; ++p) {
446 /* TODO generalize initial conditions and project into plane instead of assuming x-y */
447 const PetscReal r0 = DMPlex_NormD_Internal(dim, &coords[p * dim]);
448 const PetscReal th0 = PetscAtan2Real(coords[p * dim + 1], coords[p * dim + 0]);
449 const PetscReal v0 = DMPlex_NormD_Internal(dim, &vel[p * dim]);
450 const PetscReal omega = v0 / r0;
451 const PetscReal ct = PetscCosReal(omega * t + th0);
452 const PetscReal st = PetscSinReal(omega * t + th0);
453 const PetscScalar *x = &u[(p * 2 + 0) * dim];
454 const PetscScalar *v = &u[(p * 2 + 1) * dim];
455 const PetscReal xe[3] = {r0 * ct, r0 * st, 0.0};
456 const PetscReal ve[3] = {-v0 * st, v0 * ct, 0.0};
457 PetscInt d;
458
459 for (d = 0; d < dim; ++d) {
460 e[(p * 2 + 0) * dim + d] = x[d] - xe[d];
461 e[(p * 2 + 1) * dim + d] = v[d] - ve[d];
462 }
463 if (user->error) {
464 const PetscReal en = 0.5 * DMPlex_DotRealD_Internal(dim, v, v);
465 const PetscReal exen = 0.5 * PetscSqr(v0);
466 PetscCall(PetscPrintf(comm, "t %.4g: p%" PetscInt_FMT " error [%.2f %.2f] sol [(%.6lf %.6lf) (%.6lf %.6lf)] exact [(%.6lf %.6lf) (%.6lf %.6lf)] energy/exact energy %g / %g (%.10lf%%)\n", (double)t, p, (double)DMPlex_NormD_Internal(dim, &e[(p * 2 + 0) * dim]), (double)DMPlex_NormD_Internal(dim, &e[(p * 2 + 1) * dim]), (double)x[0], (double)x[1], (double)v[0], (double)v[1], (double)xe[0], (double)xe[1], (double)ve[0], (double)ve[1], (double)en, (double)exen, (double)(PetscAbsReal(exen - en) * 100. / exen)));
467 }
468 }
469 PetscCall(DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords));
470 PetscCall(DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel));
471 PetscCall(VecRestoreArrayRead(U, &u));
472 PetscCall(VecRestoreArray(E, &e));
473 PetscFunctionReturn(PETSC_SUCCESS);
474 }
475
EnergyMonitor(TS ts,PetscInt step,PetscReal t,Vec U,PetscCtx ctx)476 static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, PetscCtx ctx)
477 {
478 const PetscInt ostep = ((AppCtx *)ctx)->ostep;
479 DM sw;
480 const PetscScalar *u;
481 PetscInt dim, Np, p;
482
483 PetscFunctionBeginUser;
484 if (step % ostep == 0) {
485 PetscCall(TSGetDM(ts, &sw));
486 PetscCall(DMGetDimension(sw, &dim));
487 PetscCall(VecGetArrayRead(U, &u));
488 PetscCall(VecGetLocalSize(U, &Np));
489 Np /= 2 * dim;
490 if (!step) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Time Step Part Energy\n"));
491 for (p = 0; p < Np; ++p) {
492 const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);
493
494 PetscCall(PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " %4" PetscInt_FMT " %10.4lf\n", (double)t, step, p, (double)(0.5 * v2)));
495 }
496 PetscCall(PetscSynchronizedFlush(PetscObjectComm((PetscObject)ts), NULL));
497 PetscCall(VecRestoreArrayRead(U, &u));
498 }
499 PetscFunctionReturn(PETSC_SUCCESS);
500 }
501
MigrateParticles(TS ts)502 static PetscErrorCode MigrateParticles(TS ts)
503 {
504 DM sw;
505
506 PetscFunctionBeginUser;
507 PetscCall(TSGetDM(ts, &sw));
508 PetscCall(DMViewFromOptions(sw, NULL, "-migrate_view_pre"));
509 {
510 Vec u, gc, gv;
511 IS isx, isv;
512
513 PetscCall(TSGetSolution(ts, &u));
514 PetscCall(TSRHSSplitGetIS(ts, "position", &isx));
515 PetscCall(TSRHSSplitGetIS(ts, "momentum", &isv));
516 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
517 PetscCall(VecISCopy(u, isx, SCATTER_REVERSE, gc));
518 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc));
519 PetscCall(DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv));
520 PetscCall(VecISCopy(u, isv, SCATTER_REVERSE, gv));
521 PetscCall(DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv));
522 }
523 PetscCall(DMSwarmMigrate(sw, PETSC_TRUE));
524 PetscCall(DMSwarmTSRedistribute(ts));
525 PetscCall(InitializeSolveAndSwarm(ts, PETSC_FALSE));
526 PetscFunctionReturn(PETSC_SUCCESS);
527 }
528
main(int argc,char ** argv)529 int main(int argc, char **argv)
530 {
531 DM dm, sw;
532 TS ts;
533 Vec u;
534 AppCtx user;
535
536 PetscFunctionBeginUser;
537 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
538 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
539 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
540 PetscCall(CreateSwarm(dm, &user, &sw));
541 PetscCall(DMSetApplicationContext(sw, &user));
542
543 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
544 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
545 PetscCall(TSSetDM(ts, sw));
546 PetscCall(TSSetMaxTime(ts, 0.1));
547 PetscCall(TSSetTimeStep(ts, 0.00001));
548 PetscCall(TSSetMaxSteps(ts, 100));
549 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
550 PetscCall(TSMonitorSet(ts, EnergyMonitor, &user, NULL));
551 PetscCall(TSSetFromOptions(ts));
552 PetscCall(TSSetComputeInitialCondition(ts, InitializeSolve));
553 PetscCall(TSSetComputeExactError(ts, ComputeError));
554 PetscCall(TSSetPostStep(ts, MigrateParticles));
555
556 PetscCall(CreateSolution(ts));
557 PetscCall(TSGetSolution(ts, &u));
558 PetscCall(TSComputeInitialCondition(ts, u));
559 PetscCall(TSSolve(ts, NULL));
560
561 PetscCall(TSDestroy(&ts));
562 PetscCall(DMDestroy(&sw));
563 PetscCall(DMDestroy(&dm));
564 PetscCall(PetscFinalize());
565 return 0;
566 }
567
568 /*TEST
569
570 build:
571 requires: double !complex
572
573 testset:
574 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
575 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
576 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
577 -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
578 -dm_view -output_step 50 -error -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10
579 test:
580 suffix: bsi_2d_1
581 args: -ts_basicsymplectic_type 1
582 test:
583 suffix: bsi_2d_2
584 args: -ts_basicsymplectic_type 2
585 test:
586 suffix: bsi_2d_3
587 args: -ts_basicsymplectic_type 3
588 test:
589 suffix: bsi_2d_4
590 args: -ts_basicsymplectic_type 4 -ts_time_step 0.0001
591
592 testset:
593 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
594 args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
595 -ts_type theta -ts_theta_theta 0.5 -ts_convergence_estimate -convest_num_refine 2 \
596 -mat_type baij -ksp_error_if_not_converged -pc_type lu \
597 -dm_view -output_step 50 -error -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10
598 test:
599 suffix: im_2d_0
600 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5
601
602 testset:
603 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
604 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 10,10 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 -petscpartitioner_type simple \
605 -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
606 -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
607 -dm_view -output_step 50 -ts_time_step 0.01 -ts_max_time 10.0 -ts_max_steps 10
608 test:
609 suffix: bsi_2d_mesh_1
610 args: -ts_basicsymplectic_type 1 -error
611 test:
612 suffix: bsi_2d_mesh_1_par_2
613 nsize: 2
614 args: -ts_basicsymplectic_type 1
615 test:
616 suffix: bsi_2d_mesh_1_par_3
617 nsize: 3
618 args: -ts_basicsymplectic_type 1
619 test:
620 suffix: bsi_2d_mesh_1_par_4
621 nsize: 4
622 args: -ts_basicsymplectic_type 1 -dm_swarm_num_particles 0,0,2,0
623
624 testset:
625 requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
626 args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 10,10 -dm_plex_box_lower -5,-5 -dm_plex_box_upper 5,5 \
627 -dm_swarm_num_particles 10 -dm_swarm_coordinate_function circleMultipleX -dm_swarm_velocity_function circleMultipleV \
628 -ts_convergence_estimate -convest_num_refine 2 \
629 -dm_view -output_step 50 -error
630 test:
631 suffix: bsi_2d_multiple_1
632 args: -ts_type basicsymplectic -ts_basicsymplectic_type 1
633 test:
634 suffix: bsi_2d_multiple_2
635 args: -ts_type basicsymplectic -ts_basicsymplectic_type 2
636 test:
637 suffix: bsi_2d_multiple_3
638 args: -ts_type basicsymplectic -ts_basicsymplectic_type 3 -ts_time_step 0.001
639 test:
640 suffix: im_2d_multiple_0
641 args: -ts_type theta -ts_theta_theta 0.5 \
642 -mat_type baij -ksp_error_if_not_converged -pc_type lu
643
644 TEST*/
645