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: {
 31:   options->error = PETSC_FALSE;
 32:   options->ostep = 100;

 34:   PetscOptionsBegin(comm, "", "Central Orbit Options", "DMSWARM");
 35:   PetscOptionsBool("-error", "Flag to print the error", "ex5.c", options->error, &options->error, NULL);
 36:   PetscOptionsInt("-output_step", "Number of time steps between output", "ex5.c", options->ostep, &options->ostep, PETSC_NULL);
 37:   PetscOptionsEnd();
 38:   return 0;
 39: }

 41: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 42: {
 44:   DMCreate(comm, dm);
 45:   DMSetType(*dm, DMPLEX);
 46:   DMSetFromOptions(*dm);
 47:   DMViewFromOptions(*dm, NULL, "-dm_view");
 48:   return 0;
 49: }

 51: static PetscErrorCode CreateSwarm(DM dm, AppCtx *user, DM *sw)
 52: {
 53:   PetscReal v0[1] = {1.};
 54:   PetscInt  dim;

 57:   DMGetDimension(dm, &dim);
 58:   DMCreate(PetscObjectComm((PetscObject)dm), sw);
 59:   DMSetType(*sw, DMSWARM);
 60:   DMSetDimension(*sw, dim);
 61:   DMSwarmSetType(*sw, DMSWARM_PIC);
 62:   DMSwarmSetCellDM(*sw, dm);
 63:   DMSwarmRegisterPetscDatatypeField(*sw, "w_q", 1, PETSC_SCALAR);
 64:   DMSwarmRegisterPetscDatatypeField(*sw, "velocity", dim, PETSC_REAL);
 65:   DMSwarmRegisterPetscDatatypeField(*sw, "species", 1, PETSC_INT);
 66:   DMSwarmRegisterPetscDatatypeField(*sw, "initCoordinates", dim, PETSC_REAL);
 67:   DMSwarmRegisterPetscDatatypeField(*sw, "initVelocity", dim, PETSC_REAL);
 68:   DMSwarmFinalizeFieldRegister(*sw);
 69:   DMSwarmComputeLocalSizeFromOptions(*sw);
 70:   DMSwarmInitializeCoordinates(*sw);
 71:   DMSwarmInitializeVelocitiesFromOptions(*sw, v0);
 72:   DMSetFromOptions(*sw);
 73:   DMSetApplicationContext(*sw, user);
 74:   PetscObjectSetName((PetscObject)*sw, "Particles");
 75:   DMViewFromOptions(*sw, NULL, "-sw_view");
 76:   {
 77:     Vec gc, gc0, gv, gv0;

 79:     DMSwarmCreateGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc);
 80:     DMSwarmCreateGlobalVectorFromField(*sw, "initCoordinates", &gc0);
 81:     VecCopy(gc, gc0);
 82:     DMSwarmDestroyGlobalVectorFromField(*sw, DMSwarmPICField_coor, &gc);
 83:     DMSwarmDestroyGlobalVectorFromField(*sw, "initCoordinates", &gc0);
 84:     DMSwarmCreateGlobalVectorFromField(*sw, "velocity", &gv);
 85:     DMSwarmCreateGlobalVectorFromField(*sw, "initVelocity", &gv0);
 86:     VecCopy(gv, gv0);
 87:     DMSwarmDestroyGlobalVectorFromField(*sw, "velocity", &gv);
 88:     DMSwarmDestroyGlobalVectorFromField(*sw, "initVelocity", &gv0);
 89:   }
 90:   return 0;
 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;

102:   TSGetDM(ts, &sw);
103:   DMGetDimension(sw, &dim);
104:   DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
105:   DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel);
106:   VecGetLocalSize(U, &Np);
107:   VecGetArrayRead(U, &u);
108:   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:   DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
121:   DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel);
122:   VecRestoreArrayRead(U, &u);
123:   VecRestoreArray(G, &g);
124:   return 0;
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;

138:   TSGetDM(ts, &sw);
139:   DMGetDimension(sw, &dim);
140:   VecGetLocalSize(U, &Np);
141:   MatGetOwnershipRange(J, &rStart, NULL);
142:   DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
143:   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:       MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES);
154:     }
155:   }
156:   DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
157:   DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel);
158:   MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
159:   MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
160:   return 0;
161: }

163: static PetscErrorCode RHSFunctionX(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
164: {
165:   const PetscScalar *v;
166:   PetscScalar       *xres;
167:   PetscInt           Np, p;

170:   VecGetLocalSize(Xres, &Np);
171:   VecGetArrayRead(V, &v);
172:   VecGetArray(Xres, &xres);
173:   for (p = 0; p < Np; ++p) xres[p] = v[p];
174:   VecRestoreArrayRead(V, &v);
175:   VecRestoreArray(Xres, &xres);
176:   return 0;
177: }

179: static PetscErrorCode RHSFunctionV(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
180: {
181:   DM                 sw;
182:   const PetscScalar *x;
183:   const PetscReal   *coords, *vel;
184:   PetscScalar       *vres;
185:   PetscInt           Np, p, dim, d;

188:   TSGetDM(ts, &sw);
189:   DMGetDimension(sw, &dim);
190:   DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
191:   DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel);
192:   VecGetLocalSize(Vres, &Np);
193:   VecGetArrayRead(X, &x);
194:   VecGetArray(Vres, &vres);
196:   Np /= dim;
197:   for (p = 0; p < Np; ++p) {
198:     const PetscReal x0    = coords[p * dim + 0];
199:     const PetscReal vy0   = vel[p * dim + 1];
200:     const PetscReal omega = vy0 / x0;

202:     for (d = 0; d < dim; ++d) vres[p * dim + d] = -PetscSqr(omega) * x[p * dim + d];
203:   }
204:   VecRestoreArrayRead(X, &x);
205:   VecRestoreArray(Vres, &vres);
206:   DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
207:   DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel);
208:   return 0;
209: }

211: static PetscErrorCode CreateSolution(TS ts)
212: {
213:   DM       sw;
214:   Vec      u;
215:   PetscInt dim, Np;

217:   TSGetDM(ts, &sw);
218:   DMGetDimension(sw, &dim);
219:   DMSwarmGetLocalSize(sw, &Np);
220:   VecCreate(PETSC_COMM_WORLD, &u);
221:   VecSetBlockSize(u, dim);
222:   VecSetSizes(u, 2 * Np * dim, PETSC_DECIDE);
223:   VecSetUp(u);
224:   TSSetSolution(ts, u);
225:   VecDestroy(&u);
226:   return 0;
227: }

229: static PetscErrorCode SetProblem(TS ts)
230: {
231:   AppCtx *user;
232:   DM      sw;

234:   TSGetDM(ts, &sw);
235:   DMGetApplicationContext(sw, (void **)&user);
236:   // Define unified system for (X, V)
237:   {
238:     Mat      J;
239:     PetscInt dim, Np;

241:     DMGetDimension(sw, &dim);
242:     DMSwarmGetLocalSize(sw, &Np);
243:     MatCreate(PETSC_COMM_WORLD, &J);
244:     MatSetSizes(J, 2 * Np * dim, 2 * Np * dim, PETSC_DECIDE, PETSC_DECIDE);
245:     MatSetBlockSize(J, 2 * dim);
246:     MatSetFromOptions(J);
247:     MatSetUp(J);
248:     TSSetRHSFunction(ts, NULL, RHSFunction, &user);
249:     TSSetRHSJacobian(ts, J, J, RHSJacobian, &user);
250:     MatDestroy(&J);
251:   }
252:   // Define split system for X and V
253:   {
254:     Vec             u;
255:     IS              isx, isv, istmp;
256:     const PetscInt *idx;
257:     PetscInt        dim, Np, rstart;

259:     TSGetSolution(ts, &u);
260:     DMGetDimension(sw, &dim);
261:     DMSwarmGetLocalSize(sw, &Np);
262:     VecGetOwnershipRange(u, &rstart, NULL);
263:     ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 0, 2, &istmp);
264:     ISGetIndices(istmp, &idx);
265:     ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isx);
266:     ISRestoreIndices(istmp, &idx);
267:     ISDestroy(&istmp);
268:     ISCreateStride(PETSC_COMM_WORLD, Np, (rstart / dim) + 1, 2, &istmp);
269:     ISGetIndices(istmp, &idx);
270:     ISCreateBlock(PETSC_COMM_WORLD, dim, Np, idx, PETSC_COPY_VALUES, &isv);
271:     ISRestoreIndices(istmp, &idx);
272:     ISDestroy(&istmp);
273:     TSRHSSplitSetIS(ts, "position", isx);
274:     TSRHSSplitSetIS(ts, "momentum", isv);
275:     ISDestroy(&isx);
276:     ISDestroy(&isv);
277:     TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunctionX, &user);
278:     TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunctionV, &user);
279:   }
280:   // Define symplectic formulation U_t = S . G, where G = grad F
281:   {
282:     //TSDiscGradSetFormulation(ts, RHSJacobianS, RHSObjectiveF, RHSFunctionG, &user);
283:   }
284:   return 0;
285: }

287: static PetscErrorCode DMSwarmTSRedistribute(TS ts)
288: {
289:   DM        sw;
290:   Vec       u;
291:   PetscReal t, maxt, dt;
292:   PetscInt  n, maxn;

294:   TSGetDM(ts, &sw);
295:   TSGetTime(ts, &t);
296:   TSGetMaxTime(ts, &maxt);
297:   TSGetTimeStep(ts, &dt);
298:   TSGetStepNumber(ts, &n);
299:   TSGetMaxSteps(ts, &maxn);

301:   TSReset(ts);
302:   TSSetDM(ts, sw);
303:   /* TODO Check whether TS was set from options */
304:   TSSetFromOptions(ts);
305:   TSSetTime(ts, t);
306:   TSSetMaxTime(ts, maxt);
307:   TSSetTimeStep(ts, dt);
308:   TSSetStepNumber(ts, n);
309:   TSSetMaxSteps(ts, maxn);

311:   CreateSolution(ts);
312:   SetProblem(ts);
313:   TSGetSolution(ts, &u);
314:   return 0;
315: }

317: PetscErrorCode circleSingleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
318: {
319:   x[0] = p + 1.;
320:   x[1] = 0.;
321:   return 0;
322: }

324: PetscErrorCode circleSingleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx)
325: {
326:   v[0] = 0.;
327:   v[1] = PetscSqrtReal(1000. / (p + 1.));
328:   return 0;
329: }

331: /* Put 5 particles into each circle */
332: PetscErrorCode circleMultipleX(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar x[], void *ctx)
333: {
334:   const PetscInt  n   = 5;
335:   const PetscReal r0  = (p / n) + 1.;
336:   const PetscReal th0 = (2. * PETSC_PI * (p % n)) / n;

338:   x[0] = r0 * PetscCosReal(th0);
339:   x[1] = r0 * PetscSinReal(th0);
340:   return 0;
341: }

343: /* Put 5 particles into each circle */
344: PetscErrorCode circleMultipleV(PetscInt dim, PetscReal time, const PetscReal dummy[], PetscInt p, PetscScalar v[], void *ctx)
345: {
346:   const PetscInt  n     = 5;
347:   const PetscReal r0    = (p / n) + 1.;
348:   const PetscReal th0   = (2. * PETSC_PI * (p % n)) / n;
349:   const PetscReal omega = PetscSqrtReal(1000. / r0) / r0;

351:   v[0] = -r0 * omega * PetscSinReal(th0);
352:   v[1] = r0 * omega * PetscCosReal(th0);
353:   return 0;
354: }

356: /*
357:   InitializeSolveAndSwarm - Set the solution values to the swarm coordinates and velocities, and also possibly set the initial values.

359:   Input Parameters:
360: + ts         - The TS
361: - useInitial - Flag to also set the initial conditions to the current coordinates and velocities and setup the problem

363:   Output Parameters:
364: . u - The initialized solution vector

366:   Level: advanced

368: .seealso: InitializeSolve()
369: */
370: static PetscErrorCode InitializeSolveAndSwarm(TS ts, PetscBool useInitial)
371: {
372:   DM      sw;
373:   Vec     u, gc, gv, gc0, gv0;
374:   IS      isx, isv;
375:   AppCtx *user;

378:   TSGetDM(ts, &sw);
379:   DMGetApplicationContext(sw, &user);
380:   if (useInitial) {
381:     PetscReal v0[1] = {1.};

383:     DMSwarmInitializeCoordinates(sw);
384:     DMSwarmInitializeVelocitiesFromOptions(sw, v0);
385:     DMSwarmMigrate(sw, PETSC_TRUE);
386:     DMSwarmTSRedistribute(ts);
387:   }
388:   TSGetSolution(ts, &u);
389:   TSRHSSplitGetIS(ts, "position", &isx);
390:   TSRHSSplitGetIS(ts, "momentum", &isv);
391:   DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc);
392:   DMSwarmCreateGlobalVectorFromField(sw, "initCoordinates", &gc0);
393:   if (useInitial) VecCopy(gc, gc0);
394:   VecISCopy(u, isx, SCATTER_FORWARD, gc);
395:   DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc);
396:   DMSwarmDestroyGlobalVectorFromField(sw, "initCoordinates", &gc0);
397:   DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv);
398:   DMSwarmCreateGlobalVectorFromField(sw, "initVelocity", &gv0);
399:   if (useInitial) VecCopy(gv, gv0);
400:   VecISCopy(u, isv, SCATTER_FORWARD, gv);
401:   DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv);
402:   DMSwarmDestroyGlobalVectorFromField(sw, "initVelocity", &gv0);
403:   return 0;
404: }

406: static PetscErrorCode InitializeSolve(TS ts, Vec u)
407: {
408:   TSSetSolution(ts, u);
409:   InitializeSolveAndSwarm(ts, PETSC_TRUE);
410:   return 0;
411: }

413: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
414: {
415:   MPI_Comm           comm;
416:   DM                 sw;
417:   AppCtx            *user;
418:   const PetscScalar *u;
419:   const PetscReal   *coords, *vel;
420:   PetscScalar       *e;
421:   PetscReal          t;
422:   PetscInt           dim, Np, p;

425:   PetscObjectGetComm((PetscObject)ts, &comm);
426:   TSGetDM(ts, &sw);
427:   DMGetApplicationContext(sw, &user);
428:   DMGetDimension(sw, &dim);
429:   TSGetSolveTime(ts, &t);
430:   VecGetArray(E, &e);
431:   VecGetArrayRead(U, &u);
432:   VecGetLocalSize(U, &Np);
433:   DMSwarmGetField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
434:   DMSwarmGetField(sw, "initVelocity", NULL, NULL, (void **)&vel);
435:   Np /= 2 * dim;
436:   for (p = 0; p < Np; ++p) {
437:     /* TODO generalize initial conditions and project into plane instead of assuming x-y */
438:     const PetscReal    r0    = DMPlex_NormD_Internal(dim, &coords[p * dim]);
439:     const PetscReal    th0   = PetscAtan2Real(coords[p * dim + 1], coords[p * dim + 0]);
440:     const PetscReal    v0    = DMPlex_NormD_Internal(dim, &vel[p * dim]);
441:     const PetscReal    omega = v0 / r0;
442:     const PetscReal    ct    = PetscCosReal(omega * t + th0);
443:     const PetscReal    st    = PetscSinReal(omega * t + th0);
444:     const PetscScalar *x     = &u[(p * 2 + 0) * dim];
445:     const PetscScalar *v     = &u[(p * 2 + 1) * dim];
446:     const PetscReal    xe[3] = {r0 * ct, r0 * st, 0.0};
447:     const PetscReal    ve[3] = {-v0 * st, v0 * ct, 0.0};
448:     PetscInt           d;

450:     for (d = 0; d < dim; ++d) {
451:       e[(p * 2 + 0) * dim + d] = x[d] - xe[d];
452:       e[(p * 2 + 1) * dim + d] = v[d] - ve[d];
453:     }
454:     if (user->error) {
455:       const PetscReal en   = 0.5 * DMPlex_DotRealD_Internal(dim, v, v);
456:       const PetscReal exen = 0.5 * PetscSqr(v0);
457:       PetscPrintf(comm, "t %.4g: p%" PetscInt_FMT " error [%.2g %.2g] 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));
458:     }
459:   }
460:   DMSwarmRestoreField(sw, "initCoordinates", NULL, NULL, (void **)&coords);
461:   DMSwarmRestoreField(sw, "initVelocity", NULL, NULL, (void **)&vel);
462:   VecRestoreArrayRead(U, &u);
463:   VecRestoreArray(E, &e);
464:   return 0;
465: }

467: static PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
468: {
469:   const PetscInt     ostep = ((AppCtx *)ctx)->ostep;
470:   DM                 sw;
471:   const PetscScalar *u;
472:   PetscInt           dim, Np, p;

475:   if (step % ostep == 0) {
476:     TSGetDM(ts, &sw);
477:     DMGetDimension(sw, &dim);
478:     VecGetArrayRead(U, &u);
479:     VecGetLocalSize(U, &Np);
480:     Np /= 2 * dim;
481:     if (!step) PetscPrintf(PetscObjectComm((PetscObject)ts), "Time     Step Part     Energy\n");
482:     for (p = 0; p < Np; ++p) {
483:       const PetscReal v2 = DMPlex_DotRealD_Internal(dim, &u[(p * 2 + 1) * dim], &u[(p * 2 + 1) * dim]);

485:       PetscSynchronizedPrintf(PetscObjectComm((PetscObject)ts), "%.6lf %4" PetscInt_FMT " %4" PetscInt_FMT " %10.4lf\n", (double)t, step, p, (double)(0.5 * v2));
486:     }
487:     PetscSynchronizedFlush(PetscObjectComm((PetscObject)ts), NULL);
488:     VecRestoreArrayRead(U, &u);
489:   }
490:   return 0;
491: }

493: static PetscErrorCode MigrateParticles(TS ts)
494: {
495:   DM sw;

498:   TSGetDM(ts, &sw);
499:   DMViewFromOptions(sw, NULL, "-migrate_view_pre");
500:   {
501:     Vec u, gc, gv;
502:     IS  isx, isv;

504:     TSGetSolution(ts, &u);
505:     TSRHSSplitGetIS(ts, "position", &isx);
506:     TSRHSSplitGetIS(ts, "momentum", &isv);
507:     DMSwarmCreateGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc);
508:     VecISCopy(u, isx, SCATTER_REVERSE, gc);
509:     DMSwarmDestroyGlobalVectorFromField(sw, DMSwarmPICField_coor, &gc);
510:     DMSwarmCreateGlobalVectorFromField(sw, "velocity", &gv);
511:     VecISCopy(u, isv, SCATTER_REVERSE, gv);
512:     DMSwarmDestroyGlobalVectorFromField(sw, "velocity", &gv);
513:   }
514:   DMSwarmMigrate(sw, PETSC_TRUE);
515:   DMSwarmTSRedistribute(ts);
516:   InitializeSolveAndSwarm(ts, PETSC_FALSE);
517:   //VecViewFromOptions(u, NULL, "-sol_migrate_view");
518:   return 0;
519: }

521: int main(int argc, char **argv)
522: {
523:   DM     dm, sw;
524:   TS     ts;
525:   Vec    u;
526:   AppCtx user;

529:   PetscInitialize(&argc, &argv, NULL, help);
530:   ProcessOptions(PETSC_COMM_WORLD, &user);
531:   CreateMesh(PETSC_COMM_WORLD, &user, &dm);
532:   CreateSwarm(dm, &user, &sw);
533:   DMSetApplicationContext(sw, &user);

535:   TSCreate(PETSC_COMM_WORLD, &ts);
536:   TSSetProblemType(ts, TS_NONLINEAR);
537:   TSSetDM(ts, sw);
538:   TSSetMaxTime(ts, 0.1);
539:   TSSetTimeStep(ts, 0.00001);
540:   TSSetMaxSteps(ts, 100);
541:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
542:   TSMonitorSet(ts, EnergyMonitor, &user, NULL);
543:   TSSetFromOptions(ts);
544:   TSSetComputeInitialCondition(ts, InitializeSolve);
545:   TSSetComputeExactError(ts, ComputeError);
546:   TSSetPostStep(ts, MigrateParticles);

548:   CreateSolution(ts);
549:   TSGetSolution(ts, &u);
550:   TSComputeInitialCondition(ts, u);
551:   TSSolve(ts, NULL);

553:   TSDestroy(&ts);
554:   DMDestroy(&sw);
555:   DMDestroy(&dm);
556:   PetscFinalize();
557:   return 0;
558: }

560: /*TEST

562:    build:
563:      requires: double !complex

565:    testset:
566:      requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
567:      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 \
568:            -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
569:            -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
570:            -dm_view -output_step 50 -error
571:      test:
572:        suffix: bsi_2d_1
573:        args: -ts_basicsymplectic_type 1
574:      test:
575:        suffix: bsi_2d_2
576:        args: -ts_basicsymplectic_type 2
577:      test:
578:        suffix: bsi_2d_3
579:        args: -ts_basicsymplectic_type 3
580:      test:
581:        suffix: bsi_2d_4
582:        args: -ts_basicsymplectic_type 4 -ts_dt 0.0001

584:    testset:
585:      requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
586:      args: -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
587:            -ts_type theta -ts_theta_theta 0.5 -ts_convergence_estimate -convest_num_refine 2 \
588:              -mat_type baij -ksp_error_if_not_converged -pc_type lu \
589:            -dm_view -output_step 50 -error
590:      test:
591:        suffix: im_2d_0
592:        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

594:    testset:
595:      requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
596:      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 \
597:            -dm_swarm_num_particles 2 -dm_swarm_coordinate_function circleSingleX -dm_swarm_velocity_function circleSingleV \
598:            -ts_type basicsymplectic -ts_convergence_estimate -convest_num_refine 2 \
599:            -dm_view -output_step 50
600:      test:
601:        suffix: bsi_2d_mesh_1
602:        args: -ts_basicsymplectic_type 1 -error
603:      test:
604:        suffix: bsi_2d_mesh_1_par_2
605:        nsize: 2
606:        args: -ts_basicsymplectic_type 1
607:      test:
608:        suffix: bsi_2d_mesh_1_par_3
609:        nsize: 3
610:        args: -ts_basicsymplectic_type 1
611:      test:
612:        suffix: bsi_2d_mesh_1_par_4
613:        nsize: 4
614:        args: -ts_basicsymplectic_type 1 -dm_swarm_num_particles 0,0,2,0

616:    testset:
617:      requires: defined(PETSC_HAVE_EXECUTABLE_EXPORT)
618:      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 \
619:            -dm_swarm_num_particles 10 -dm_swarm_coordinate_function circleMultipleX -dm_swarm_velocity_function circleMultipleV \
620:            -ts_convergence_estimate -convest_num_refine 2 \
621:            -dm_view -output_step 50 -error
622:      test:
623:        suffix: bsi_2d_multiple_1
624:        args: -ts_type basicsymplectic -ts_basicsymplectic_type 1
625:      test:
626:        suffix: bsi_2d_multiple_2
627:        args: -ts_type basicsymplectic -ts_basicsymplectic_type 2
628:      test:
629:        suffix: bsi_2d_multiple_3
630:        args: -ts_type basicsymplectic -ts_basicsymplectic_type 3 -ts_dt 0.001
631:      test:
632:        suffix: im_2d_multiple_0
633:        args: -ts_type theta -ts_theta_theta 0.5 \
634:                -mat_type baij -ksp_error_if_not_converged -pc_type lu

636: TEST*/