Actual source code: ex5.c
1: static char help[] = "Vlasov example of central orbits\n";
3: /*
4: To visualize the orbit, we can used
6: -ts_monitor_sp_swarm -ts_monitor_sp_swarm_retain -1 -ts_monitor_sp_swarm_phase 0 -draw_size 500,500
8: and we probably want it to run fast and not check convergence
10: -convest_num_refine 0 -ts_dt 0.01 -ts_max_steps 100 -ts_max_time 100 -output_step 10
11: */
13: #include <petscts.h>
14: #include <petscdmplex.h>
15: #include <petscdmswarm.h>
16: #include <petsc/private/dmpleximpl.h>
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 *);
23: typedef struct {
24: PetscBool error; /* Flag for printing the error */
25: PetscInt ostep; /* print the energy at each ostep time steps */
26: } AppCtx;
28: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
29: {
30: PetscFunctionBeginUser;
31: options->error = PETSC_FALSE;
32: options->ostep = 100;
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: }
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: }
51: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
52: {
53: PetscReal v0[1] = {1.};
54: PetscInt dim;
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;
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: }
93: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
94: {
95: DM sw;
96: const PetscReal *coords, *vel;
97: const PetscScalar *u;
98: PetscScalar *g;
99: PetscInt dim, d, Np, p;
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;
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: }
127: /* J_{ij} = dF_i/dx_j
128: J_p = ( 0 1)
129: (-w^2 0)
130: */
131: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat P, void *ctx)
132: {
133: DM sw;
134: const PetscReal *coords, *vel;
135: PetscInt dim, d, Np, p, rStart;
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.};
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: }
163: static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
164: {
165: DM sw;
166: const PetscScalar *v;
167: PetscScalar *xres;
168: PetscInt Np, p, dim, d;
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: }
184: static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
185: {
186: DM sw;
187: const PetscScalar *x;
188: const PetscReal *coords, *vel;
189: PetscScalar *vres;
190: PetscInt Np, p, dim, d;
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;
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: }
216: static PetscErrorCode CreateSolution(TS ts)
217: {
218: DM sw;
219: Vec u;
220: PetscInt dim, Np;
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: }
235: static PetscErrorCode SetProblem(TS ts)
236: {
237: AppCtx *user;
238: DM sw;
240: PetscFunctionBegin;
241: PetscCall(TSGetDM(ts, &sw));
242: PetscCall(DMGetApplicationContext(sw, (void **)&user));
243: // Define unified system for (X, V)
244: {
245: Mat J;
246: PetscInt dim, Np;
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;
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: }
294: static PetscErrorCode DMSwarmTSRedistribute(TS ts)
295: {
296: DM sw;
297: Vec u;
298: PetscReal t, maxt, dt;
299: PetscInt n, maxn;
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));
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));
319: PetscCall(CreateSolution(ts));
320: PetscCall(SetProblem(ts));
321: PetscCall(TSGetSolution(ts, &u));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: PetscErrorCode circleSingleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
326: {
327: x[0] = p + 1.;
328: x[1] = 0.;
329: return PETSC_SUCCESS;
330: }
332: PetscErrorCode circleSingleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx)
333: {
334: v[0] = 0.;
335: v[1] = PetscSqrtReal(1000. / (p + 1.));
336: return PETSC_SUCCESS;
337: }
339: /* Put 5 particles into each circle */
340: PetscErrorCode circleMultipleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
341: {
342: const PetscInt n = 5;
343: const PetscReal r0 = (p / n) + 1.;
344: const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;
346: x[0] = r0 * PetscCosReal(th0);
347: x[1] = r0 * PetscSinReal(th0);
348: return PETSC_SUCCESS;
349: }
351: /* Put 5 particles into each circle */
352: PetscErrorCode circleMultipleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *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;
359: v[0] = -r0 * omega * PetscSinReal(th0);
360: v[1] = r0 * omega * PetscCosReal(th0);
361: return PETSC_SUCCESS;
362: }
364: /*
365: InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.
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
371: Output Parameter:
372: . u - The initialized solution vector
374: Level: advanced
376: .seealso: InitializeSolve()
377: */
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;
385: PetscFunctionBeginUser;
386: PetscCall(TSGetDM(ts, &sw));
387: PetscCall(DMGetApplicationContext(sw, &user));
388: if (useInitial) {
389: PetscReal v0[1] = {1.};
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: }
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: }
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;
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;
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: }
476: static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
477: {
478: const PetscInt ostep = ((AppCtx *)ctx)->ostep;
479: DM sw;
480: const PetscScalar *u;
481: PetscInt dim, Np, p;
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]);
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: }
502: static PetscErrorCode MigrateParticles(TS ts)
503: {
504: DM sw;
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;
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: }
529: int main(int argc, char **argv)
530: {
531: DM dm, sw;
532: TS ts;
533: Vec u;
534: AppCtx user;
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));
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));
556: PetscCall(CreateSolution(ts));
557: PetscCall(TSGetSolution(ts, &u));
558: PetscCall(TSComputeInitialCondition(ts, u));
559: PetscCall(TSSolve(ts, NULL));
561: PetscCall(TSDestroy(&ts));
562: PetscCall(DMDestroy(&sw));
563: PetscCall(DMDestroy(&dm));
564: PetscCall(PetscFinalize());
565: return 0;
566: }
568: /*TEST
570: build:
571: requires: double !complex
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_dt 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_dt 0.0001
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_dt 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
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_dt 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
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_dt 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
644: TEST*/