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