Actual source code: ex12.c
petsc-3.5.4 2015-05-23
1: static char help[] = "Poisson Problem in 2d and 3d with simplicial finite elements.\n\
2: We solve the Poisson problem in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
4: This example supports discretized auxiliary fields (conductivity) as well as\n\
5: multilevel nonlienar solvers.\n\n\n";
7: #include <petscdmplex.h>
8: #include <petscsnes.h>
9: #include <petscds.h>
10: #include <petscviewerhdf5.h>
12: PetscInt spatialDim = 0;
14: typedef enum {NEUMANN, DIRICHLET, NONE} BCType;
15: typedef enum {RUN_FULL, RUN_TEST, RUN_PERF} RunType;
16: typedef enum {COEFF_NONE, COEFF_ANALYTIC, COEFF_FIELD, COEFF_NONLINEAR} CoeffType;
18: typedef struct {
19: PetscInt debug; /* The debugging level */
20: RunType runType; /* Whether to run tests, or solve the full problem */
21: PetscBool jacobianMF; /* Whether to calculate the Jacobian action on the fly */
22: PetscLogEvent createMeshEvent;
23: PetscBool showInitial, showSolution, restart, check;
24: PetscViewer checkpoint;
25: /* Domain and mesh definition */
26: PetscInt dim; /* The topological mesh dimension */
27: char filename[2048]; /* The optional ExodusII file */
28: PetscBool interpolate; /* Generate intermediate mesh elements */
29: PetscReal refinementLimit; /* The largest allowable cell volume */
30: PetscBool refinementUniform; /* Uniformly refine the mesh */
31: PetscInt refinementRounds; /* The number of uniform refinements */
32: char partitioner[2048]; /* The graph partitioner */
33: /* Problem definition */
34: BCType bcType;
35: CoeffType variableCoefficient;
36: void (**exactFuncs)(const PetscReal x[], PetscScalar *u, void *ctx);
37: } AppCtx;
39: void zero(const PetscReal coords[], PetscScalar *u, void *ctx)
40: {
41: *u = 0.0;
42: }
44: /*
45: In 2D for Dirichlet conditions, we use exact solution:
47: u = x^2 + y^2
48: f = 4
50: so that
52: -\Delta u + f = -4 + 4 = 0
54: For Neumann conditions, we have
56: -\nabla u \cdot -\hat y |_{y=0} = (2y)|_{y=0} = 0 (bottom)
57: -\nabla u \cdot \hat y |_{y=1} = -(2y)|_{y=1} = -2 (top)
58: -\nabla u \cdot -\hat x |_{x=0} = (2x)|_{x=0} = 0 (left)
59: -\nabla u \cdot \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)
61: Which we can express as
63: \nabla u \cdot \hat n|_\Gamma = {2 x, 2 y} \cdot \hat n = 2 (x + y)
64: */
65: void quadratic_u_2d(const PetscReal x[], PetscScalar *u, void *ctx)
66: {
67: *u = x[0]*x[0] + x[1]*x[1];
68: }
70: 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[])
71: {
72: f0[0] = 4.0;
73: }
75: void f0_bd_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[], const PetscReal n[], PetscScalar f0[])
76: {
77: PetscInt d;
78: for (d = 0, f0[0] = 0.0; d < spatialDim; ++d) f0[0] += -n[d]*2.0*x[d];
79: }
81: void f0_bd_zero(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[], const PetscReal n[], PetscScalar f0[])
82: {
83: f0[0] = 0.0;
84: }
86: void f1_bd_zero(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[], const PetscReal n[], PetscScalar f1[])
87: {
88: PetscInt comp;
89: for (comp = 0; comp < spatialDim; ++comp) f1[comp] = 0.0;
90: }
92: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
93: 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[])
94: {
95: PetscInt d;
97: for (d = 0; d < spatialDim; ++d) f1[d] = u_x[d];
98: }
100: /* < \nabla v, \nabla u + {\nabla u}^T >
101: This just gives \nabla u, give the perdiagonal for the transpose */
102: 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[])
103: {
104: PetscInt d;
106: for (d = 0; d < spatialDim; ++d) g3[d*spatialDim+d] = 1.0;
107: }
109: /*
110: In 2D for Dirichlet conditions with a variable coefficient, we use exact solution:
112: u = x^2 + y^2
113: f = 6 (x + y)
114: nu = (x + y)
116: so that
118: -\div \nu \grad u + f = -6 (x + y) + 6 (x + y) = 0
119: */
120: void nu_2d(const PetscReal x[], PetscScalar *u, void *ctx)
121: {
122: *u = x[0] + x[1];
123: }
125: void f0_analytic_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[])
126: {
127: f0[0] = 6.0*(x[0] + x[1]);
128: }
130: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
131: void f1_analytic_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[])
132: {
133: PetscInt d;
135: for (d = 0; d < spatialDim; ++d) f1[d] = (x[0] + x[1])*u_x[d];
136: }
137: void f1_field_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[])
138: {
139: PetscInt d;
141: for (d = 0; d < spatialDim; ++d) f1[d] = a[0]*u_x[d];
142: }
144: /* < \nabla v, \nabla u + {\nabla u}^T >
145: This just gives \nabla u, give the perdiagonal for the transpose */
146: void g3_analytic_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[])
147: {
148: PetscInt d;
150: for (d = 0; d < spatialDim; ++d) g3[d*spatialDim+d] = x[0] + x[1];
151: }
152: void g3_field_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[])
153: {
154: PetscInt d;
156: for (d = 0; d < spatialDim; ++d) g3[d*spatialDim+d] = a[0];
157: }
159: /*
160: In 2D for Dirichlet conditions with a nonlinear coefficient (p-Laplacian with p = 4), we use exact solution:
162: u = x^2 + y^2
163: f = 16 (x^2 + y^2)
164: nu = 1/2 |grad u|^2
166: so that
168: -\div \nu \grad u + f = -16 (x^2 + y^2) + 16 (x^2 + y^2) = 0
169: */
170: void f0_analytic_nonlinear_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[])
171: {
172: f0[0] = 16.0*(x[0]*x[0] + x[1]*x[1]);
173: }
175: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
176: void f1_analytic_nonlinear_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[])
177: {
178: PetscScalar nu = 0.0;
179: PetscInt d;
180: for (d = 0; d < spatialDim; ++d) nu += u_x[d]*u_x[d];
181: for (d = 0; d < spatialDim; ++d) f1[d] = 0.5*nu*u_x[d];
182: }
184: /*
185: grad (u + eps w) - grad u = eps grad w
187: 1/2 |grad (u + eps w)|^2 grad (u + eps w) - 1/2 |grad u|^2 grad u
188: = 1/2 (|grad u|^2 + 2 eps <grad u,grad w>) (grad u + eps grad w) - 1/2 |grad u|^2 grad u
189: = 1/2 (eps |grad u|^2 grad w + 2 eps <grad u,grad w> grad u)
190: = eps (1/2 |grad u|^2 grad w + grad u <grad u,grad w>)
191: */
192: void g3_analytic_nonlinear_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[])
193: {
194: PetscScalar nu = 0.0;
195: PetscInt d, e;
196: for (d = 0; d < spatialDim; ++d) nu += u_x[d]*u_x[d];
197: for (d = 0; d < spatialDim; ++d) {
198: g3[d*spatialDim+d] = 0.5*nu;
199: for (e = 0; e < spatialDim; ++e) {
200: g3[d*spatialDim+e] += u_x[d]*u_x[e];
201: }
202: }
203: }
205: /*
206: In 3D for Dirichlet conditions we use exact solution:
208: u = x^2 + y^2 + z^2
209: f = 6
211: so that
213: -\Delta u + f = -6 + 6 = 0
215: For Neumann conditions, we have
217: -\nabla u \cdot -\hat z |_{z=0} = (2z)|_{z=0} = 0 (bottom)
218: -\nabla u \cdot \hat z |_{z=1} = -(2z)|_{z=1} = -2 (top)
219: -\nabla u \cdot -\hat y |_{y=0} = (2y)|_{y=0} = 0 (front)
220: -\nabla u \cdot \hat y |_{y=1} = -(2y)|_{y=1} = -2 (back)
221: -\nabla u \cdot -\hat x |_{x=0} = (2x)|_{x=0} = 0 (left)
222: -\nabla u \cdot \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)
224: Which we can express as
226: \nabla u \cdot \hat n|_\Gamma = {2 x, 2 y, 2z} \cdot \hat n = 2 (x + y + z)
227: */
228: void quadratic_u_3d(const PetscReal x[], PetscScalar *u, void *ctx)
229: {
230: *u = x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
231: }
235: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
236: {
237: const char *bcTypes[3] = {"neumann", "dirichlet", "none"};
238: const char *runTypes[3] = {"full", "test", "perf"};
239: const char *coeffTypes[4] = {"none", "analytic", "field", "nonlinear"};
240: PetscInt bc, run, coeff;
241: PetscBool flg;
245: options->debug = 0;
246: options->runType = RUN_FULL;
247: options->dim = 2;
248: options->filename[0] = '\0';
249: options->interpolate = PETSC_FALSE;
250: options->refinementLimit = 0.0;
251: options->refinementUniform = PETSC_FALSE;
252: options->refinementRounds = 1;
253: options->bcType = DIRICHLET;
254: options->variableCoefficient = COEFF_NONE;
255: options->jacobianMF = PETSC_FALSE;
256: options->showInitial = PETSC_FALSE;
257: options->showSolution = PETSC_FALSE;
258: options->restart = PETSC_FALSE;
259: options->check = PETSC_FALSE;
260: options->checkpoint = NULL;
262: PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
263: PetscOptionsInt("-debug", "The debugging level", "ex12.c", options->debug, &options->debug, NULL);
264: run = options->runType;
265: PetscOptionsEList("-run_type", "The run type", "ex12.c", runTypes, 3, runTypes[options->runType], &run, NULL);
267: options->runType = (RunType) run;
269: PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);
270: spatialDim = options->dim;
271: PetscOptionsString("-f", "Exodus.II filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);
272: #if !defined(PETSC_HAVE_EXODUSII)
273: if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "This option requires ExodusII support. Reconfigure using --download-exodusii");
274: #endif
275: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);
276: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, NULL);
277: PetscOptionsBool("-refinement_uniform", "Uniformly refine the mesh", "ex52.c", options->refinementUniform, &options->refinementUniform, NULL);
278: PetscOptionsInt("-refinement_rounds", "The number of uniform refinements", "ex52.c", options->refinementRounds, &options->refinementRounds, NULL);
279: PetscStrcpy(options->partitioner, "chaco");
280: PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, NULL);
281: bc = options->bcType;
282: PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,3,bcTypes[options->bcType],&bc,NULL);
283: options->bcType = (BCType) bc;
284: coeff = options->variableCoefficient;
285: PetscOptionsEList("-variable_coefficient","Type of variable coefficent","ex12.c",coeffTypes,4,coeffTypes[options->variableCoefficient],&coeff,NULL);
286: options->variableCoefficient = (CoeffType) coeff;
288: PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex12.c", options->jacobianMF, &options->jacobianMF, NULL);
289: PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex12.c", options->showInitial, &options->showInitial, NULL);
290: PetscOptionsBool("-show_solution", "Output the solution for verification", "ex12.c", options->showSolution, &options->showSolution, NULL);
291: PetscOptionsBool("-restart", "Read in the mesh and solution from a file", "ex12.c", options->restart, &options->restart, NULL);
292: PetscOptionsBool("-check", "Compare with default integration routines", "ex12.c", options->check, &options->check, NULL);
293: PetscOptionsEnd();
295: PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
297: if (options->restart) {
298: PetscViewerCreate(comm, &options->checkpoint);
299: PetscViewerSetType(options->checkpoint, PETSCVIEWERHDF5);
300: PetscViewerFileSetMode(options->checkpoint, FILE_MODE_READ);
301: PetscViewerFileSetName(options->checkpoint, options->filename);
302: }
303: return(0);
304: }
308: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
309: {
310: PetscInt dim = user->dim;
311: const char *filename = user->filename;
312: PetscBool interpolate = user->interpolate;
313: PetscReal refinementLimit = user->refinementLimit;
314: PetscBool refinementUniform = user->refinementUniform;
315: PetscInt refinementRounds = user->refinementRounds;
316: const char *partitioner = user->partitioner;
317: size_t len;
321: PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
322: PetscStrlen(filename, &len);
323: if (!len) {
324: DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
325: PetscObjectSetName((PetscObject) *dm, "Mesh");
326: } else if (user->checkpoint) {
327: DMCreate(comm, dm);
328: DMSetType(*dm, DMPLEX);
329: DMLoad(*dm, user->checkpoint);
330: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
331: } else {
332: PetscMPIInt rank;
334: MPI_Comm_rank(comm, &rank);
335: DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);
336: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
337: /* Must have boundary marker for Dirichlet conditions */
338: }
339: {
340: DM refinedMesh = NULL;
341: DM distributedMesh = NULL;
343: /* Refine mesh using a volume constraint */
344: DMPlexSetRefinementLimit(*dm, refinementLimit);
345: DMRefine(*dm, comm, &refinedMesh);
346: if (refinedMesh) {
347: const char *name;
349: PetscObjectGetName((PetscObject) *dm, &name);
350: PetscObjectSetName((PetscObject) refinedMesh, name);
351: DMDestroy(dm);
352: *dm = refinedMesh;
353: }
354: /* Distribute mesh over processes */
355: DMPlexDistribute(*dm, partitioner, 0, NULL, &distributedMesh);
356: if (distributedMesh) {
357: DMDestroy(dm);
358: *dm = distributedMesh;
359: }
360: /* Use regular refinement in parallel */
361: if (refinementUniform) {
362: PetscInt r;
364: DMPlexSetRefinementUniform(*dm, refinementUniform);
365: for (r = 0; r < refinementRounds; ++r) {
366: DMRefine(*dm, comm, &refinedMesh);
367: if (refinedMesh) {
368: DMDestroy(dm);
369: *dm = refinedMesh;
370: }
371: }
372: }
373: }
374: DMSetFromOptions(*dm);
375: if (user->bcType == NEUMANN) {
376: DMLabel label;
378: DMPlexCreateLabel(*dm, "boundary");
379: DMPlexGetLabel(*dm, "boundary", &label);
380: DMPlexMarkBoundaryFaces(*dm, label);
381: }
382: DMViewFromOptions(*dm, NULL, "-dm_view");
383: PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
384: return(0);
385: }
389: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
390: {
391: PetscDS prob;
395: DMGetDS(dm, &prob);
396: switch (user->variableCoefficient) {
397: case COEFF_NONE:
398: PetscDSSetResidual(prob, 0, f0_u, f1_u);
399: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
400: break;
401: case COEFF_ANALYTIC:
402: PetscDSSetResidual(prob, 0, f0_analytic_u, f1_analytic_u);
403: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_uu);
404: break;
405: case COEFF_FIELD:
406: PetscDSSetResidual(prob, 0, f0_analytic_u, f1_field_u);
407: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);
408: break;
409: case COEFF_NONLINEAR:
410: PetscDSSetResidual(prob, 0, f0_analytic_nonlinear_u, f1_analytic_nonlinear_u);
411: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_nonlinear_uu);
412: break;
413: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid variable coefficient type %d", user->variableCoefficient);
414: }
415: switch (user->dim) {
416: case 2:
417: user->exactFuncs[0] = quadratic_u_2d;
418: if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
419: break;
420: case 3:
421: user->exactFuncs[0] = quadratic_u_3d;
422: if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
423: break;
424: default:
425: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
426: }
427: return(0);
428: }
432: PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
433: {
434: void (*matFuncs[1])(const PetscReal x[], PetscScalar *u, void *ctx) = {nu_2d};
435: Vec nu;
439: DMCreateLocalVector(dmAux, &nu);
440: DMPlexProjectFunctionLocal(dmAux, matFuncs, NULL, INSERT_ALL_VALUES, nu);
441: PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
442: VecDestroy(&nu);
443: return(0);
444: }
448: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
449: {
450: DM cdm = dm;
451: const PetscInt dim = user->dim;
452: const PetscInt id = 1;
453: PetscFE feAux = NULL;
454: PetscFE feBd = NULL;
455: PetscFE feCh = NULL;
456: PetscFE fe;
457: PetscDS prob;
461: /* Create finite element */
462: PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, NULL, -1, &fe);
463: PetscObjectSetName((PetscObject) fe, "potential");
464: if (user->bcType == NEUMANN) {
465: PetscFECreateDefault(dm, dim-1, 1, PETSC_TRUE, "bd_", -1, &feBd);
466: PetscObjectSetName((PetscObject) feBd, "potential");
467: }
468: if (user->variableCoefficient == COEFF_FIELD) {
469: PetscQuadrature q;
471: PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, "mat_", -1, &feAux);
472: PetscFEGetQuadrature(fe, &q);
473: PetscFESetQuadrature(feAux, q);
474: }
475: if (user->check) {PetscFECreateDefault(dm, dim, 1, PETSC_TRUE, "ch_", -1, &feCh);}
476: /* Set discretization and boundary conditions for each mesh */
477: while (cdm) {
478: DMGetDS(cdm, &prob);
479: PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
480: PetscDSSetBdDiscretization(prob, 0, (PetscObject) feBd);
481: if (feAux) {
482: DM dmAux;
483: PetscDS probAux;
485: DMClone(cdm, &dmAux);
486: DMPlexCopyCoordinates(cdm, dmAux);
487: DMGetDS(dmAux, &probAux);
488: PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
489: PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
490: SetupMaterial(cdm, dmAux, user);
491: DMDestroy(&dmAux);
492: }
493: if (feCh) {
494: DM dmCh;
495: PetscDS probCh;
497: DMClone(cdm, &dmCh);
498: DMPlexCopyCoordinates(cdm, dmCh);
499: DMGetDS(dmCh, &probCh);
500: PetscDSSetDiscretization(probCh, 0, (PetscObject) feCh);
501: PetscObjectCompose((PetscObject) dm, "dmCh", (PetscObject) dmCh);
502: DMDestroy(&dmCh);
503: }
504: SetupProblem(cdm, user);
505: DMPlexAddBoundary(cdm, user->bcType == DIRICHLET, "wall", user->bcType == NEUMANN ? "boundary" : "marker", 0, user->exactFuncs[0], 1, &id, user);
506: DMPlexGetCoarseDM(cdm, &cdm);
507: }
508: PetscFEDestroy(&fe);
509: PetscFEDestroy(&feBd);
510: PetscFEDestroy(&feAux);
511: PetscFEDestroy(&feCh);
512: return(0);
513: }
517: int main(int argc, char **argv)
518: {
519: DM dm; /* Problem specification */
520: SNES snes; /* nonlinear solver */
521: Vec u,r; /* solution, residual vectors */
522: Mat A,J; /* Jacobian matrix */
523: MatNullSpace nullSpace; /* May be necessary for Neumann conditions */
524: AppCtx user; /* user-defined work context */
525: JacActionCtx userJ; /* context for Jacobian MF action */
526: PetscInt its; /* iterations for convergence */
527: PetscReal error = 0.0; /* L_2 error in the solution */
530: PetscInitialize(&argc, &argv, NULL, help);
531: ProcessOptions(PETSC_COMM_WORLD, &user);
532: SNESCreate(PETSC_COMM_WORLD, &snes);
533: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
534: SNESSetDM(snes, dm);
535: DMSetApplicationContext(dm, &user);
537: PetscMalloc(1 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs);
538: SetupDiscretization(dm, &user);
540: DMCreateGlobalVector(dm, &u);
541: PetscObjectSetName((PetscObject) u, "potential");
542: VecDuplicate(u, &r);
544: DMSetMatType(dm,MATAIJ);
545: DMCreateMatrix(dm, &J);
546: if (user.jacobianMF) {
547: PetscInt M, m, N, n;
549: MatGetSize(J, &M, &N);
550: MatGetLocalSize(J, &m, &n);
551: MatCreate(PETSC_COMM_WORLD, &A);
552: MatSetSizes(A, m, n, M, N);
553: MatSetType(A, MATSHELL);
554: MatSetUp(A);
555: #if 0
556: MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);
557: #endif
559: userJ.dm = dm;
560: userJ.J = J;
561: userJ.user = &user;
563: DMCreateLocalVector(dm, &userJ.u);
564: DMPlexProjectFunctionLocal(dm, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);
565: MatShellSetContext(A, &userJ);
566: } else {
567: A = J;
568: }
569: if (user.bcType == NEUMANN) {
570: MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);
571: MatSetNullSpace(J, nullSpace);
572: if (A != J) {
573: MatSetNullSpace(A, nullSpace);
574: }
575: }
577: DMSNESSetFunctionLocal(dm, (PetscErrorCode (*)(DM,Vec,Vec,void*)) DMPlexSNESComputeResidualFEM, &user);
578: DMSNESSetJacobianLocal(dm, (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*)) DMPlexSNESComputeJacobianFEM, &user);
579: SNESSetJacobian(snes, A, J, NULL, NULL);
581: SNESSetFromOptions(snes);
583: DMPlexProjectFunction(dm, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
584: if (user.checkpoint) {
585: #if defined(PETSC_HAVE_HDF5)
586: PetscViewerHDF5PushGroup(user.checkpoint, "/fields");
587: VecLoad(u, user.checkpoint);
588: PetscViewerHDF5PopGroup(user.checkpoint);
589: #endif
590: }
591: PetscViewerDestroy(&user.checkpoint);
592: if (user.showInitial) {
593: Vec lv;
594: DMGetLocalVector(dm, &lv);
595: DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);
596: DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);
597: DMPrintLocalVec(dm, "Local function", 1.0e-10, lv);
598: DMRestoreLocalVector(dm, &lv);
599: }
600: if (user.runType == RUN_FULL) {
601: void (*initialGuess[1])(const PetscReal x[], PetscScalar *, void *ctx) = {zero};
603: DMPlexProjectFunction(dm, initialGuess, NULL, INSERT_VALUES, u);
604: if (user.debug) {
605: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
606: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
607: }
608: SNESSolve(snes, NULL, u);
609: SNESGetIterationNumber(snes, &its);
610: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
611: DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
612: if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
613: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
614: if (user.showSolution) {
615: PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
616: VecChop(u, 3.0e-9);
617: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
618: }
619: } else if (user.runType == RUN_PERF) {
620: PetscReal res = 0.0;
622: SNESComputeFunction(snes, u, r);
623: PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
624: VecChop(r, 1.0e-10);
625: VecNorm(r, NORM_2, &res);
626: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
627: } else {
628: PetscReal res = 0.0;
630: /* Check discretization error */
631: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
632: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
633: DMPlexComputeL2Diff(dm, user.exactFuncs, NULL, u, &error);
634: if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
635: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
636: /* Check residual */
637: SNESComputeFunction(snes, u, r);
638: PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
639: VecChop(r, 1.0e-10);
640: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
641: VecNorm(r, NORM_2, &res);
642: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
643: /* Check Jacobian */
644: {
645: Vec b;
647: SNESComputeJacobian(snes, u, A, A);
648: VecDuplicate(u, &b);
649: VecSet(r, 0.0);
650: SNESComputeFunction(snes, r, b);
651: MatMult(A, u, r);
652: VecAXPY(r, 1.0, b);
653: VecDestroy(&b);
654: PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
655: VecChop(r, 1.0e-10);
656: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
657: VecNorm(r, NORM_2, &res);
658: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
659: }
660: }
661: VecViewFromOptions(u, NULL, "-vec_view");
663: if (user.bcType == NEUMANN) {MatNullSpaceDestroy(&nullSpace);}
664: if (user.jacobianMF) {VecDestroy(&userJ.u);}
665: if (A != J) {MatDestroy(&A);}
666: MatDestroy(&J);
667: VecDestroy(&u);
668: VecDestroy(&r);
669: SNESDestroy(&snes);
670: DMDestroy(&dm);
671: PetscFree(user.exactFuncs);
672: PetscFinalize();
673: return 0;
674: }