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*/