Actual source code: ex4.c

  1: static char help[] = "Example of simple hamiltonian system (harmonic oscillator) with particles and a basic symplectic integrator\n";

  3: #include <petscdmplex.h>
  4: #include <petsc/private/dmpleximpl.h>
  5: #include <petsc/private/petscfeimpl.h>
  6: #include <petscdmswarm.h>
  7: #include <petscts.h>

  9: typedef struct {
 10:   PetscInt  dim;
 11:   PetscReal omega;            /* Oscillation frequency omega */
 12:   PetscInt  particlesPerCell; /* The number of partices per cell */
 13:   PetscReal momentTol;        /* Tolerance for checking moment conservation */
 14:   PetscBool monitor;          /* Flag for using the TS monitor */
 15:   PetscBool error;            /* Flag for printing the error */
 16:   PetscInt  ostep;            /* print the energy at each ostep time steps */
 17: } AppCtx;

 19: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 20: {

 24:   options->monitor          = PETSC_FALSE;
 25:   options->error            = PETSC_FALSE;
 26:   options->particlesPerCell = 1;
 27:   options->momentTol        = 100.0*PETSC_MACHINE_EPSILON;
 28:   options->omega            = 64.0;
 29:   options->ostep            = 100;

 31:   PetscOptionsBegin(comm, "", "Harmonic Oscillator Options", "DMPLEX");
 32:   PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL);
 33:   PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL);
 34:   PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL);
 35:   PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL);
 36:   PetscOptionsReal("-omega", "Oscillator frequency", "ex4.c", options->omega, &options->omega, PETSC_NULL);
 37:   PetscOptionsEnd();

 39:   return(0);
 40: }

 42: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
 43: {

 47:   DMCreate(comm, dm);
 48:   DMSetType(*dm, DMPLEX);
 49:   DMSetFromOptions(*dm);
 50:   DMViewFromOptions(*dm, NULL, "-dm_view");
 51:   DMGetDimension(*dm, &user->dim);
 52:   return(0);
 53: }

 55: static PetscErrorCode SetInitialCoordinates(DM dmSw)
 56: {
 57:   DM             dm;
 58:   AppCtx        *user;
 59:   PetscRandom    rnd;
 60:   PetscBool      simplex;
 61:   PetscReal     *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
 62:   PetscInt       dim, d, cStart, cEnd, c, Np, p;

 66:   PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd);
 67:   PetscRandomSetInterval(rnd, -1.0, 1.0);
 68:   PetscRandomSetFromOptions(rnd);

 70:   DMGetApplicationContext(dmSw, &user);
 71:   Np   = user->particlesPerCell;
 72:   DMSwarmGetCellDM(dmSw, &dm);
 73:   DMGetDimension(dm, &dim);
 74:   DMPlexIsSimplex(dm, &simplex);
 75:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
 76:   PetscMalloc5(dim, &centroid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
 77:   for (d = 0; d < dim; ++d) xi0[d] = -1.0;
 78:   DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
 79:   for (c = cStart; c < cEnd; ++c) {
 80:     if (Np == 1) {
 81:       DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
 82:       for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
 83:     } else {
 84:       DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
 85:       for (p = 0; p < Np; ++p) {
 86:         const PetscInt n   = c*Np + p;
 87:         PetscReal      sum = 0.0, refcoords[3];

 89:         for (d = 0; d < dim; ++d) {
 90:           PetscRandomGetValueReal(rnd, &refcoords[d]);
 91:           sum += refcoords[d];
 92:         }
 93:         if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
 94:         CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
 95:       }
 96:     }
 97:   }
 98:   DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
 99:   PetscFree5(centroid, xi0, v0, J, invJ);
100:   PetscRandomDestroy(&rnd);
101:   return(0);
102: }

104: static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
105: {
106:   DM             dm;
107:   AppCtx        *user;
108:   PetscReal     *coords;
109:   PetscScalar   *initialConditions;
110:   PetscInt       dim, cStart, cEnd, c, Np, p;

114:   DMGetApplicationContext(dmSw, &user);
115:   Np   = user->particlesPerCell;
116:   DMSwarmGetCellDM(dmSw, &dm);
117:   DMGetDimension(dm, &dim);
118:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
119:   DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
120:   VecGetArray(u, &initialConditions);
121:   for (c = cStart; c < cEnd; ++c) {
122:     for (p = 0; p < Np; ++p) {
123:       const PetscInt n = c*Np + p;

125:       initialConditions[n*2+0] = DMPlex_NormD_Internal(dim, &coords[n*dim]);
126:       initialConditions[n*2+1] = 0.0;
127:     }
128:   }
129:   VecRestoreArray(u, &initialConditions);
130:   DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);

132:   VecView(u, PETSC_VIEWER_STDOUT_WORLD);
133:   return(0);
134: }

136: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
137: {
138:   PetscInt      *cellid;
139:   PetscInt       dim, cStart, cEnd, c, Np = user->particlesPerCell, p;

143:   DMGetDimension(dm, &dim);
144:   DMCreate(PetscObjectComm((PetscObject) dm), sw);
145:   DMSetType(*sw, DMSWARM);
146:   DMSetDimension(*sw, dim);

148:   DMSwarmSetType(*sw, DMSWARM_PIC);
149:   DMSwarmSetCellDM(*sw, dm);
150:   DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2, PETSC_REAL);
151:   DMSwarmFinalizeFieldRegister(*sw);
152:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
153:   DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);
154:   DMSetFromOptions(*sw);
155:   DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
156:   for (c = cStart; c < cEnd; ++c) {
157:     for (p = 0; p < Np; ++p) {
158:       const PetscInt n = c*Np + p;

160:       cellid[n] = c;
161:     }
162:   }
163:   DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
164:   PetscObjectSetName((PetscObject) *sw, "Particles");
165:   DMViewFromOptions(*sw, NULL, "-sw_view");
166:   return(0);
167: }

169: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
170: {
171:   AppCtx            *user  = (AppCtx *) ctx;
172:   const PetscReal    omega = user->omega;
173:   const PetscScalar *u;
174:   MPI_Comm           comm;
175:   PetscReal          dt;
176:   PetscInt           Np, p;
177:   PetscErrorCode     ierr;

180:   if (step%user->ostep == 0) {
181:     PetscObjectGetComm((PetscObject) ts, &comm);
182:     if (!step) {PetscPrintf(comm, "Time     Step Part     Energy Mod Energy\n");}
183:     TSGetTimeStep(ts, &dt);
184:     VecGetArrayRead(U, &u);
185:     VecGetLocalSize(U, &Np);
186:     Np /= 2;
187:     for (p = 0; p < Np; ++p) {
188:       const PetscReal x  = PetscRealPart(u[p*2+0]);
189:       const PetscReal v  = PetscRealPart(u[p*2+1]);
190:       const PetscReal E  = 0.5*(v*v + PetscSqr(omega)*x*x);
191:       const PetscReal mE = 0.5*(v*v + PetscSqr(omega)*x*x - PetscSqr(omega)*dt*x*v);

193:       PetscPrintf(comm, "%.6lf %4D %4D %10.4lf %10.4lf\n", t, step, p, (double) E, (double) mE);
194:     }
195:     VecRestoreArrayRead(U, &u);
196:   }
197:   return(0);
198: }

200: static PetscErrorCode InitializeSolve(TS ts, Vec u)
201: {
202:   DM             dm;
203:   AppCtx        *user;

207:   TSGetDM(ts, &dm);
208:   DMGetApplicationContext(dm, &user);
209:   SetInitialCoordinates(dm);
210:   SetInitialConditions(dm, u);
211:   return(0);
212: }

214: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
215: {
216:   MPI_Comm           comm;
217:   DM                 sdm;
218:   AppCtx            *user;
219:   const PetscScalar *u, *coords;
220:   PetscScalar       *e;
221:   PetscReal          t, omega;
222:   PetscInt           dim, Np, p;
223:   PetscErrorCode     ierr;

226:   PetscObjectGetComm((PetscObject) ts, &comm);
227:   TSGetDM(ts, &sdm);
228:   DMGetApplicationContext(sdm, &user);
229:   omega = user->omega;
230:   DMGetDimension(sdm, &dim);
231:   TSGetSolveTime(ts, &t);
232:   VecGetArray(E, &e);
233:   VecGetArrayRead(U, &u);
234:   VecGetLocalSize(U, &Np);
235:   Np  /= 2;
236:   DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
237:   for (p = 0; p < Np; ++p) {
238:     const PetscReal x  = PetscRealPart(u[p*2+0]);
239:     const PetscReal v  = PetscRealPart(u[p*2+1]);
240:     const PetscReal x0 = DMPlex_NormD_Internal(dim, &coords[p*dim]);
241:     const PetscReal ex =  x0*PetscCosReal(omega*t);
242:     const PetscReal ev = -x0*omega*PetscSinReal(omega*t);

244:     if (user->error) {PetscPrintf(comm, "p%D error [%.2g %.2g] sol [%.6lf %.6lf] exact [%.6lf %.6lf] energy/exact energy %g / %g\n", p, (double) PetscAbsReal(x-ex), (double) PetscAbsReal(v-ev), (double) x, (double) v, (double) ex, (double) ev, 0.5*(v*v + PetscSqr(omega)*x*x), (double) 0.5*PetscSqr(omega*x0));}
245:     e[p*2+0] = x - ex;
246:     e[p*2+1] = v - ev;
247:   }
248:   DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);

250:   VecRestoreArrayRead(U, &u);
251:   VecRestoreArray(E, &e);
252:   return(0);
253: }

255: /*---------------------Create particle RHS Functions--------------------------*/
256: static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
257: {
258:   const PetscScalar *v;
259:   PetscScalar       *xres;
260:   PetscInt           Np, p;
261:   PetscErrorCode     ierr;

264:   VecGetArray(Xres, &xres);
265:   VecGetArrayRead(V, &v);
266:   VecGetLocalSize(Xres, &Np);
267:   for (p = 0; p < Np; ++p) {
268:      xres[p] = v[p];
269:   }
270:   VecRestoreArrayRead(V, &v);
271:   VecRestoreArray(Xres, &xres);
272:   return(0);
273: }

275: static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
276: {
277:   AppCtx            *user = (AppCtx *)ctx;
278:   const PetscScalar *x;
279:   PetscScalar       *vres;
280:   PetscInt           Np, p;
281:   PetscErrorCode     ierr;

284:   VecGetArray(Vres, &vres);
285:   VecGetArrayRead(X, &x);
286:   VecGetLocalSize(Vres, &Np);
287:   for (p = 0; p < Np; ++p) {
288:     vres[p] = -PetscSqr(user->omega)*x[p];
289:   }
290:   VecRestoreArrayRead(X, &x);
291:   VecRestoreArray(Vres, &vres);
292:   return(0);
293: }

295: // static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
296: // {
297: //   AppCtx            *user = (AppCtx *) ctx;
298: //   DM                 dm;
299: //   const PetscScalar *u;
300: //   PetscScalar       *r;
301: //   PetscInt           Np, p;
302: //   PetscErrorCode     ierr;
303: //
305: //   TSGetDM(ts, &dm);
306: //   VecGetArray(R, &r);
307: //   VecGetArrayRead(U, &u);
308: //   VecGetLocalSize(U, &Np);
309: //   Np  /= 2;
310: //   for (p = 0; p < Np; ++p) {
311: //     r[p*2+0] = u[p*2+1];
312: //     r[p*2+1] = -PetscSqr(user->omega)*u[p*2+0];
313: //   }
314: //   VecRestoreArrayRead(U, &u);
315: //   VecRestoreArray(R, &r);
316: //   return(0);
317: // }
318: /*----------------------------------------------------------------------------*/

320: /*--------------------Define RHSFunction, RHSJacobian (Theta)-----------------*/
321: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
322: {
323:   AppCtx            *user = (AppCtx *) ctx;
324:   DM                 dm;
325:   const PetscScalar *u;
326:   PetscScalar       *g;
327:   PetscInt           Np, p;
328:   PetscErrorCode     ierr;

331:   TSGetDM(ts, &dm);
332:   VecGetArray(G, &g);
333:   VecGetArrayRead(U, &u);
334:   VecGetLocalSize(U, &Np);
335:   Np  /= 2;
336:   for (p = 0; p < Np; ++p) {
337:     g[p*2+0] = u[p*2+1];
338:     g[p*2+1] = -PetscSqr(user->omega)*u[p*2+0];
339:   }
340:   VecRestoreArrayRead(U, &u);
341:   VecRestoreArray(G, &g);
342:   return(0);
343: }

345: /*Ji = dFi/dxj
346: J= (0    1)
347:    (-w^2 0)
348: */
349: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U , Mat J, Mat P, void *ctx)
350: {
351:   AppCtx            *user = (AppCtx *) ctx;
352:   PetscInt           Np = user->dim * user->particlesPerCell;
353:   PetscInt           i, m, n;
354:   const PetscScalar *u;
355:   PetscScalar        vals[4] = {0., 1., -PetscSqr(user->omega), 0.};
356:   PetscErrorCode     ierr;

359:   VecGetArrayRead(U, &u);
360:   //PetscPrintf(PETSC_COMM_WORLD, "# Particles (Np) = %d \n" , Np);

362:   MatGetOwnershipRange(J, &m, &n);
363:   for (i = 0; i < Np; ++i) {
364:     const PetscInt rows[2] = {2*i, 2*i+1};
365:     MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES);
366:   }
367:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
368:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
369:   //MatView(J,PETSC_VIEWER_STDOUT_WORLD);

371:   return(0);

373: }
374: /*----------------------------------------------------------------------------*/

376: /*----------------Define S, F, G Functions (Discrete Gradients)---------------*/
377: /*
378:   u_t = S * gradF
379:     --or--
380:   u_t = S * G

382:   + Sfunc - constructor for the S matrix from the formulation
383:   . Ffunc - functional F from the formulation
384:   - Gfunc - constructor for the gradient of F from the formulation
385: */

387: PetscErrorCode Sfunc(TS ts, PetscReal t, Vec U, Mat S, void *ctx)
388: {
389:   AppCtx            *user = (AppCtx *) ctx;
390:   PetscInt           Np = user->dim * user->particlesPerCell;
391:   PetscInt           i, m, n;
392:   const PetscScalar *u;
393:   PetscScalar        vals[4] = {0., 1., -1, 0.};
394:   PetscErrorCode     ierr;

397:   VecGetArrayRead(U, &u);
398:   MatGetOwnershipRange(S, &m, &n);
399:   for (i = 0; i < Np; ++i) {
400:     const PetscInt rows[2] = {2*i, 2*i+1};
401:     MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES);
402:   }
403:   MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);
404:   MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);

406:   return(0);
407: }

409: PetscErrorCode Ffunc(TS ts, PetscReal t, Vec U, PetscScalar *F, void *ctx)
410: {
411:   AppCtx            *user = (AppCtx *) ctx;
412:   DM                 dm;
413:   const PetscScalar *u;
414:   PetscInt           Np = user->dim * user->particlesPerCell;
415:   PetscInt           p;
416:   PetscErrorCode     ierr;

419:   TSGetDM(ts, &dm);
420:   //Define F
421:   VecGetArrayRead(U, &u);
422:   for (p = 0; p < Np; ++p) {
423:     *F += (1/2)*PetscSqr(user->omega)*PetscSqr(u[p*2+0]) + (1/2)*PetscSqr(u[p*2+1]);
424:   }
425:   VecRestoreArrayRead(U, &u);

427:   return(0);
428: }

430: PetscErrorCode gradFfunc(TS ts, PetscReal t, Vec U, Vec gradF, void *ctx)
431: {
432:   AppCtx            *user = (AppCtx *) ctx;
433:   DM                 dm;
434:   const PetscScalar *u;
435:   PetscScalar       *g;
436:   PetscInt           Np = user->dim * user->particlesPerCell;
437:   PetscInt           p;
438:   PetscErrorCode     ierr;

441:   TSGetDM(ts, &dm);
442:   VecGetArrayRead(U, &u);

444:   //Define gradF
445:   VecGetArray(gradF, &g);
446:   for (p = 0; p < Np; ++p) {
447:     g[p*2+0] = PetscSqr(user->omega)*u[p*2+0]; /*dF/dx*/
448:     g[p*2+1] = u[p*2+1]; /*dF/dv*/
449:   }
450:   VecRestoreArrayRead(U, &u);
451:   VecRestoreArray(gradF, &g);

453:   return(0);
454: }
455: /*----------------------------------------------------------------------------*/

457: int main(int argc,char **argv)
458: {
459:   TS             ts;     /* nonlinear solver                 */
460:   DM             dm, sw; /* Mesh and particle managers       */
461:   Vec            u;      /* swarm vector                     */
462:   PetscInt       n;      /* swarm vector size                */
463:   IS             is1, is2;
464:   MPI_Comm       comm;
465:   Mat            J;      /* Jacobian matrix                  */
466:   AppCtx         user;

469:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
470:      Initialize program
471:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
472:   PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
473:   comm = PETSC_COMM_WORLD;
474:   ProcessOptions(comm, &user);

476:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
477:      Create Particle-Mesh
478:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
479:   CreateMesh(comm, &dm, &user);
480:   CreateParticles(dm, &sw, &user);
481:   DMSetApplicationContext(sw, &user);

483:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
484:      Setup timestepping solver
485:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
486:   TSCreate(comm, &ts);
487:   TSSetProblemType(ts,TS_NONLINEAR);

489:   TSSetDM(ts, sw);
490:   TSSetMaxTime(ts, 0.1);
491:   TSSetTimeStep(ts, 0.00001);
492:   TSSetMaxSteps(ts, 100);
493:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
494:   if (user.monitor) {TSMonitorSet(ts, Monitor, &user, NULL);}

496:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
497:      Prepare to solve
498:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
499:   DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);
500:   VecGetLocalSize(u, &n);
501:   TSSetFromOptions(ts);
502:   TSSetComputeInitialCondition(ts, InitializeSolve);
503:   TSSetComputeExactError(ts, ComputeError);

505:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
506:      Define function F(U, Udot , x , t) = G(U, x, t)
507:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

509:   /* - - - - - - - Basic Symplectic - - - - - - - - - - - - - - - - - - - - - -*/
510:   ISCreateStride(comm, n/2, 0, 2, &is1);
511:   ISCreateStride(comm, n/2, 1, 2, &is2);
512:   TSRHSSplitSetIS(ts, "position", is1);
513:   TSRHSSplitSetIS(ts, "momentum", is2);
514:   ISDestroy(&is1);
515:   ISDestroy(&is2);
516:   TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunction1, &user);
517:   TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunction2, &user);

519:   /* - - - - - - - Theta (Implicit Midpoint) - - - - - - - - - - - - - - - - -*/
520:   TSSetRHSFunction(ts, NULL, RHSFunction, &user);

522:   MatCreate(PETSC_COMM_WORLD,&J);
523:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
524:   MatSetFromOptions(J);
525:   MatSetUp(J);
526:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
527:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
528:   PetscPrintf(comm, "n = %d\n",n);//Check number of particles
529:   TSSetRHSJacobian(ts,J,J,RHSJacobian,&user);

531:   /* - - - - - - - Discrete Gradients - - - - - - - - - - - - - - - - - - - - */
532:   TSDiscGradSetFormulation(ts, Sfunc, Ffunc, gradFfunc, &user);

534:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
535:      Solve
536:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

538:   TSComputeInitialCondition(ts, u);
539:   TSSolve(ts, u);

541:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
542:      Clean up workspace
543:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
544:   MatDestroy(&J);
545:   DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);
546:   TSDestroy(&ts);
547:   DMDestroy(&sw);
548:   DMDestroy(&dm);
549:   PetscFinalize();
550:   return ierr;
551: }

553: /*TEST

555:    build:
556:      requires: triangle !single !complex
557:    test:
558:      suffix: 1
559:      args: -dm_plex_box_faces 1,1 -ts_type basicsymplectic -ts_basicsymplectic_type 1 -ts_convergence_estimate -convest_num_refine 2 -dm_view -sw_view -monitor -output_step 50 -error

561:    test:
562:      suffix: 2
563:      args: -dm_plex_box_faces 1,1 -ts_type basicsymplectic -ts_basicsymplectic_type 2 -ts_convergence_estimate -convest_num_refine 2 -dm_view -sw_view -monitor -output_step 50 -error

565:    test:
566:      suffix: 3
567:      args: -dm_plex_box_faces 1,1 -ts_type basicsymplectic -ts_basicsymplectic_type 3 -ts_convergence_estimate -convest_num_refine 2 -ts_dt 0.0001 -dm_view -sw_view -monitor -output_step 50 -error

569:    test:
570:      suffix: 4
571:      args: -dm_plex_box_faces 1,1 -ts_type basicsymplectic -ts_basicsymplectic_type 4 -ts_convergence_estimate -convest_num_refine 2 -ts_dt 0.0001 -dm_view -sw_view -monitor -output_step 50 -error

573:    test:
574:      suffix: 5
575:      args: -dm_plex_box_faces 1,1 -ts_type theta -ts_theta_theta 0.5 -ts_convergence_estimate -convest_num_refine 2 -dm_view -sw_view -monitor -output_step 50 -error

577:    test:
578:      suffix: 6
579:      args: -dm_plex_box_faces 1,1 -ts_type discgrad -monitor -output_step 50 -error -ts_convergence_estimate -convest_num_refine 2

581: TEST*/