xref: /petsc/src/dm/impls/swarm/tests/ex5.c (revision 6a98f8dc3f2c9149905a87dc2e9d0fedaf64e09a)
1 static char help[] = "Vlasov example of particles orbiting around a central massive point.\n";
2 
3 #include <petscdmplex.h>
4 #include <petsc/private/dmpleximpl.h>  /* For norm */
5 #include <petsc/private/petscfeimpl.h> /* For CoordinatesRefToReal() */
6 #include <petscdmswarm.h>
7 #include <petscts.h>
8 
9 typedef struct {
10   PetscInt  dim;                          /* The topological mesh dimension */
11   PetscBool simplex;                      /* Flag for simplices or tensor cells */
12   char      filename[PETSC_MAX_PATH_LEN]; /* Name of the mesh filename if any */
13   PetscInt  particlesPerCell;             /* The number of partices per cell */
14   PetscReal momentTol;                    /* Tolerance for checking moment conservation */
15   PetscBool monitor;                      /* Flag for using the TS monitor */
16   PetscBool error;                        /* Flag for printing the error */
17   PetscInt  ostep;                        /* print the energy at each ostep time steps */
18 } AppCtx;
19 
20 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
21 {
22   PetscErrorCode ierr;
23 
24   PetscFunctionBeginUser;
25   options->dim              = 2;
26   options->simplex          = PETSC_TRUE;
27   options->monitor          = PETSC_FALSE;
28   options->error            = PETSC_FALSE;
29   options->particlesPerCell = 1;
30   options->momentTol        = 100.0*PETSC_MACHINE_EPSILON;
31   options->ostep            = 100;
32 
33   ierr = PetscStrcpy(options->filename, "");CHKERRQ(ierr);
34 
35   ierr = PetscOptionsBegin(comm, "", "Vlasov Options", "DMPLEX");CHKERRQ(ierr);
36   ierr = PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL);CHKERRQ(ierr);
37   ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex4.c", options->dim, &options->dim, NULL);CHKERRQ(ierr);
38   ierr = PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL);CHKERRQ(ierr);
39   ierr = PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL);CHKERRQ(ierr);
40   ierr = PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex4.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr);
41   ierr = PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex4.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr);
42   ierr = PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL);CHKERRQ(ierr);
43   ierr = PetscOptionsEnd();CHKERRQ(ierr);
44 
45   PetscFunctionReturn(0);
46 }
47 
48 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
49 {
50   PetscBool      flg;
51   PetscErrorCode ierr;
52 
53   PetscFunctionBeginUser;
54   ierr = PetscStrcmp(user->filename, "", &flg);CHKERRQ(ierr);
55   if (flg) {
56     ierr = DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);
57   } else {
58     ierr = DMPlexCreateFromFile(comm, user->filename, PETSC_TRUE, dm);CHKERRQ(ierr);
59     ierr = DMGetDimension(*dm, &user->dim);CHKERRQ(ierr);
60   }
61   {
62     DM distributedMesh = NULL;
63 
64     ierr = DMPlexDistribute(*dm, 0, NULL, &distributedMesh);CHKERRQ(ierr);
65     if (distributedMesh) {
66       ierr = DMDestroy(dm);CHKERRQ(ierr);
67       *dm  = distributedMesh;
68     }
69   }
70   ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); /* needed for periodic */
71   ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
72   ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
73   ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
74   PetscFunctionReturn(0);
75 }
76 
77 static PetscErrorCode SetInitialCoordinates(DM dmSw)
78 {
79   DM             dm;
80   AppCtx        *user;
81   PetscRandom    rnd;
82   PetscBool      simplex;
83   PetscReal     *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
84   PetscInt       dim, d, cStart, cEnd, c, Np, p;
85   PetscErrorCode ierr;
86 
87   PetscFunctionBeginUser;
88   ierr = PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd);CHKERRQ(ierr);
89   ierr = PetscRandomSetInterval(rnd, -1.0, 1.0);CHKERRQ(ierr);
90   ierr = PetscRandomSetFromOptions(rnd);CHKERRQ(ierr);
91 
92   ierr = DMGetApplicationContext(dmSw, (void **) &user);CHKERRQ(ierr);
93   simplex = user->simplex;
94   Np   = user->particlesPerCell;
95   ierr = DMSwarmGetCellDM(dmSw, &dm);CHKERRQ(ierr);
96   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
97   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
98   ierr = PetscMalloc5(dim, &centroid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);CHKERRQ(ierr);
99   for (d = 0; d < dim; ++d) xi0[d] = -1.0;
100   ierr = DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr);
101   for (c = cStart; c < cEnd; ++c) {
102     if (Np == 1) {
103       ierr = DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);CHKERRQ(ierr);
104       for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
105     } else {
106       ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); /* affine */
107       for (p = 0; p < Np; ++p) {
108         const PetscInt n   = c*Np + p;
109         PetscReal      sum = 0.0, refcoords[3];
110 
111         for (d = 0; d < dim; ++d) {
112           ierr = PetscRandomGetValueReal(rnd, &refcoords[d]);CHKERRQ(ierr);
113           sum += refcoords[d];
114         }
115         if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
116         CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
117       }
118     }
119   }
120   ierr = DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr);
121   ierr = PetscFree5(centroid, xi0, v0, J, invJ);CHKERRQ(ierr);
122   ierr = PetscRandomDestroy(&rnd);CHKERRQ(ierr);
123   PetscFunctionReturn(0);
124 }
125 
126 static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
127 {
128   DM             dm;
129   AppCtx        *user;
130   PetscScalar   *initialConditions;
131   PetscInt       dim, cStart, cEnd, c, Np, p;
132   PetscErrorCode ierr;
133 
134   PetscFunctionBeginUser;
135   ierr = DMGetApplicationContext(dmSw, (void **) &user);CHKERRQ(ierr);
136   Np   = user->particlesPerCell;
137   ierr = DMSwarmGetCellDM(dmSw, &dm);CHKERRQ(ierr);
138   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
139   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
140   ierr = VecGetArray(u, &initialConditions);CHKERRQ(ierr);
141   for (c = cStart; c < cEnd; ++c) {
142     for (p = 0; p < Np; ++p) {
143       const PetscInt n = c*Np + p;
144 
145       initialConditions[(n*2 + 0)*dim + 0] = n+1;
146       initialConditions[(n*2 + 0)*dim + 1] = 0;
147       initialConditions[(n*2 + 1)*dim + 0] = 0;
148       initialConditions[(n*2 + 1)*dim + 1] = PetscSqrtReal(1000./(n+1.));
149     }
150   }
151   ierr = VecRestoreArray(u, &initialConditions);CHKERRQ(ierr);
152   PetscFunctionReturn(0);
153 }
154 
155 static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
156 {
157   PetscInt      *cellid;
158   PetscInt       dim, cStart, cEnd, c, Np = user->particlesPerCell, p;
159   PetscErrorCode ierr;
160 
161   PetscFunctionBeginUser;
162   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
163   ierr = DMCreate(PetscObjectComm((PetscObject) dm), sw);CHKERRQ(ierr);
164   ierr = DMSetType(*sw, DMSWARM);CHKERRQ(ierr);
165   ierr = DMSetDimension(*sw, dim);CHKERRQ(ierr);
166 
167   ierr = DMSwarmSetType(*sw, DMSWARM_PIC);CHKERRQ(ierr);
168   ierr = DMSwarmSetCellDM(*sw, dm);CHKERRQ(ierr);
169   ierr = DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2*dim, PETSC_REAL);CHKERRQ(ierr);
170   ierr = DMSwarmFinalizeFieldRegister(*sw);CHKERRQ(ierr);
171   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
172   ierr = DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);CHKERRQ(ierr);
173   ierr = DMSetFromOptions(*sw);CHKERRQ(ierr);
174   ierr = DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr);
175   for (c = cStart; c < cEnd; ++c) {
176     for (p = 0; p < Np; ++p) {
177       const PetscInt n = c*Np + p;
178 
179       cellid[n] = c;
180     }
181   }
182   ierr = DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);CHKERRQ(ierr);
183   ierr = PetscObjectSetName((PetscObject) *sw, "Particles");CHKERRQ(ierr);
184   ierr = DMViewFromOptions(*sw, NULL, "-sw_view");CHKERRQ(ierr);
185   PetscFunctionReturn(0);
186 }
187 
188 /* Create particle RHS Functions for gravity with G = 1 for simplification */
189 static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
190 {
191   DM                dm;
192   const PetscScalar *v;
193   PetscScalar       *xres;
194   PetscInt          Np, p, dim, d;
195   PetscErrorCode    ierr;
196 
197   PetscFunctionBeginUser;
198   ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
199   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
200   ierr = VecGetLocalSize(Xres, &Np);CHKERRQ(ierr);
201   Np  /= dim;
202   ierr = VecGetArray(Xres, &xres);CHKERRQ(ierr);
203   ierr = VecGetArrayRead(V, &v);CHKERRQ(ierr);
204   for (p = 0; p < Np; ++p) {
205      for(d = 0; d < dim; ++d) {
206        xres[p*dim+d] = v[p*dim+d];
207      }
208   }
209   ierr = VecRestoreArrayRead(V,& v);CHKERRQ(ierr);
210   ierr = VecRestoreArray(Xres, &xres);CHKERRQ(ierr);
211   PetscFunctionReturn(0);
212 }
213 
214 static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
215 {
216   DM                dm;
217   const PetscScalar *x;
218   PetscScalar       *vres;
219   PetscInt          Np, p, dim, d;
220   PetscErrorCode    ierr;
221 
222 
223   PetscFunctionBeginUser;
224   ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
225   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
226   ierr = VecGetLocalSize(Vres, &Np);CHKERRQ(ierr);
227   Np  /= dim;
228   ierr = VecGetArray(Vres, &vres);CHKERRQ(ierr);
229   ierr = VecGetArrayRead(X, &x);CHKERRQ(ierr);
230   for (p = 0; p < Np; ++p) {
231     const PetscScalar rsqr = DMPlex_NormD_Internal(dim, &x[p*dim]);
232 
233     for (d = 0; d < dim; ++d) {
234       vres[p*dim+d] = -(1000./(p+1.)) * x[p*dim+d]/PetscSqr(rsqr);
235     }
236   }
237   ierr = VecRestoreArray(Vres, &vres);CHKERRQ(ierr);
238   ierr = VecRestoreArrayRead(X, &x);CHKERRQ(ierr);
239   PetscFunctionReturn(0);
240 }
241 
242 static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t , Vec U, Vec R, void *ctx)
243 {
244   DM                dm;
245   const PetscScalar *u;
246   PetscScalar       *r;
247   PetscInt          Np, p, dim, d;
248   PetscErrorCode    ierr;
249 
250   PetscFunctionBeginUser;
251   ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
252   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
253   ierr = VecGetLocalSize(U, &Np);CHKERRQ(ierr);
254   Np  /= 2*dim;
255   ierr = VecGetArrayRead(U, &u);CHKERRQ(ierr);
256   ierr = VecGetArray(R, &r);CHKERRQ(ierr);
257   for( p = 0; p < Np; ++p) {
258     const PetscScalar rsqr = DMPlex_NormD_Internal(dim, &u[p*2*dim]);
259 
260     for (d = 0; d < dim; ++d) {
261         r[p*2*dim+d]   = u[p*2*dim+d+2];
262         r[p*2*dim+d+2] = -(1000./(1.+p)) * u[p*2*dim+d]/PetscSqr(rsqr);
263     }
264   }
265   ierr = VecRestoreArrayRead(U, &u);CHKERRQ(ierr);
266   ierr = VecRestoreArray(R, &r);CHKERRQ(ierr);
267   PetscFunctionReturn(0);
268 }
269 
270 static PetscErrorCode InitializeSolve(TS ts, Vec u)
271 {
272   DM             dm;
273   AppCtx        *user;
274   PetscErrorCode ierr;
275 
276   PetscFunctionBeginUser;
277   ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
278   ierr = DMGetApplicationContext(dm, (void **) &user);CHKERRQ(ierr);
279   ierr = SetInitialCoordinates(dm);CHKERRQ(ierr);
280   ierr = SetInitialConditions(dm, u);CHKERRQ(ierr);
281   PetscFunctionReturn(0);
282 }
283 
284 static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
285 {
286   MPI_Comm           comm;
287   DM                 sdm;
288   AppCtx            *user;
289   const PetscScalar *u, *coords;
290   PetscScalar       *e;
291   PetscReal          t;
292   PetscInt           dim, Np, p;
293   PetscErrorCode     ierr;
294 
295   PetscFunctionBeginUser;
296   ierr = PetscObjectGetComm((PetscObject) ts, &comm);CHKERRQ(ierr);
297   ierr = TSGetDM(ts, &sdm);CHKERRQ(ierr);
298   ierr = DMGetApplicationContext(sdm, (void **) &user);CHKERRQ(ierr);
299   ierr = DMGetDimension(sdm, &dim);CHKERRQ(ierr);
300   ierr = TSGetSolveTime(ts, &t);CHKERRQ(ierr);
301   ierr = VecGetArray(E, &e);CHKERRQ(ierr);
302   ierr = VecGetArrayRead(U, &u);CHKERRQ(ierr);
303   ierr = VecGetLocalSize(U, &Np);CHKERRQ(ierr);
304   Np  /= 2*dim;
305   ierr = DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr);
306   for (p = 0; p < Np; ++p) {
307     const PetscScalar *x     = &u[(p*2+0)*dim];
308     const PetscScalar *v     = &u[(p*2+1)*dim];
309     const PetscReal    x0    = p+1.;
310     const PetscReal    omega = PetscSqrtReal(1000./(p+1.))/x0;
311     const PetscReal    xe[3] = { x0*PetscCosReal(omega*t),       x0*PetscSinReal(omega*t),       0.0};
312     const PetscReal    ve[3] = {-x0*omega*PetscSinReal(omega*t), x0*omega*PetscCosReal(omega*t), 0.0};
313     PetscInt           d;
314 
315     for (d = 0; d < dim; ++d) {
316       e[(p*2+0)*dim+d] = x[d] - xe[d];
317       e[(p*2+1)*dim+d] = v[d] - ve[d];
318     }
319     if (user->error) {ierr = PetscPrintf(comm, "p%D error [%.2g %.2g] sol [(%.6lf %.6lf) (%.6lf %.6lf)] exact [(%.6lf %.6lf) (%.6lf %.6lf)] energy/exact energy %g / %g\n", 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], 0.5*DMPlex_NormD_Internal(dim, v), (double) (0.5*(1000./(p+1))));}
320   }
321   ierr = DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);CHKERRQ(ierr);
322   ierr = VecRestoreArrayRead(U, &u);CHKERRQ(ierr);
323   ierr = VecRestoreArray(E, &e);CHKERRQ(ierr);
324   PetscFunctionReturn(0);
325 }
326 
327 int main(int argc, char **argv)
328 {
329   TS                 ts;
330   TSConvergedReason  reason;
331   DM                 dm, sw;
332   Vec                u;
333   IS                 is1, is2;
334   PetscInt          *idx1, *idx2;
335   MPI_Comm           comm;
336   AppCtx             user;
337   const PetscScalar *endVals;
338   PetscReal          ftime   = .1;
339   PetscInt           locSize, dim, d, Np, p, steps;
340   PetscErrorCode     ierr;
341 
342   ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
343   comm = PETSC_COMM_WORLD;
344   ierr = ProcessOptions(comm, &user);CHKERRQ(ierr);
345 
346   ierr = CreateMesh(comm, &dm, &user);CHKERRQ(ierr);
347   ierr = CreateParticles(dm, &sw, &user);CHKERRQ(ierr);
348   ierr = DMSetApplicationContext(sw, &user);CHKERRQ(ierr);
349 
350   ierr = TSCreate(comm, &ts);CHKERRQ(ierr);
351   ierr = TSSetType(ts, TSBASICSYMPLECTIC);CHKERRQ(ierr);
352   ierr = TSSetDM(ts, sw);CHKERRQ(ierr);
353   ierr = TSSetMaxTime(ts, ftime);CHKERRQ(ierr);
354   ierr = TSSetTimeStep(ts, 0.0001);CHKERRQ(ierr);
355   ierr = TSSetMaxSteps(ts, 10);CHKERRQ(ierr);
356   ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
357   ierr = TSSetTime(ts, 0.0);CHKERRQ(ierr);
358   ierr = TSSetRHSFunction(ts, NULL, RHSFunctionParticles, &user);CHKERRQ(ierr);
359 
360   ierr = DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);CHKERRQ(ierr);
361   ierr = DMGetDimension(sw, &dim);CHKERRQ(ierr);
362   ierr = VecGetLocalSize(u, &locSize);CHKERRQ(ierr);
363   Np   = locSize/(2*dim);
364   ierr = PetscMalloc1(locSize/2, &idx1);CHKERRQ(ierr);
365   ierr = PetscMalloc1(locSize/2, &idx2);CHKERRQ(ierr);
366   for (p = 0; p < Np; ++p) {
367     for (d = 0; d < dim; ++d) {
368       idx1[p*dim+d] = (p*2+0)*dim + d;
369       idx2[p*dim+d] = (p*2+1)*dim + d;
370     }
371   }
372   ierr = ISCreateGeneral(comm, locSize/2, idx1, PETSC_OWN_POINTER, &is1);CHKERRQ(ierr);
373   ierr = ISCreateGeneral(comm, locSize/2, idx2, PETSC_OWN_POINTER, &is2);CHKERRQ(ierr);
374   ierr = TSRHSSplitSetIS(ts, "position", is1);CHKERRQ(ierr);
375   ierr = TSRHSSplitSetIS(ts, "momentum", is2);CHKERRQ(ierr);
376   ierr = ISDestroy(&is1);CHKERRQ(ierr);
377   ierr = ISDestroy(&is2);CHKERRQ(ierr);
378   ierr = TSRHSSplitSetRHSFunction(ts,"position",NULL,RHSFunction1,&user);CHKERRQ(ierr);
379   ierr = TSRHSSplitSetRHSFunction(ts,"momentum",NULL,RHSFunction2,&user);CHKERRQ(ierr);
380 
381   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
382   ierr = TSSetComputeInitialCondition(ts, InitializeSolve);CHKERRQ(ierr);
383   ierr = TSSetComputeExactError(ts, ComputeError);CHKERRQ(ierr);
384   ierr = TSComputeInitialCondition(ts, u);CHKERRQ(ierr);
385   ierr = VecViewFromOptions(u, NULL, "-init_view");CHKERRQ(ierr);
386   ierr = TSSolve(ts, u);CHKERRQ(ierr);
387   ierr = TSGetSolveTime(ts, &ftime);CHKERRQ(ierr);
388   ierr = TSGetConvergedReason(ts, &reason);CHKERRQ(ierr);
389   ierr = TSGetStepNumber(ts, &steps);CHKERRQ(ierr);
390   ierr = PetscPrintf(comm,"%s at time %g after %D steps\n", TSConvergedReasons[reason], (double) ftime, steps);CHKERRQ(ierr);
391 
392   ierr = VecGetArrayRead(u, &endVals);CHKERRQ(ierr);
393   for (p = 0; p < Np; ++p) {
394     const PetscReal norm = DMPlex_NormD_Internal(dim, &endVals[(p*2 + 1)*dim]);
395     ierr = PetscPrintf(comm, "Particle %D initial Energy: %g  Final Energy: %g\n", p, (double) (0.5*(1000./(p+1))), (double) 0.5*PetscSqr(norm));CHKERRQ(ierr);
396   }
397   ierr = VecRestoreArrayRead(u, &endVals);CHKERRQ(ierr);
398   ierr = DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);CHKERRQ(ierr);
399   ierr = TSDestroy(&ts);CHKERRQ(ierr);
400   ierr = DMDestroy(&sw);CHKERRQ(ierr);
401   ierr = DMDestroy(&dm);CHKERRQ(ierr);
402   ierr = PetscFinalize();
403   return ierr;
404 }
405 
406 
407 /*TEST
408 
409    build:
410      requires: triangle !single !complex
411    test:
412      suffix: bsi1
413      args: -dim 2 -dm_plex_box_faces 1,1 -dm_view -sw_view -ts_basicsymplectic_type 1 -ts_max_time 0.1 -ts_monitor_sp_swarm -ts_convergence_estimate -convest_num_refine 2
414    test:
415      suffix: bsi2
416      args: -dim 2 -dm_plex_box_faces 1,1 -dm_view -sw_view -ts_basicsymplectic_type 2 -ts_max_time 0.1 -ts_monitor_sp_swarm -ts_convergence_estimate -convest_num_refine 2
417    test:
418      suffix: euler
419      args: -dim 2 -dm_plex_box_faces 1,1 -dm_view -sw_view -ts_type euler -ts_max_time 0.1 -ts_monitor_sp_swarm -ts_convergence_estimate -convest_num_refine 2
420 
421 TEST*/
422