Actual source code: ex4.c
petsc-3.14.6 2021-03-30
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; /* 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: PetscReal omega; /* Oscillation frequency omega */
14: PetscInt particlesPerCell; /* The number of partices per cell */
15: PetscReal momentTol; /* Tolerance for checking moment conservation */
16: PetscBool monitor; /* Flag for using the TS monitor */
17: PetscBool error; /* Flag for printing the error */
18: PetscInt ostep; /* print the energy at each ostep time steps */
19: } AppCtx;
21: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
22: {
26: options->dim = 2;
27: options->simplex = PETSC_TRUE;
28: options->monitor = PETSC_FALSE;
29: options->error = PETSC_FALSE;
30: options->particlesPerCell = 1;
31: options->momentTol = 100.0*PETSC_MACHINE_EPSILON;
32: options->omega = 64.0;
33: options->ostep = 100;
35: PetscStrcpy(options->filename, "");
37: PetscOptionsBegin(comm, "", "Harmonic Oscillator Options", "DMPLEX");
38: PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL);
39: PetscOptionsInt("-dim", "The topological mesh dimension", "ex4.c", options->dim, &options->dim, NULL);
40: PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL);
41: PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL);
42: PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex4.c", options->simplex, &options->simplex, NULL);
43: PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex4.c", options->filename, options->filename, sizeof(options->filename), NULL);
44: PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL);
45: PetscOptionsReal("-omega", "Oscillator frequency", "ex4.c", options->omega, &options->omega, PETSC_NULL);
46: PetscOptionsEnd();
48: return(0);
49: }
51: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
52: {
53: PetscBool flg;
57: PetscStrcmp(user->filename, "", &flg);
58: if (flg) {
59: DMPlexCreateBoxMesh(comm, user->dim, user->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);
60: } else {
61: DMPlexCreateFromFile(comm, user->filename, PETSC_TRUE, dm);
62: DMGetDimension(*dm, &user->dim);
63: }
64: {
65: DM distributedMesh = NULL;
67: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
68: if (distributedMesh) {
69: DMDestroy(dm);
70: *dm = distributedMesh;
71: }
72: }
73: DMLocalizeCoordinates(*dm); /* needed for periodic */
74: DMSetFromOptions(*dm);
75: PetscObjectSetName((PetscObject) *dm, "Mesh");
76: DMViewFromOptions(*dm, NULL, "-dm_view");
77: return(0);
78: }
80: static PetscErrorCode SetInitialCoordinates(DM dmSw)
81: {
82: DM dm;
83: AppCtx *user;
84: PetscRandom rnd;
85: PetscBool simplex;
86: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
87: PetscInt dim, d, cStart, cEnd, c, Np, p;
91: PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd);
92: PetscRandomSetInterval(rnd, -1.0, 1.0);
93: PetscRandomSetFromOptions(rnd);
95: DMGetApplicationContext(dmSw, (void **) &user);
96: simplex = user->simplex;
97: Np = user->particlesPerCell;
98: DMSwarmGetCellDM(dmSw, &dm);
99: DMGetDimension(dm, &dim);
100: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
101: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
102: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
103: DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
104: for (c = cStart; c < cEnd; ++c) {
105: if (Np == 1) {
106: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
107: for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
108: } else {
109: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
110: for (p = 0; p < Np; ++p) {
111: const PetscInt n = c*Np + p;
112: PetscReal sum = 0.0, refcoords[3];
114: for (d = 0; d < dim; ++d) {
115: PetscRandomGetValueReal(rnd, &refcoords[d]);
116: sum += refcoords[d];
117: }
118: if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
119: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
120: }
121: }
122: }
123: DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
124: PetscFree5(centroid, xi0, v0, J, invJ);
125: PetscRandomDestroy(&rnd);
126: return(0);
127: }
129: static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
130: {
131: DM dm;
132: AppCtx *user;
133: PetscReal *coords;
134: PetscScalar *initialConditions;
135: PetscInt dim, cStart, cEnd, c, Np, p;
139: DMGetApplicationContext(dmSw, (void **) &user);
140: Np = user->particlesPerCell;
141: DMSwarmGetCellDM(dmSw, &dm);
142: DMGetDimension(dm, &dim);
143: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
144: DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
145: VecGetArray(u, &initialConditions);
146: for (c = cStart; c < cEnd; ++c) {
147: for (p = 0; p < Np; ++p) {
148: const PetscInt n = c*Np + p;
150: initialConditions[n*2+0] = DMPlex_NormD_Internal(dim, &coords[n*dim]);
151: initialConditions[n*2+1] = 0.0;
152: }
153: }
154: VecRestoreArray(u, &initialConditions);
155: DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
157: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
158: return(0);
159: }
161: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
162: {
163: PetscInt *cellid;
164: PetscInt dim, cStart, cEnd, c, Np = user->particlesPerCell, p;
168: DMGetDimension(dm, &dim);
169: DMCreate(PetscObjectComm((PetscObject) dm), sw);
170: DMSetType(*sw, DMSWARM);
171: DMSetDimension(*sw, dim);
173: DMSwarmSetType(*sw, DMSWARM_PIC);
174: DMSwarmSetCellDM(*sw, dm);
175: DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2, PETSC_REAL);
176: DMSwarmFinalizeFieldRegister(*sw);
177: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
178: DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);
179: DMSetFromOptions(*sw);
180: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
181: for (c = cStart; c < cEnd; ++c) {
182: for (p = 0; p < Np; ++p) {
183: const PetscInt n = c*Np + p;
185: cellid[n] = c;
186: }
187: }
188: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
189: PetscObjectSetName((PetscObject) *sw, "Particles");
190: DMViewFromOptions(*sw, NULL, "-sw_view");
191: return(0);
192: }
194: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
195: {
196: AppCtx *user = (AppCtx *) ctx;
197: const PetscReal omega = user->omega;
198: const PetscScalar *u;
199: MPI_Comm comm;
200: PetscReal dt;
201: PetscInt Np, p;
202: PetscErrorCode ierr;
205: if (step%user->ostep == 0) {
206: PetscObjectGetComm((PetscObject) ts, &comm);
207: if (!step) {PetscPrintf(comm, "Time Step Part Energy Mod Energy\n");}
208: TSGetTimeStep(ts, &dt);
209: VecGetArrayRead(U, &u);
210: VecGetLocalSize(U, &Np);
211: Np /= 2;
212: for (p = 0; p < Np; ++p) {
213: const PetscReal x = PetscRealPart(u[p*2+0]);
214: const PetscReal v = PetscRealPart(u[p*2+1]);
215: const PetscReal E = 0.5*(v*v + PetscSqr(omega)*x*x);
216: const PetscReal mE = 0.5*(v*v + PetscSqr(omega)*x*x - PetscSqr(omega)*dt*x*v);
218: PetscPrintf(comm, "%.6lf %4D %4D %10.4lf %10.4lf\n", t, step, p, (double) E, (double) mE);
219: }
220: VecRestoreArrayRead(U, &u);
221: }
222: return(0);
223: }
225: static PetscErrorCode InitializeSolve(TS ts, Vec u)
226: {
227: DM dm;
228: AppCtx *user;
232: TSGetDM(ts, &dm);
233: DMGetApplicationContext(dm, (void **) &user);
234: SetInitialCoordinates(dm);
235: SetInitialConditions(dm, u);
236: return(0);
237: }
239: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
240: {
241: MPI_Comm comm;
242: DM sdm;
243: AppCtx *user;
244: const PetscScalar *u, *coords;
245: PetscScalar *e;
246: PetscReal t, omega;
247: PetscInt dim, Np, p;
248: PetscErrorCode ierr;
252: PetscObjectGetComm((PetscObject) ts, &comm);
253: TSGetDM(ts, &sdm);
254: DMGetApplicationContext(sdm, (void **) &user);
255: omega = user->omega;
256: DMGetDimension(sdm, &dim);
257: TSGetSolveTime(ts, &t);
258: VecGetArray(E, &e);
259: VecGetArrayRead(U, &u);
260: VecGetLocalSize(U, &Np);
261: Np /= 2;
262: DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
263: for (p = 0; p < Np; ++p) {
264: const PetscReal x = PetscRealPart(u[p*2+0]);
265: const PetscReal v = PetscRealPart(u[p*2+1]);
266: const PetscReal x0 = DMPlex_NormD_Internal(dim, &coords[p*dim]);
267: const PetscReal ex = x0*PetscCosReal(omega*t);
268: const PetscReal ev = -x0*omega*PetscSinReal(omega*t);
270: 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));}
271: e[p*2+0] = x - ex;
272: e[p*2+1] = v - ev;
273: }
274: DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
276: VecRestoreArrayRead(U, &u);
277: VecRestoreArray(E, &e);
278: return(0);
279: }
282: /*---------------------Create particle RHS Functions--------------------------*/
283: static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
284: {
285: const PetscScalar *v;
286: PetscScalar *xres;
287: PetscInt Np, p;
288: PetscErrorCode ierr;
291: VecGetArray(Xres, &xres);
292: VecGetArrayRead(V, &v);
293: VecGetLocalSize(Xres, &Np);
294: for (p = 0; p < Np; ++p) {
295: xres[p] = v[p];
296: }
297: VecRestoreArrayRead(V, &v);
298: VecRestoreArray(Xres, &xres);
299: return(0);
300: }
302: static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
303: {
304: AppCtx *user = (AppCtx *)ctx;
305: const PetscScalar *x;
306: PetscScalar *vres;
307: PetscInt Np, p;
308: PetscErrorCode ierr;
311: VecGetArray(Vres, &vres);
312: VecGetArrayRead(X, &x);
313: VecGetLocalSize(Vres, &Np);
314: for (p = 0; p < Np; ++p) {
315: vres[p] = -PetscSqr(user->omega)*x[p];
316: }
317: VecRestoreArrayRead(X, &x);
318: VecRestoreArray(Vres, &vres);
319: return(0);
320: }
322: // static PetscErrorCode RHSFunctionParticles(TS ts, PetscReal t, Vec U, Vec R, void *ctx)
323: // {
324: // AppCtx *user = (AppCtx *) ctx;
325: // DM dm;
326: // const PetscScalar *u;
327: // PetscScalar *r;
328: // PetscInt Np, p;
329: // PetscErrorCode ierr;
330: //
332: // TSGetDM(ts, &dm);
333: // VecGetArray(R, &r);
334: // VecGetArrayRead(U, &u);
335: // VecGetLocalSize(U, &Np);
336: // Np /= 2;
337: // for (p = 0; p < Np; ++p) {
338: // r[p*2+0] = u[p*2+1];
339: // r[p*2+1] = -PetscSqr(user->omega)*u[p*2+0];
340: // }
341: // VecRestoreArrayRead(U, &u);
342: // VecRestoreArray(R, &r);
343: // return(0);
344: // }
345: /*----------------------------------------------------------------------------*/
347: /*--------------------Define RHSFunction, RHSJacobian (Theta)-----------------*/
348: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
349: {
350: AppCtx *user = (AppCtx *) ctx;
351: DM dm;
352: const PetscScalar *u;
353: PetscScalar *g;
354: PetscInt Np, p;
355: PetscErrorCode ierr;
358: TSGetDM(ts, &dm);
359: VecGetArray(G, &g);
360: VecGetArrayRead(U, &u);
361: VecGetLocalSize(U, &Np);
362: Np /= 2;
363: for (p = 0; p < Np; ++p) {
364: g[p*2+0] = u[p*2+1];
365: g[p*2+1] = -PetscSqr(user->omega)*u[p*2+0];
366: }
367: VecRestoreArrayRead(U, &u);
368: VecRestoreArray(G, &g);
369: return(0);
370: }
372: /*Ji = dFi/dxj
373: J= (0 1)
374: (-w^2 0)
375: */
376: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U , Mat J, Mat P, void *ctx)
377: {
378: AppCtx *user = (AppCtx *) ctx;
379: PetscInt Np = user->dim * user->particlesPerCell;
380: PetscInt i, m, n;
381: const PetscScalar *u;
382: PetscScalar vals[4] = {0., 1., -PetscSqr(user->omega), 0.};
383: PetscErrorCode ierr;
386: VecGetArrayRead(U, &u);
387: //PetscPrintf(PETSC_COMM_WORLD, "# Particles (Np) = %d \n" , Np);
389: MatGetOwnershipRange(J, &m, &n);
390: for (i = 0; i < Np; ++i) {
391: const PetscInt rows[2] = {2*i, 2*i+1};
392: MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES);
393: }
394: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
395: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
396: //MatView(J,PETSC_VIEWER_STDOUT_WORLD);
398: return(0);
400: }
401: /*----------------------------------------------------------------------------*/
403: /*----------------Define S, F, G Functions (Discrete Gradients)---------------*/
404: /*
405: u_t = S * gradF
406: --or--
407: u_t = S * G
409: + Sfunc - constructor for the S matrix from the formulation
410: . Ffunc - functional F from the formulation
411: - Gfunc - constructor for the gradient of F from the formulation
412: */
414: PetscErrorCode Sfunc(TS ts, PetscReal t, Vec U, Mat S, void *ctx)
415: {
416: AppCtx *user = (AppCtx *) ctx;
417: PetscInt Np = user->dim * user->particlesPerCell;
418: PetscInt i, m, n;
419: const PetscScalar *u;
420: PetscScalar vals[4] = {0., 1., -1, 0.};
421: PetscErrorCode ierr;
424: VecGetArrayRead(U, &u);
425: MatGetOwnershipRange(S, &m, &n);
426: for (i = 0; i < Np; ++i) {
427: const PetscInt rows[2] = {2*i, 2*i+1};
428: MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES);
429: }
430: MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);
431: MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);
433: return(0);
434: }
436: PetscErrorCode Ffunc(TS ts, PetscReal t, Vec U, PetscScalar *F, void *ctx)
437: {
438: AppCtx *user = (AppCtx *) ctx;
439: DM dm;
440: const PetscScalar *u;
441: PetscInt Np = user->dim * user->particlesPerCell;
442: PetscInt p;
443: PetscErrorCode ierr;
446: TSGetDM(ts, &dm);
447: //Define F
448: VecGetArrayRead(U, &u);
449: for (p = 0; p < Np; ++p) {
450: *F += (1/2)*PetscSqr(user->omega)*PetscSqr(u[p*2+0]) + (1/2)*PetscSqr(u[p*2+1]);
451: }
452: VecRestoreArrayRead(U, &u);
454: return(0);
455: }
457: PetscErrorCode gradFfunc(TS ts, PetscReal t, Vec U, Vec gradF, void *ctx)
458: {
459: AppCtx *user = (AppCtx *) ctx;
460: DM dm;
461: const PetscScalar *u;
462: PetscScalar *g;
463: PetscInt Np = user->dim * user->particlesPerCell;
464: PetscInt p;
465: PetscErrorCode ierr;
468: TSGetDM(ts, &dm);
469: VecGetArrayRead(U, &u);
471: //Define gradF
472: VecGetArray(gradF, &g);
473: for (p = 0; p < Np; ++p) {
474: g[p*2+0] = PetscSqr(user->omega)*u[p*2+0]; /*dF/dx*/
475: g[p*2+1] = u[p*2+1]; /*dF/dv*/
476: }
477: VecRestoreArrayRead(U, &u);
478: VecRestoreArray(gradF, &g);
480: return(0);
481: }
482: /*----------------------------------------------------------------------------*/
484: int main(int argc,char **argv)
485: {
486: TS ts; /* nonlinear solver */
487: DM dm, sw; /* Mesh and particle managers */
488: Vec u; /* swarm vector */
489: PetscInt n; /* swarm vector size */
490: IS is1, is2;
491: MPI_Comm comm;
492: Mat J; /* Jacobian matrix */
493: AppCtx user;
496: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
497: Initialize program
498: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
499: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
500: comm = PETSC_COMM_WORLD;
501: ProcessOptions(comm, &user);
504: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
505: Create Particle-Mesh
506: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
507: CreateMesh(comm, &dm, &user);
508: CreateParticles(dm, &sw, &user);
509: DMSetApplicationContext(sw, &user);
512: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
513: Setup timestepping solver
514: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
515: TSCreate(comm, &ts);
516: TSSetProblemType(ts,TS_NONLINEAR);
518: TSSetDM(ts, sw);
519: TSSetMaxTime(ts, 0.1);
520: TSSetTimeStep(ts, 0.00001);
521: TSSetMaxSteps(ts, 100);
522: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
523: if (user.monitor) {TSMonitorSet(ts, Monitor, &user, NULL);}
526: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
527: Prepare to solve
528: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
529: DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);
530: VecGetLocalSize(u, &n);
531: TSSetFromOptions(ts);
532: TSSetComputeInitialCondition(ts, InitializeSolve);
533: TSSetComputeExactError(ts, ComputeError);
535: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
536: Define function F(U, Udot , x , t) = G(U, x, t)
537: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
539: /* - - - - - - - Basic Symplectic - - - - - - - - - - - - - - - - - - - - - -*/
540: ISCreateStride(comm, n/2, 0, 2, &is1);
541: ISCreateStride(comm, n/2, 1, 2, &is2);
542: TSRHSSplitSetIS(ts, "position", is1);
543: TSRHSSplitSetIS(ts, "momentum", is2);
544: ISDestroy(&is1);
545: ISDestroy(&is2);
546: TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunction1, &user);
547: TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunction2, &user);
549: /* - - - - - - - Theta (Implicit Midpoint) - - - - - - - - - - - - - - - - -*/
550: TSSetRHSFunction(ts, NULL, RHSFunction, &user);
552: MatCreate(PETSC_COMM_WORLD,&J);
553: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
554: MatSetFromOptions(J);
555: MatSetUp(J);
556: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
557: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
558: PetscPrintf(comm, "n = %d\n",n);//Check number of particles
559: TSSetRHSJacobian(ts,J,J,RHSJacobian,&user);
561: /* - - - - - - - Discrete Gradients - - - - - - - - - - - - - - - - - - - - */
562: TSDiscGradSetFormulation(ts, Sfunc, Ffunc, gradFfunc, &user);
564: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
565: Solve
566: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
568: TSComputeInitialCondition(ts, u);
569: TSSolve(ts, u);
572: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
573: Clean up workspace
574: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
575: MatDestroy(&J);
576: DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);
577: TSDestroy(&ts);
578: DMDestroy(&sw);
579: DMDestroy(&dm);
580: PetscFinalize();
581: return ierr;
582: }
584: /*TEST
586: build:
587: requires: triangle !single !complex
588: test:
589: suffix: 1
590: 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
592: test:
593: suffix: 2
594: 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
596: test:
597: suffix: 3
598: 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
600: test:
601: suffix: 4
602: 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
604: test:
605: suffix: 5
606: 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
608: test:
609: suffix: 6
610: args: -dm_plex_box_faces 1,1 -ts_type discgrad -monitor -output_step 50 -error -ts_convergence_estimate -convest_num_refine 2
612: TEST*/