Actual source code: ex4.c
1: static char help[] = "Comparing basic symplectic, theta and discrete gradients interators on a simple hamiltonian system (harmonic oscillator) with particles\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: PetscInt numberOfCells; /* Number of cells in mesh */
16: PetscReal momentTol; /* Tolerance for checking moment conservation */
17: PetscBool monitor; /* Flag for using the TS monitor */
18: PetscBool error; /* Flag for printing the error */
19: PetscInt ostep; /* print the energy at each ostep time steps */
20: } AppCtx;
22: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
23: {
27: options->dim = 2;
28: options->simplex = PETSC_TRUE;
29: options->monitor = PETSC_FALSE;
30: options->error = PETSC_FALSE;
31: options->particlesPerCell = 1;
32: options->numberOfCells = 2;
33: options->momentTol = 100.0*PETSC_MACHINE_EPSILON;
34: options->omega = 64.0;
35: options->ostep = 100;
37: PetscStrcpy(options->filename, "");
39: PetscOptionsBegin(comm, "", "Harmonic Oscillator Options", "DMPLEX");
40: PetscOptionsInt("-output_step", "Number of time steps between output", "ex4.c", options->ostep, &options->ostep, PETSC_NULL);
41: PetscOptionsInt("-dim", "The topological mesh dimension", "ex4.c", options->dim, &options->dim, NULL);
42: PetscOptionsBool("-monitor", "Flag to use the TS monitor", "ex4.c", options->monitor, &options->monitor, NULL);
43: PetscOptionsBool("-error", "Flag to print the error", "ex4.c", options->error, &options->error, NULL);
44: PetscOptionsBool("-simplex", "The flag for simplices or tensor cells", "ex4.c", options->simplex, &options->simplex, NULL);
45: PetscOptionsString("-mesh", "Name of the mesh filename if any", "ex4.c", options->filename, options->filename, sizeof(options->filename), NULL);
46: PetscOptionsInt("-particles_per_cell", "Number of particles per cell", "ex4.c", options->particlesPerCell, &options->particlesPerCell, NULL);
47: PetscOptionsReal("-omega", "Oscillator frequency", "ex4.c", options->omega, &options->omega, PETSC_NULL);
48: PetscOptionsEnd();
50: return 0;
52: }
54: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *user)
55: {
57: DMCreate(comm, dm);
58: DMSetType(*dm, DMPLEX);
59: DMSetFromOptions(*dm);
60: DMViewFromOptions(*dm, NULL, "-dm_view");
62: return 0;
63: }
65: static PetscErrorCode SetInitialCoordinates(DM dmSw)
66: {
67: DM dm;
68: AppCtx *user;
69: PetscRandom rnd;
70: PetscBool simplex;
71: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
72: PetscInt dim, d, cStart, cEnd, c, Np, p;
75: PetscRandomCreate(PetscObjectComm((PetscObject) dmSw), &rnd);
76: PetscRandomSetInterval(rnd, -1.0, 1.0);
77: PetscRandomSetFromOptions(rnd);
79: DMGetApplicationContext(dmSw, &user);
80: Np = user->particlesPerCell;
81: DMSwarmGetCellDM(dmSw, &dm);
82: DMGetDimension(dm, &dim);
83: DMPlexIsSimplex(dm, &simplex);
84: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
85: user->numberOfCells = cEnd - cStart;
86: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
87: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
88: DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
89: for (c = cStart; c < cEnd; ++c) {
90: if (Np == 1) {
91: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
92: for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
93: } else {
94: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ); /* affine */
95: for (p = 0; p < Np; ++p) {
96: const PetscInt n = c*Np + p;
97: PetscReal sum = 0.0, refcoords[3];
99: for (d = 0; d < dim; ++d) {
100: PetscRandomGetValueReal(rnd, &refcoords[d]);
101: sum += refcoords[d];
102: }
103: if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
104: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
105: }
106: }
107: }
109: DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
110: PetscFree5(centroid, xi0, v0, J, invJ);
111: PetscRandomDestroy(&rnd);
112: return 0;
113: }
115: static PetscErrorCode SetInitialConditions(DM dmSw, Vec u)
116: {
117: DM dm;
118: AppCtx *user;
119: PetscReal *coords;
120: PetscScalar *initialConditions;
121: PetscInt dim, cStart, cEnd, c, Np, p;
124: DMGetApplicationContext(dmSw, &user);
125: Np = user->particlesPerCell;
126: DMSwarmGetCellDM(dmSw, &dm);
127: DMGetDimension(dm, &dim);
128: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
129: DMSwarmGetField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
130: VecGetArray(u, &initialConditions);
131: for (c = cStart; c < cEnd; ++c) {
132: for (p = 0; p < Np; ++p) {
133: const PetscInt n = c*Np + p;
134: initialConditions[n*2+0] = DMPlex_NormD_Internal(dim, &coords[n*dim]);
135: initialConditions[n*2+1] = 0.0;
136: }
137: }
138: VecRestoreArray(u, &initialConditions);
139: DMSwarmRestoreField(dmSw, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
140: return 0;
141: }
143: static PetscErrorCode CreateParticles(DM dm, DM *sw, AppCtx *user)
144: {
145: PetscInt *cellid;
146: PetscInt dim, cStart, cEnd, c, Np, p;
149: DMGetDimension(dm, &dim);
150: DMCreate(PetscObjectComm((PetscObject) dm), sw);
151: DMSetType(*sw, DMSWARM);
152: DMSetDimension(*sw, dim);
153: Np = user->particlesPerCell;
154: DMSwarmSetType(*sw, DMSWARM_PIC);
155: DMSwarmSetCellDM(*sw, dm);
156: DMSwarmRegisterPetscDatatypeField(*sw, "kinematics", 2, PETSC_REAL);
157: DMSwarmFinalizeFieldRegister(*sw);
158: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
159: DMSwarmSetLocalSizes(*sw, (cEnd - cStart) * Np, 0);
160: DMSetFromOptions(*sw);
161: DMSwarmGetField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
162: for (c = cStart; c < cEnd; ++c) {
163: for (p = 0; p < Np; ++p) {
164: const PetscInt n = c*Np + p;
166: cellid[n] = c;
167: }
168: }
169: DMSwarmRestoreField(*sw, DMSwarmPICField_cellid, NULL, NULL, (void **) &cellid);
170: PetscObjectSetName((PetscObject) *sw, "Particles");
171: DMViewFromOptions(*sw, NULL, "-sw_view");
172: return 0;
173: }
175: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, void *ctx)
176: {
177: AppCtx *user = (AppCtx *) ctx;
178: const PetscReal omega = user->omega;
179: const PetscScalar *u;
180: MPI_Comm comm;
181: PetscReal dt;
182: PetscInt Np, p;
185: if (step%user->ostep == 0) {
186: PetscObjectGetComm((PetscObject) ts, &comm);
187: if (!step) PetscPrintf(comm, "Time Step Part Energy Mod Energy\n");
188: TSGetTimeStep(ts, &dt);
189: VecGetArrayRead(U, &u);
190: VecGetLocalSize(U, &Np);
191: Np /= 2;
192: for (p = 0; p < Np; ++p) {
193: const PetscReal x = PetscRealPart(u[p*2+0]);
194: const PetscReal v = PetscRealPart(u[p*2+1]);
195: const PetscReal E = 0.5*(v*v + PetscSqr(omega)*x*x);
196: const PetscReal mE = 0.5*(v*v + PetscSqr(omega)*x*x - PetscSqr(omega)*dt*x*v);
197: PetscPrintf(comm, "%.6lf %4D %4D %10.4lf %10.4lf\n", t, step, p, (double) E, (double) mE);
198: }
199: VecRestoreArrayRead(U, &u);
200: }
201: return 0;
202: }
204: static PetscErrorCode InitializeSolve(TS ts, Vec u)
205: {
206: DM dm;
207: AppCtx *user;
210: TSGetDM(ts, &dm);
211: DMGetApplicationContext(dm, &user);
212: SetInitialCoordinates(dm);
213: SetInitialConditions(dm, u);
214: return 0;
215: }
217: static PetscErrorCode ComputeError(TS ts, Vec U, Vec E)
218: {
219: MPI_Comm comm;
220: DM sdm;
221: AppCtx *user;
222: const PetscScalar *u, *coords;
223: PetscScalar *e;
224: PetscReal t, omega;
225: PetscInt dim, Np, p;
228: PetscObjectGetComm((PetscObject) ts, &comm);
229: TSGetDM(ts, &sdm);
230: DMGetApplicationContext(sdm, &user);
231: omega = user->omega;
232: DMGetDimension(sdm, &dim);
233: TSGetSolveTime(ts, &t);
234: VecGetArray(E, &e);
235: VecGetArrayRead(U, &u);
236: VecGetLocalSize(U, &Np);
237: Np /= 2;
238: DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
239: for (p = 0; p < Np; ++p) {
240: const PetscReal x = PetscRealPart(u[p*2+0]);
241: const PetscReal v = PetscRealPart(u[p*2+1]);
242: const PetscReal x0 = DMPlex_NormD_Internal(dim, &coords[p*dim]);
243: const PetscReal ex = x0*PetscCosReal(omega*t);
244: const PetscReal ev = -x0*omega*PetscSinReal(omega*t);
246: 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));
247: e[p*2+0] = x - ex;
248: e[p*2+1] = v - ev;
249: }
250: DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
252: VecRestoreArrayRead(U, &u);
253: VecRestoreArray(E, &e);
254: return 0;
255: }
257: /*---------------------Create particle RHS Functions--------------------------*/
258: static PetscErrorCode RHSFunction1(TS ts, PetscReal t, Vec V, Vec Xres, void *ctx)
259: {
260: const PetscScalar *v;
261: PetscScalar *xres;
262: PetscInt Np, p;
265: VecGetArray(Xres, &xres);
266: VecGetArrayRead(V, &v);
267: VecGetLocalSize(Xres, &Np);
268: for (p = 0; p < Np; ++p) {
269: xres[p] = v[p];
270: }
271: VecRestoreArrayRead(V, &v);
272: VecRestoreArray(Xres, &xres);
273: return 0;
274: }
276: static PetscErrorCode RHSFunction2(TS ts, PetscReal t, Vec X, Vec Vres, void *ctx)
277: {
278: AppCtx *user = (AppCtx *)ctx;
279: const PetscScalar *x;
280: PetscScalar *vres;
281: PetscInt Np, p;
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: /*----------------------------------------------------------------------------*/
297: /*--------------------Define RHSFunction, RHSJacobian (Theta)-----------------*/
298: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec G, void *ctx)
299: {
300: AppCtx *user = (AppCtx *) ctx;
301: DM dm;
302: const PetscScalar *u;
303: PetscScalar *g;
304: PetscInt Np, p;
307: TSGetDM(ts, &dm);
308: VecGetArray(G, &g);
309: VecGetArrayRead(U, &u);
310: VecGetLocalSize(U, &Np);
311: Np /= 2;
312: for (p = 0; p < Np; ++p) {
313: g[p*2+0] = u[p*2+1];
314: g[p*2+1] = -PetscSqr(user->omega)*u[p*2+0];
315: }
316: VecRestoreArrayRead(U, &u);
317: VecRestoreArray(G, &g);
318: return 0;
319: }
321: /*Ji = dFi/dxj
322: J= (0 1)
323: (-w^2 0)
324: */
325: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U , Mat J, Mat P, void *ctx)
326: {
327: AppCtx *user = (AppCtx *) ctx;
328: PetscInt Np;
329: PetscInt i, m, n;
330: const PetscScalar *u;
331: PetscScalar vals[4] = {0., 1., -PetscSqr(user->omega), 0.};
334: VecGetArrayRead(U, &u);
335: VecGetLocalSize(U, &Np);
336: Np /= 2;
337: MatGetOwnershipRange(J, &m, &n);
338: for (i = 0; i < Np; ++i) {
339: const PetscInt rows[2] = {2*i, 2*i+1};
340: MatSetValues(J, 2, rows, 2, rows, vals, INSERT_VALUES);
341: }
342: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
343: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
344: return 0;
346: }
348: /*----------------------------------------------------------------------------*/
350: /*----------------Define S, F, G Functions (Discrete Gradients)---------------*/
351: /*
352: u_t = S * gradF
353: --or--
354: u_t = S * G
356: + Sfunc - constructor for the S matrix from the formulation
357: . Ffunc - functional F from the formulation
358: - Gfunc - constructor for the gradient of F from the formulation
359: */
361: PetscErrorCode Sfunc(TS ts, PetscReal t, Vec U, Mat S, void *ctx)
362: {
363: PetscInt Np;
364: PetscInt i, m, n;
365: const PetscScalar *u;
366: PetscScalar vals[4] = {0., 1., -1, 0.};
369: VecGetArrayRead(U, &u);
370: VecGetLocalSize(U, &Np);
371: Np /= 2;
372: MatGetOwnershipRange(S, &m, &n);
373: for (i = 0; i < Np; ++i) {
374: const PetscInt rows[2] = {2*i, 2*i+1};
375: MatSetValues(S, 2, rows, 2, rows, vals, INSERT_VALUES);
376: }
377: MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY);
378: MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY);
379: return 0;
381: }
383: PetscErrorCode Ffunc(TS ts, PetscReal t, Vec U, PetscScalar *F, void *ctx)
384: {
385: AppCtx *user = (AppCtx *) ctx;
386: DM dm;
387: const PetscScalar *u;
388: PetscInt Np;
389: PetscInt p;
392: TSGetDM(ts, &dm);
394: /* Define F */
395: VecGetArrayRead(U, &u);
396: VecGetLocalSize(U, &Np);
397: Np /= 2;
398: for (p = 0; p < Np; ++p) {
399: *F += 0.5*PetscSqr(user->omega)*PetscSqr(u[p*2+0]) + 0.5*PetscSqr(u[p*2+1]);
400: }
401: VecRestoreArrayRead(U, &u);
402: return 0;
403: }
405: PetscErrorCode gradFfunc(TS ts, PetscReal t, Vec U, Vec gradF, void *ctx)
406: {
407: AppCtx *user = (AppCtx *) ctx;
408: DM dm;
409: const PetscScalar *u;
410: PetscScalar *g;
411: PetscInt Np;
412: PetscInt p;
415: TSGetDM(ts, &dm);
416: VecGetArrayRead(U, &u);
417: VecGetLocalSize(U, &Np);
418: Np /= 2;
420: VecGetArray(gradF, &g);
421: for (p = 0; p < Np; ++p) {
422: g[p*2+0] = PetscSqr(user->omega)*u[p*2+0]; /*dF/dx*/
423: g[p*2+1] = u[p*2+1]; /*dF/dv*/
424: }
425: VecRestoreArrayRead(U, &u);
426: VecRestoreArray(gradF, &g);
427: return 0;
428: }
430: /*----------------------------------------------------------------------------*/
432: int main(int argc,char **argv)
433: {
434: TS ts; /* nonlinear solver */
435: DM dm, sw; /* Mesh and particle managers */
436: Vec u; /* swarm vector */
437: PetscInt n; /* swarm vector size */
438: IS is1, is2;
439: MPI_Comm comm;
440: Mat J; /* Jacobian matrix */
441: AppCtx user;
443: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
444: Initialize program
445: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
446: PetscInitialize(&argc, &argv, NULL, help);
447: comm = PETSC_COMM_WORLD;
448: ProcessOptions(comm, &user);
450: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
451: Create Particle-Mesh
452: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
453: CreateMesh(comm, &dm, &user);
454: CreateParticles(dm, &sw, &user);
455: DMSetApplicationContext(sw, &user);
457: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
458: Setup timestepping solver
459: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
460: TSCreate(comm, &ts);
461: TSSetProblemType(ts,TS_NONLINEAR);
463: TSSetDM(ts, sw);
464: TSSetMaxTime(ts, 0.1);
465: TSSetTimeStep(ts, 0.00001);
466: TSSetMaxSteps(ts, 100);
467: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
468: if (user.monitor) TSMonitorSet(ts, Monitor, &user, NULL);
470: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
471: Prepare to solve
472: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
473: DMSwarmCreateGlobalVectorFromField(sw, "kinematics", &u);
474: VecGetLocalSize(u, &n);
475: TSSetFromOptions(ts);
476: TSSetComputeInitialCondition(ts, InitializeSolve);
477: TSSetComputeExactError(ts, ComputeError);
479: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
480: Define function F(U, Udot , x , t) = G(U, x, t)
481: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
483: /* - - - - - - - Basic Symplectic - - - - - - - - - - - - - - - - - - - - - -*/
484: ISCreateStride(comm, n/2, 0, 2, &is1);
485: ISCreateStride(comm, n/2, 1, 2, &is2);
486: TSRHSSplitSetIS(ts, "position", is1);
487: TSRHSSplitSetIS(ts, "momentum", is2);
488: ISDestroy(&is1);
489: ISDestroy(&is2);
490: TSRHSSplitSetRHSFunction(ts, "position", NULL, RHSFunction1, &user);
491: TSRHSSplitSetRHSFunction(ts, "momentum", NULL, RHSFunction2, &user);
493: /* - - - - - - - Theta (Implicit Midpoint) - - - - - - - - - - - - - - - - -*/
495: TSSetRHSFunction(ts, NULL, RHSFunction, &user);
497: MatCreate(PETSC_COMM_WORLD,&J);
498: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
499: MatSetFromOptions(J);
500: MatSetUp(J);
501: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
502: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
503: TSSetRHSJacobian(ts,J,J,RHSJacobian,&user);
505: /* - - - - - - - Discrete Gradients - - - - - - - - - - - - - - - - - - - - */
506: TSDiscGradSetFormulation(ts, Sfunc, Ffunc, gradFfunc, &user);
508: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
509: Solve
510: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
511: TSComputeInitialCondition(ts, u);
512: TSSolve(ts, u);
514: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
515: Clean up workspace
516: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
517: MatDestroy(&J);
518: DMSwarmDestroyGlobalVectorFromField(sw, "kinematics", &u);
519: TSDestroy(&ts);
520: DMDestroy(&sw);
521: DMDestroy(&dm);
522: PetscFinalize();
523: return 0;
524: }
526: /*TEST
528: build:
529: requires: triangle !single !complex
530: test:
531: suffix: 1
532: args: -dm_plex_box_faces 1,1 -ts_type basicsymplectic -ts_basicsymplectic_type 1 -ts_convergence_estimate -convest_num_refine 2 -dm_view -monitor -output_step 50 -error
534: test:
535: suffix: 2
536: args: -dm_plex_box_faces 1,1 -ts_type basicsymplectic -ts_basicsymplectic_type 2 -ts_convergence_estimate -convest_num_refine 2 -dm_view -monitor -output_step 50 -error
538: test:
539: suffix: 3
540: 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 -monitor -output_step 50 -error
542: test:
543: suffix: 4
544: 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 -monitor -output_step 50 -error
546: test:
547: suffix: 5
548: args: -dm_plex_box_faces 1,1 -ts_type theta -ts_theta_theta 0.5 -monitor -output_step 50 -error -ts_convergence_estimate -convest_num_refine 2 -dm_view
550: test:
551: suffix: 6
552: args: -dm_plex_box_faces 1,1 -ts_type discgrad -monitor -output_step 50 -error -ts_convergence_estimate -convest_num_refine 2 -dm_view
554: test:
555: suffix: 7
556: args: -dm_plex_box_faces 1,1 -ts_type discgrad -ts_discgrad_gonzalez -monitor -output_step 50 -error -ts_convergence_estimate -convest_num_refine 2 -dm_view
558: TEST*/