Actual source code: ex62.c
petsc-3.5.4 2015-05-23
1: static char help[] = "Stokes Problem in 2d and 3d with simplicial finite elements.\n\
2: We solve the Stokes problem in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";
5: /*
6: The isoviscous Stokes problem, which we discretize using the finite
7: element method on an unstructured mesh. The weak form equations are
9: < \nabla v, \nabla u + {\nabla u}^T > - < \nabla\cdot v, p > + < v, f > = 0
10: < q, \nabla\cdot u > = 0
12: We start with homogeneous Dirichlet conditions. We will expand this as the set
13: of test problems is developed.
15: Discretization:
17: We use PetscFE to generate a tabulation of the finite element basis functions
18: at quadrature points. We can currently generate an arbitrary order Lagrange
19: element.
21: Field Data:
23: DMPLEX data is organized by point, and the closure operation just stacks up the
24: data from each sieve point in the closure. Thus, for a P_2-P_1 Stokes element, we
25: have
27: cl{e} = {f e_0 e_1 e_2 v_0 v_1 v_2}
28: x = [u_{e_0} v_{e_0} u_{e_1} v_{e_1} u_{e_2} v_{e_2} u_{v_0} v_{v_0} p_{v_0} u_{v_1} v_{v_1} p_{v_1} u_{v_2} v_{v_2} p_{v_2}]
30: The problem here is that we would like to loop over each field separately for
31: integration. Therefore, the closure visitor in DMPlexVecGetClosure() reorders
32: the data so that each field is contiguous
34: x' = [u_{e_0} v_{e_0} u_{e_1} v_{e_1} u_{e_2} v_{e_2} u_{v_0} v_{v_0} u_{v_1} v_{v_1} u_{v_2} v_{v_2} p_{v_0} p_{v_1} p_{v_2}]
36: Likewise, DMPlexVecSetClosure() takes data partitioned by field, and correctly
37: puts it into the Sieve ordering.
39: Next Steps:
41: - Refine and show convergence of correct order automatically (use femTest.py)
42: - Fix InitialGuess for arbitrary disc (means making dual application work again)
43: - Redo slides from GUCASTutorial for this new example
45: For tensor product meshes, see SNES ex67, ex72
46: */
48: #include <petscdmplex.h>
49: #include <petscsnes.h>
50: #include <petscds.h>
52: PetscInt spatialDim = 0;
54: typedef enum {NEUMANN, DIRICHLET} BCType;
55: typedef enum {RUN_FULL, RUN_TEST} RunType;
57: typedef struct {
58: PetscInt debug; /* The debugging level */
59: RunType runType; /* Whether to run tests, or solve the full problem */
60: PetscLogEvent createMeshEvent;
61: PetscBool showInitial, showSolution, showError;
62: /* Domain and mesh definition */
63: PetscInt dim; /* The topological mesh dimension */
64: PetscBool interpolate; /* Generate intermediate mesh elements */
65: PetscBool simplex; /* Use simplices or tensor product cells */
66: PetscReal refinementLimit; /* The largest allowable cell volume */
67: char partitioner[2048]; /* The graph partitioner */
68: /* Problem definition */
69: BCType bcType;
70: void (**exactFuncs)(const PetscReal x[], PetscScalar *u, void *ctx);
71: } AppCtx;
73: void zero_scalar(const PetscReal coords[], PetscScalar *u, void *ctx)
74: {
75: u[0] = 0.0;
76: }
77: void zero_vector(const PetscReal coords[], PetscScalar *u, void *ctx)
78: {
79: PetscInt d;
80: for (d = 0; d < spatialDim; ++d) u[d] = 0.0;
81: }
83: /*
84: In 2D we use exact solution:
86: u = x^2 + y^2
87: v = 2 x^2 - 2xy
88: p = x + y - 1
89: f_x = f_y = 3
91: so that
93: -\Delta u + \nabla p + f = <-4, -4> + <1, 1> + <3, 3> = 0
94: \nabla \cdot u = 2x - 2x = 0
95: */
96: void quadratic_u_2d(const PetscReal x[], PetscScalar *u, void *ctx)
97: {
98: u[0] = x[0]*x[0] + x[1]*x[1];
99: u[1] = 2.0*x[0]*x[0] - 2.0*x[0]*x[1];
100: }
102: void linear_p_2d(const PetscReal x[], PetscScalar *p, void *ctx)
103: {
104: *p = x[0] + x[1] - 1.0;
105: }
106: void constant_p(const PetscReal x[], PetscScalar *p, void *ctx)
107: {
108: *p = 1.0;
109: }
111: void f0_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[])
112: {
113: PetscInt c;
114: for (c = 0; c < spatialDim; ++c) f0[c] = 3.0;
115: }
117: /* gradU[comp*dim+d] = {u_x, u_y, v_x, v_y} or {u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z}
118: u[Ncomp] = {p} */
119: void f1_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[])
120: {
121: const PetscInt dim = spatialDim;
122: const PetscInt Ncomp = spatialDim;
123: PetscInt comp, d;
125: for (comp = 0; comp < Ncomp; ++comp) {
126: for (d = 0; d < dim; ++d) {
127: /* f1[comp*dim+d] = 0.5*(gradU[comp*dim+d] + gradU[d*dim+comp]); */
128: f1[comp*dim+d] = u_x[comp*dim+d];
129: }
130: f1[comp*dim+comp] -= u[Ncomp];
131: }
132: }
134: /* gradU[comp*dim+d] = {u_x, u_y, v_x, v_y} or {u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z} */
135: void f0_p(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[])
136: {
137: const PetscInt dim = spatialDim;
138: PetscInt d;
140: f0[0] = 0.0;
141: for (d = 0; d < dim; ++d) f0[0] += u_x[d*dim+d];
142: }
144: void f1_p(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[])
145: {
146: PetscInt d;
147: for (d = 0; d < spatialDim; ++d) f1[d] = 0.0;
148: }
150: /* < q, \nabla\cdot u >
151: NcompI = 1, NcompJ = dim */
152: void g1_pu(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g1[])
153: {
154: const PetscInt dim = spatialDim;
155: PetscInt d;
157: for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
158: }
160: /* -< \nabla\cdot v, p >
161: NcompI = dim, NcompJ = 1 */
162: void g2_up(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g2[])
163: {
164: const PetscInt dim = spatialDim;
165: PetscInt d;
167: for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */
168: }
170: /* < \nabla v, \nabla u + {\nabla u}^T >
171: This just gives \nabla u, give the perdiagonal for the transpose */
172: void g3_uu(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[])
173: {
174: const PetscInt dim = spatialDim;
175: const PetscInt Ncomp = spatialDim;
176: PetscInt compI, d;
178: for (compI = 0; compI < Ncomp; ++compI) {
179: for (d = 0; d < dim; ++d) {
180: g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
181: }
182: }
183: }
185: /*
186: In 3D we use exact solution:
188: u = x^2 + y^2
189: v = y^2 + z^2
190: w = x^2 + y^2 - 2(x+y)z
191: p = x + y + z - 3/2
192: f_x = f_y = f_z = 3
194: so that
196: -\Delta u + \nabla p + f = <-4, -4, -4> + <1, 1, 1> + <3, 3, 3> = 0
197: \nabla \cdot u = 2x + 2y - 2(x + y) = 0
198: */
199: void quadratic_u_3d(const PetscReal x[], PetscScalar *u, void *ctx)
200: {
201: u[0] = x[0]*x[0] + x[1]*x[1];
202: u[1] = x[1]*x[1] + x[2]*x[2];
203: u[2] = x[0]*x[0] + x[1]*x[1] - 2.0*(x[0] + x[1])*x[2];
204: }
206: void linear_p_3d(const PetscReal x[], PetscScalar *p, void *ctx)
207: {
208: *p = x[0] + x[1] + x[2] - 1.5;
209: }
213: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
214: {
215: const char *bcTypes[2] = {"neumann", "dirichlet"};
216: const char *runTypes[2] = {"full", "test"};
217: PetscInt bc, run;
221: options->debug = 0;
222: options->runType = RUN_FULL;
223: options->dim = 2;
224: options->interpolate = PETSC_FALSE;
225: options->simplex = PETSC_TRUE;
226: options->refinementLimit = 0.0;
227: options->bcType = DIRICHLET;
228: options->showInitial = PETSC_FALSE;
229: options->showSolution = PETSC_TRUE;
230: options->showError = PETSC_FALSE;
232: PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
233: PetscOptionsInt("-debug", "The debugging level", "ex62.c", options->debug, &options->debug, NULL);
234: run = options->runType;
235: PetscOptionsEList("-run_type", "The run type", "ex62.c", runTypes, 2, runTypes[options->runType], &run, NULL);
237: options->runType = (RunType) run;
239: PetscOptionsInt("-dim", "The topological mesh dimension", "ex62.c", options->dim, &options->dim, NULL);
240: spatialDim = options->dim;
241: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex62.c", options->interpolate, &options->interpolate, NULL);
242: PetscOptionsBool("-simplex", "Use simplices or tensor product cells", "ex62.c", options->simplex, &options->simplex, NULL);
243: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex62.c", options->refinementLimit, &options->refinementLimit, NULL);
244: PetscStrcpy(options->partitioner, "chaco");
245: PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, NULL);
246: bc = options->bcType;
247: PetscOptionsEList("-bc_type","Type of boundary condition","ex62.c",bcTypes,2,bcTypes[options->bcType],&bc,NULL);
249: options->bcType = (BCType) bc;
251: PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex62.c", options->showInitial, &options->showInitial, NULL);
252: PetscOptionsBool("-show_solution", "Output the solution for verification", "ex62.c", options->showSolution, &options->showSolution, NULL);
253: PetscOptionsBool("-show_error", "Output the error for verification", "ex62.c", options->showError, &options->showError, NULL);
254: PetscOptionsEnd();
256: PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
257: return(0);
258: }
262: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
263: {
264: Vec lv;
265: PetscInt p;
266: PetscMPIInt rank, numProcs;
270: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
271: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
272: DMGetLocalVector(dm, &lv);
273: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
274: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
275: PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
276: for (p = 0; p < numProcs; ++p) {
277: if (p == rank) {VecView(lv, PETSC_VIEWER_STDOUT_SELF);}
278: PetscBarrier((PetscObject) dm);
279: }
280: DMRestoreLocalVector(dm, &lv);
281: return(0);
282: }
286: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
287: {
288: DMLabel label;
289: PetscInt dim = user->dim;
290: PetscBool interpolate = user->interpolate;
291: PetscReal refinementLimit = user->refinementLimit;
292: const char *partitioner = user->partitioner;
293: const PetscInt cells[3] = {3, 3, 3};
297: PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
298: if (user->simplex) {DMPlexCreateBoxMesh(comm, dim, interpolate, dm);}
299: else {DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);}
300: DMPlexGetLabel(*dm, "marker", &label);
301: if (label) {DMPlexLabelComplete(*dm, label);}
302: {
303: DM refinedMesh = NULL;
304: DM distributedMesh = NULL;
306: /* Refine mesh using a volume constraint */
307: DMPlexSetRefinementLimit(*dm, refinementLimit);
308: if (user->simplex) {DMRefine(*dm, comm, &refinedMesh);}
309: if (refinedMesh) {
310: DMDestroy(dm);
311: *dm = refinedMesh;
312: }
313: /* Distribute mesh over processes */
314: DMPlexDistribute(*dm, partitioner, 0, NULL, &distributedMesh);
315: if (distributedMesh) {
316: DMDestroy(dm);
317: *dm = distributedMesh;
318: }
319: }
320: DMSetFromOptions(*dm);
321: DMViewFromOptions(*dm, NULL, "-dm_view");
322: PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
323: return(0);
324: }
328: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
329: {
330: PetscDS prob;
334: DMGetDS(dm, &prob);
335: PetscDSSetResidual(prob, 0, f0_u, f1_u);
336: PetscDSSetResidual(prob, 1, f0_p, f1_p);
337: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
338: PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_up, NULL);
339: PetscDSSetJacobian(prob, 1, 0, NULL, g1_pu, NULL, NULL);
340: switch (user->dim) {
341: case 2:
342: user->exactFuncs[0] = quadratic_u_2d;
343: user->exactFuncs[1] = linear_p_2d;
344: break;
345: case 3:
346: user->exactFuncs[0] = quadratic_u_3d;
347: user->exactFuncs[1] = linear_p_3d;
348: break;
349: default:
350: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
351: }
352: return(0);
353: }
357: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
358: {
359: DM cdm = dm;
360: const PetscInt dim = user->dim;
361: const PetscInt id = 1;
362: PetscFE fe[2];
363: PetscSpace P;
364: PetscDS prob;
365: PetscInt order;
369: /* Create finite element */
370: PetscFECreateDefault(dm, dim, dim, user->simplex, "vel_", -1, &fe[0]);
371: PetscObjectSetName((PetscObject) fe[0], "velocity");
372: PetscFEGetBasisSpace(fe[0], &P);
373: PetscSpaceGetOrder(P, &order);
374: PetscFECreateDefault(dm, dim, 1, user->simplex, "pres_", order, &fe[1]);
375: PetscObjectSetName((PetscObject) fe[1], "pressure");
376: /* Set discretization and boundary conditions for each mesh */
377: while (cdm) {
378: DMGetDS(cdm, &prob);
379: PetscDSSetDiscretization(prob, 0, (PetscObject) fe[0]);
380: PetscDSSetDiscretization(prob, 1, (PetscObject) fe[1]);
381: SetupProblem(cdm, user);
382: DMPlexAddBoundary(cdm, user->bcType == DIRICHLET, "wall", user->bcType == NEUMANN ? "boundary" : "marker", 0, user->exactFuncs[0], 1, &id, user);
383: DMPlexGetCoarseDM(cdm, &cdm);
384: }
385: PetscFEDestroy(&fe[0]);
386: PetscFEDestroy(&fe[1]);
387: return(0);
388: }
392: PetscErrorCode CreatePressureNullSpace(DM dm, AppCtx *user, Vec *v, MatNullSpace *nullSpace)
393: {
394: Vec vec;
395: void (*funcs[2])(const PetscReal x[], PetscScalar *u, void* ctx) = {zero_vector, constant_p};
399: DMGetGlobalVector(dm, &vec);
400: DMPlexProjectFunction(dm, funcs, NULL, INSERT_ALL_VALUES, vec);
401: VecNormalize(vec, NULL);
402: if (user->debug) {
403: PetscPrintf(PetscObjectComm((PetscObject)dm), "Pressure Null Space\n");
404: VecView(vec, PETSC_VIEWER_STDOUT_WORLD);
405: }
406: MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace);
407: if (v) {
408: DMCreateGlobalVector(dm, v);
409: VecCopy(vec, *v);
410: }
411: DMRestoreGlobalVector(dm, &vec);
412: /* New style for field null spaces */
413: {
414: PetscObject pressure;
415: MatNullSpace nullSpacePres;
417: DMGetField(dm, 1, &pressure);
418: MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullSpacePres);
419: PetscObjectCompose(pressure, "nullspace", (PetscObject) nullSpacePres);
420: MatNullSpaceDestroy(&nullSpacePres);
421: }
422: return(0);
423: }
427: int main(int argc, char **argv)
428: {
429: SNES snes; /* nonlinear solver */
430: DM dm; /* problem definition */
431: Vec u,r; /* solution, residual vectors */
432: Mat A,J; /* Jacobian matrix */
433: MatNullSpace nullSpace; /* May be necessary for pressure */
434: AppCtx user; /* user-defined work context */
435: PetscInt its; /* iterations for convergence */
436: PetscReal error = 0.0; /* L_2 error in the solution */
437: PetscReal ferrors[2];
440: PetscInitialize(&argc, &argv, NULL, help);
441: ProcessOptions(PETSC_COMM_WORLD, &user);
442: SNESCreate(PETSC_COMM_WORLD, &snes);
443: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
444: SNESSetDM(snes, dm);
445: DMSetApplicationContext(dm, &user);
447: PetscMalloc(2 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs);
448: SetupDiscretization(dm, &user);
449: DMPlexCreateClosureIndex(dm, NULL);
451: DMCreateGlobalVector(dm, &u);
452: VecDuplicate(u, &r);
454: DMSetMatType(dm,MATAIJ);
455: DMCreateMatrix(dm, &J);
456: A = J;
457: CreatePressureNullSpace(dm, &user, NULL, &nullSpace);
458: MatSetNullSpace(J, nullSpace);
459: if (A != J) {
460: MatSetNullSpace(A, nullSpace);
461: }
463: DMSNESSetFunctionLocal(dm, (PetscErrorCode (*)(DM,Vec,Vec,void*))DMPlexSNESComputeResidualFEM,&user);
464: DMSNESSetJacobianLocal(dm, (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*))DMPlexSNESComputeJacobianFEM,&user);
465: SNESSetJacobian(snes, A, J, NULL, NULL);
467: SNESSetFromOptions(snes);
469: DMPlexProjectFunction(dm, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
470: if (user.showInitial) {DMVecViewLocal(dm, u, PETSC_VIEWER_STDOUT_SELF);}
471: if (user.runType == RUN_FULL) {
472: void (*initialGuess[2])(const PetscReal x[], PetscScalar *u, void* ctx) = {zero_vector, zero_scalar};
474: DMPlexProjectFunction(dm, initialGuess, NULL, INSERT_VALUES, u);
475: if (user.debug) {
476: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
477: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
478: }
479: SNESSolve(snes, NULL, u);
480: SNESGetIterationNumber(snes, &its);
481: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
482: DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
483: DMPlexComputeL2FieldDiff(dm, user.exactFuncs, NULL, u, ferrors);
484: PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %.3g [%.3g, %.3g]\n", error, ferrors[0], ferrors[1]);
485: if (user.showError) {
486: Vec r;
487: DMGetGlobalVector(dm, &r);
488: DMPlexProjectFunction(dm, user.exactFuncs, NULL, INSERT_ALL_VALUES, r);
489: VecAXPY(r, -1.0, u);
490: PetscPrintf(PETSC_COMM_WORLD, "Solution Error\n");
491: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
492: DMRestoreGlobalVector(dm, &r);
493: }
494: if (user.showSolution) {
495: PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
496: VecChop(u, 3.0e-9);
497: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
498: }
499: } else {
500: PetscReal res = 0.0;
502: /* Check discretization error */
503: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
504: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
505: DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
506: if (error >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
507: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
508: /* Check residual */
509: SNESComputeFunction(snes, u, r);
510: PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
511: VecChop(r, 1.0e-10);
512: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
513: VecNorm(r, NORM_2, &res);
514: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
515: /* Check Jacobian */
516: {
517: Vec b;
518: PetscBool isNull;
520: SNESComputeJacobian(snes, u, A, A);
521: MatNullSpaceTest(nullSpace, J, &isNull);
522: //if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
523: VecDuplicate(u, &b);
524: VecSet(r, 0.0);
525: SNESComputeFunction(snes, r, b);
526: MatMult(A, u, r);
527: VecAXPY(r, 1.0, b);
528: VecDestroy(&b);
529: PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
530: VecChop(r, 1.0e-10);
531: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
532: VecNorm(r, NORM_2, &res);
533: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
534: }
535: }
536: VecViewFromOptions(u, NULL, "-sol_vec_view");
538: MatNullSpaceDestroy(&nullSpace);
539: if (A != J) {MatDestroy(&A);}
540: MatDestroy(&J);
541: VecDestroy(&u);
542: VecDestroy(&r);
543: SNESDestroy(&snes);
544: DMDestroy(&dm);
545: PetscFree(user.exactFuncs);
546: PetscFinalize();
547: return 0;
548: }