Actual source code: ex12.c
petsc-3.7.5 2017-01-01
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_EXACT, 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: /* Domain and mesh definition */
23: PetscInt dim; /* The topological mesh dimension */
24: char filename[2048]; /* The optional ExodusII file */
25: PetscBool interpolate; /* Generate intermediate mesh elements */
26: PetscReal refinementLimit; /* The largest allowable cell volume */
27: PetscBool viewHierarchy; /* Whether to view the hierarchy */
28: PetscBool simplex; /* Simplicial mesh */
29: /* Problem definition */
30: BCType bcType;
31: CoeffType variableCoefficient;
32: PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
33: /* Solver */
34: PC pcmg; /* This is needed for error monitoring */
35: } AppCtx;
37: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
38: {
39: u[0] = 0.0;
40: return 0;
41: }
43: /*
44: In 2D for Dirichlet conditions, we use exact solution:
46: u = x^2 + y^2
47: f = 4
49: so that
51: -\Delta u + f = -4 + 4 = 0
53: For Neumann conditions, we have
55: -\nabla u \cdot -\hat y |_{y=0} = (2y)|_{y=0} = 0 (bottom)
56: -\nabla u \cdot \hat y |_{y=1} = -(2y)|_{y=1} = -2 (top)
57: -\nabla u \cdot -\hat x |_{x=0} = (2x)|_{x=0} = 0 (left)
58: -\nabla u \cdot \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)
60: Which we can express as
62: \nabla u \cdot \hat n|_\Gamma = {2 x, 2 y} \cdot \hat n = 2 (x + y)
63: */
64: static PetscErrorCode quadratic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
65: {
66: *u = x[0]*x[0] + x[1]*x[1];
67: return 0;
68: }
70: void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
71: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
72: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
73: PetscReal t, const PetscReal x[], PetscScalar f0[])
74: {
75: f0[0] = 4.0;
76: }
78: void f0_bd_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
79: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
80: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
81: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
82: {
83: PetscInt d;
84: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += -n[d]*2.0*x[d];
85: }
87: void f0_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux,
88: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
89: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
90: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])
91: {
92: f0[0] = 0.0;
93: }
95: void f1_bd_zero(PetscInt dim, PetscInt Nf, PetscInt NfAux,
96: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
97: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
98: PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f1[])
99: {
100: PetscInt comp;
101: for (comp = 0; comp < dim; ++comp) f1[comp] = 0.0;
102: }
104: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
105: void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
106: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
107: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
108: PetscReal t, const PetscReal x[], PetscScalar f1[])
109: {
110: PetscInt d;
111: for (d = 0; d < dim; ++d) f1[d] = u_x[d];
112: }
114: /* < \nabla v, \nabla u + {\nabla u}^T >
115: This just gives \nabla u, give the perdiagonal for the transpose */
116: void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
117: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
118: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
119: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
120: {
121: PetscInt d;
122: for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0;
123: }
125: /*
126: In 2D for Dirichlet conditions with a variable coefficient, we use exact solution:
128: u = x^2 + y^2
129: f = 6 (x + y)
130: nu = (x + y)
132: so that
134: -\div \nu \grad u + f = -6 (x + y) + 6 (x + y) = 0
135: */
136: static PetscErrorCode nu_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
137: {
138: *u = x[0] + x[1];
139: return 0;
140: }
142: void f0_analytic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
143: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
144: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
145: PetscReal t, const PetscReal x[], PetscScalar f0[])
146: {
147: f0[0] = 6.0*(x[0] + x[1]);
148: }
150: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
151: void f1_analytic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
152: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
153: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
154: PetscReal t, const PetscReal x[], PetscScalar f1[])
155: {
156: PetscInt d;
157: for (d = 0; d < dim; ++d) f1[d] = (x[0] + x[1])*u_x[d];
158: }
160: void f1_field_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
161: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
162: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
163: PetscReal t, const PetscReal x[], PetscScalar f1[])
164: {
165: PetscInt d;
166: for (d = 0; d < dim; ++d) f1[d] = a[0]*u_x[d];
167: }
169: /* < \nabla v, \nabla u + {\nabla u}^T >
170: This just gives \nabla u, give the perdiagonal for the transpose */
171: void g3_analytic_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
172: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
173: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
174: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
175: {
176: PetscInt d;
177: for (d = 0; d < dim; ++d) g3[d*dim+d] = x[0] + x[1];
178: }
180: void g3_field_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
181: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
182: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
183: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
184: {
185: PetscInt d;
186: for (d = 0; d < dim; ++d) g3[d*dim+d] = a[0];
187: }
189: /*
190: In 2D for Dirichlet conditions with a nonlinear coefficient (p-Laplacian with p = 4), we use exact solution:
192: u = x^2 + y^2
193: f = 16 (x^2 + y^2)
194: nu = 1/2 |grad u|^2
196: so that
198: -\div \nu \grad u + f = -16 (x^2 + y^2) + 16 (x^2 + y^2) = 0
199: */
200: void f0_analytic_nonlinear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
201: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
202: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
203: PetscReal t, const PetscReal x[], PetscScalar f0[])
204: {
205: f0[0] = 16.0*(x[0]*x[0] + x[1]*x[1]);
206: }
208: /* gradU[comp*dim+d] = {u_x, u_y} or {u_x, u_y, u_z} */
209: void f1_analytic_nonlinear_u(PetscInt dim, PetscInt Nf, PetscInt NfAux,
210: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
211: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
212: PetscReal t, const PetscReal x[], PetscScalar f1[])
213: {
214: PetscScalar nu = 0.0;
215: PetscInt d;
216: for (d = 0; d < dim; ++d) nu += u_x[d]*u_x[d];
217: for (d = 0; d < dim; ++d) f1[d] = 0.5*nu*u_x[d];
218: }
220: /*
221: grad (u + eps w) - grad u = eps grad w
223: 1/2 |grad (u + eps w)|^2 grad (u + eps w) - 1/2 |grad u|^2 grad u
224: = 1/2 (|grad u|^2 + 2 eps <grad u,grad w>) (grad u + eps grad w) - 1/2 |grad u|^2 grad u
225: = 1/2 (eps |grad u|^2 grad w + 2 eps <grad u,grad w> grad u)
226: = eps (1/2 |grad u|^2 grad w + grad u <grad u,grad w>)
227: */
228: void g3_analytic_nonlinear_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux,
229: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
230: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
231: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscScalar g3[])
232: {
233: PetscScalar nu = 0.0;
234: PetscInt d, e;
235: for (d = 0; d < dim; ++d) nu += u_x[d]*u_x[d];
236: for (d = 0; d < dim; ++d) {
237: g3[d*dim+d] = 0.5*nu;
238: for (e = 0; e < dim; ++e) {
239: g3[d*dim+e] += u_x[d]*u_x[e];
240: }
241: }
242: }
244: /*
245: In 3D for Dirichlet conditions we use exact solution:
247: u = 2/3 (x^2 + y^2 + z^2)
248: f = 4
250: so that
252: -\Delta u + f = -2/3 * 6 + 4 = 0
254: For Neumann conditions, we have
256: -\nabla u \cdot -\hat z |_{z=0} = (2z)|_{z=0} = 0 (bottom)
257: -\nabla u \cdot \hat z |_{z=1} = -(2z)|_{z=1} = -2 (top)
258: -\nabla u \cdot -\hat y |_{y=0} = (2y)|_{y=0} = 0 (front)
259: -\nabla u \cdot \hat y |_{y=1} = -(2y)|_{y=1} = -2 (back)
260: -\nabla u \cdot -\hat x |_{x=0} = (2x)|_{x=0} = 0 (left)
261: -\nabla u \cdot \hat x |_{x=1} = -(2x)|_{x=1} = -2 (right)
263: Which we can express as
265: \nabla u \cdot \hat n|_\Gamma = {2 x, 2 y, 2z} \cdot \hat n = 2 (x + y + z)
266: */
267: static PetscErrorCode quadratic_u_3d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
268: {
269: *u = 2.0*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2])/3.0;
270: return 0;
271: }
275: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
276: {
277: const char *bcTypes[3] = {"neumann", "dirichlet", "none"};
278: const char *runTypes[4] = {"full", "exact", "test", "perf"};
279: const char *coeffTypes[4] = {"none", "analytic", "field", "nonlinear"};
280: PetscInt bc, run, coeff;
281: PetscBool flg;
285: options->debug = 0;
286: options->runType = RUN_FULL;
287: options->dim = 2;
288: options->filename[0] = '\0';
289: options->interpolate = PETSC_FALSE;
290: options->refinementLimit = 0.0;
291: options->bcType = DIRICHLET;
292: options->variableCoefficient = COEFF_NONE;
293: options->jacobianMF = PETSC_FALSE;
294: options->showInitial = PETSC_FALSE;
295: options->showSolution = PETSC_FALSE;
296: options->restart = PETSC_FALSE;
297: options->check = PETSC_FALSE;
298: options->viewHierarchy = PETSC_FALSE;
299: options->simplex = PETSC_TRUE;
301: PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
302: PetscOptionsInt("-debug", "The debugging level", "ex12.c", options->debug, &options->debug, NULL);
303: run = options->runType;
304: PetscOptionsEList("-run_type", "The run type", "ex12.c", runTypes, 3, runTypes[options->runType], &run, NULL);
306: options->runType = (RunType) run;
308: PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, NULL);
309: PetscOptionsString("-f", "Exodus.II filename to read", "ex12.c", options->filename, options->filename, sizeof(options->filename), &flg);
310: #if !defined(PETSC_HAVE_EXODUSII)
311: if (flg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "This option requires ExodusII support. Reconfigure using --download-exodusii");
312: #endif
313: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, NULL);
314: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, NULL);
315: bc = options->bcType;
316: PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,3,bcTypes[options->bcType],&bc,NULL);
317: options->bcType = (BCType) bc;
318: coeff = options->variableCoefficient;
319: PetscOptionsEList("-variable_coefficient","Type of variable coefficent","ex12.c",coeffTypes,4,coeffTypes[options->variableCoefficient],&coeff,NULL);
320: options->variableCoefficient = (CoeffType) coeff;
322: PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex12.c", options->jacobianMF, &options->jacobianMF, NULL);
323: PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex12.c", options->showInitial, &options->showInitial, NULL);
324: PetscOptionsBool("-show_solution", "Output the solution for verification", "ex12.c", options->showSolution, &options->showSolution, NULL);
325: PetscOptionsBool("-restart", "Read in the mesh and solution from a file", "ex12.c", options->restart, &options->restart, NULL);
326: PetscOptionsBool("-check", "Compare with default integration routines", "ex12.c", options->check, &options->check, NULL);
327: PetscOptionsBool("-dm_view_hierarchy", "View the coarsened hierarchy", "ex12.c", options->viewHierarchy, &options->viewHierarchy, NULL);
328: PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex12.c", options->simplex, &options->simplex, NULL);
329: PetscOptionsEnd();
330: PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
331: return(0);
332: }
336: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
337: {
338: DMLabel label;
342: DMCreateLabel(dm, name);
343: DMGetLabel(dm, name, &label);
344: DMPlexMarkBoundaryFaces(dm, label);
345: DMPlexLabelComplete(dm, label);
346: return(0);
347: }
351: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
352: {
353: PetscInt dim = user->dim;
354: const char *filename = user->filename;
355: PetscBool interpolate = user->interpolate;
356: PetscReal refinementLimit = user->refinementLimit;
357: size_t len;
361: PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
362: PetscStrlen(filename, &len);
363: if (!len) {
364: if (user->simplex) {
365: DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
366: }
367: else {
368: PetscInt cells[3] = {1, 1, 1}; /* coarse mesh is one cell; refine from there */
370: DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);
371: }
372: PetscObjectSetName((PetscObject) *dm, "Mesh");
373: } else {
374: DMPlexCreateFromFile(comm, filename, interpolate, dm);
375: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
376: }
377: {
378: DM refinedMesh = NULL;
379: DM distributedMesh = NULL;
381: /* Refine mesh using a volume constraint */
382: DMPlexSetRefinementLimit(*dm, refinementLimit);
383: DMRefine(*dm, comm, &refinedMesh);
384: if (refinedMesh) {
385: const char *name;
387: PetscObjectGetName((PetscObject) *dm, &name);
388: PetscObjectSetName((PetscObject) refinedMesh, name);
389: DMDestroy(dm);
390: *dm = refinedMesh;
391: }
392: /* Distribute mesh over processes */
393: DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
394: if (distributedMesh) {
395: DMDestroy(dm);
396: *dm = distributedMesh;
397: }
398: }
399: if (user->bcType == NEUMANN) {
400: DMLabel label;
402: DMCreateLabel(*dm, "boundary");
403: DMGetLabel(*dm, "boundary", &label);
404: DMPlexMarkBoundaryFaces(*dm, label);
405: } else if (user->bcType == DIRICHLET) {
406: PetscBool hasLabel;
408: DMHasLabel(*dm,"marker",&hasLabel);
409: if (!hasLabel) {CreateBCLabel(*dm, "marker");}
410: }
411: {
412: char convType[256];
413: PetscBool flg;
415: PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
416: PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex12",DMList,DMPLEX,convType,256,&flg);
417: PetscOptionsEnd();
418: if (flg) {
419: DM dmConv;
421: DMConvert(*dm,convType,&dmConv);
422: if (dmConv) {
423: DMDestroy(dm);
424: *dm = dmConv;
425: }
426: }
427: }
428: DMSetFromOptions(*dm);
429: DMViewFromOptions(*dm, NULL, "-dm_view");
430: if (user->viewHierarchy) {
431: DM cdm = *dm;
432: PetscInt i = 0;
433: char buf[256];
435: while (cdm) {
436: DMSetUp(cdm);
437: DMGetCoarseDM(cdm, &cdm);
438: ++i;
439: }
440: cdm = *dm;
441: while (cdm) {
442: PetscViewer viewer;
443: PetscBool isHDF5, isVTK;
445: --i;
446: PetscViewerCreate(comm,&viewer);
447: PetscViewerSetType(viewer,PETSCVIEWERHDF5);
448: PetscViewerSetOptionsPrefix(viewer,"hierarchy_");
449: PetscViewerSetFromOptions(viewer);
450: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);
451: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);
452: if (isHDF5) {
453: PetscSNPrintf(buf, 256, "ex12-%d.h5", i);
454: }
455: else if (isVTK) {
456: PetscSNPrintf(buf, 256, "ex12-%d.vtu", i);
457: PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);
458: }
459: else {
460: PetscSNPrintf(buf, 256, "ex12-%d", i);
461: }
462: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
463: PetscViewerFileSetName(viewer,buf);
464: DMView(cdm, viewer);
465: PetscViewerDestroy(&viewer);
466: DMGetCoarseDM(cdm, &cdm);
467: }
468: }
469: PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
470: return(0);
471: }
475: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
476: {
477: PetscDS prob;
481: DMGetDS(dm, &prob);
482: switch (user->variableCoefficient) {
483: case COEFF_NONE:
484: PetscDSSetResidual(prob, 0, f0_u, f1_u);
485: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
486: break;
487: case COEFF_ANALYTIC:
488: PetscDSSetResidual(prob, 0, f0_analytic_u, f1_analytic_u);
489: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_uu);
490: break;
491: case COEFF_FIELD:
492: PetscDSSetResidual(prob, 0, f0_analytic_u, f1_field_u);
493: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_field_uu);
494: break;
495: case COEFF_NONLINEAR:
496: PetscDSSetResidual(prob, 0, f0_analytic_nonlinear_u, f1_analytic_nonlinear_u);
497: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_analytic_nonlinear_uu);
498: break;
499: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid variable coefficient type %d", user->variableCoefficient);
500: }
501: switch (user->dim) {
502: case 2:
503: user->exactFuncs[0] = quadratic_u_2d;
504: if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
505: break;
506: case 3:
507: user->exactFuncs[0] = quadratic_u_3d;
508: if (user->bcType == NEUMANN) {PetscDSSetBdResidual(prob, 0, f0_bd_u, f1_bd_zero);}
509: break;
510: default:
511: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
512: }
513: return(0);
514: }
518: static PetscErrorCode SetupMaterial(DM dm, DM dmAux, AppCtx *user)
519: {
520: PetscErrorCode (*matFuncs[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {nu_2d};
521: Vec nu;
525: DMCreateLocalVector(dmAux, &nu);
526: DMProjectFunctionLocal(dmAux, 0.0, matFuncs, NULL, INSERT_ALL_VALUES, nu);
527: PetscObjectCompose((PetscObject) dm, "A", (PetscObject) nu);
528: VecDestroy(&nu);
529: return(0);
530: }
534: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
535: {
536: DM cdm = dm;
537: const PetscInt dim = user->dim;
538: const PetscInt id = 1;
539: PetscFE feAux = NULL;
540: PetscFE feBd = NULL;
541: PetscFE feCh = NULL;
542: PetscFE fe;
543: PetscDS prob;
544: PetscBool simplex = user->simplex;
548: /* Create finite element */
549: PetscFECreateDefault(dm, dim, 1, simplex, NULL, -1, &fe);
550: PetscObjectSetName((PetscObject) fe, "potential");
551: if (user->bcType == NEUMANN) {
552: PetscFECreateDefault(dm, dim-1, 1, simplex, "bd_", -1, &feBd);
553: PetscObjectSetName((PetscObject) feBd, "potential");
554: }
555: if (user->variableCoefficient == COEFF_FIELD) {
556: PetscQuadrature q;
558: PetscFECreateDefault(dm, dim, 1, simplex, "mat_", -1, &feAux);
559: PetscFEGetQuadrature(fe, &q);
560: PetscFESetQuadrature(feAux, q);
561: }
562: if (user->check) {PetscFECreateDefault(dm, dim, 1, simplex, "ch_", -1, &feCh);}
563: /* Set discretization and boundary conditions for each mesh */
564: while (cdm) {
565: DM coordDM;
567: DMGetDS(cdm, &prob);
568: PetscDSSetDiscretization(prob, 0, (PetscObject) fe);
569: PetscDSSetBdDiscretization(prob, 0, (PetscObject) feBd);
570: DMGetCoordinateDM(cdm,&coordDM);
571: if (feAux) {
572: DM dmAux;
573: PetscDS probAux;
575: DMClone(cdm, &dmAux);
576: DMSetCoordinateDM(dmAux, coordDM);
577: DMGetDS(dmAux, &probAux);
578: PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);
579: PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);
580: SetupMaterial(cdm, dmAux, user);
581: DMDestroy(&dmAux);
582: }
583: if (feCh) {
584: DM dmCh;
585: PetscDS probCh;
587: DMClone(cdm, &dmCh);
588: DMSetCoordinateDM(dmCh, coordDM);
589: DMGetDS(dmCh, &probCh);
590: PetscDSSetDiscretization(probCh, 0, (PetscObject) feCh);
591: PetscObjectCompose((PetscObject) dm, "dmCh", (PetscObject) dmCh);
592: DMDestroy(&dmCh);
593: }
594: if (user->bcType == DIRICHLET) {
595: PetscBool hasLabel;
597: DMHasLabel(cdm, "marker", &hasLabel);
598: if (!hasLabel) {CreateBCLabel(cdm, "marker");}
599: }
600: SetupProblem(cdm, user);
601: DMAddBoundary(cdm, user->bcType == DIRICHLET ? PETSC_TRUE : PETSC_FALSE, "wall", user->bcType == DIRICHLET ? "marker" : "boundary", 0, 0, NULL, (void (*)()) user->exactFuncs[0], 1, &id, user);
602: DMGetCoarseDM(cdm, &cdm);
603: }
604: PetscFEDestroy(&fe);
605: PetscFEDestroy(&feBd);
606: PetscFEDestroy(&feAux);
607: PetscFEDestroy(&feCh);
608: return(0);
609: }
611: #include petsc/private/petscimpl.h
615: /*@C
616: KSPMonitorError - Outputs the error at each iteration of an iterative solver.
618: Collective on KSP
620: Input Parameters:
621: + ksp - the KSP
622: . its - iteration number
623: . rnorm - 2-norm, preconditioned residual value (may be estimated).
624: - ctx - monitor context
626: Level: intermediate
628: .keywords: KSP, default, monitor, residual
629: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorDefault()
630: @*/
631: static PetscErrorCode KSPMonitorError(KSP ksp, PetscInt its, PetscReal rnorm, void *ctx)
632: {
633: AppCtx *user = (AppCtx *) ctx;
634: DM dm;
635: Vec du, r;
636: PetscInt level = 0;
637: PetscBool hasLevel;
638: PetscViewer viewer;
639: char buf[256];
643: KSPGetDM(ksp, &dm);
644: /* Calculate solution */
645: {
646: PC pc = user->pcmg; /* The MG PC */
647: DM fdm, cdm;
648: KSP fksp, cksp;
649: Vec fu, cu;
650: PetscInt levels, l;
652: KSPBuildSolution(ksp, NULL, &du);
653: PetscObjectComposedDataGetInt((PetscObject) ksp, PetscMGLevelId, level, hasLevel);
654: PCMGGetLevels(pc, &levels);
655: PCMGGetSmoother(pc, levels-1, &fksp);
656: KSPBuildSolution(fksp, NULL, &fu);
657: for (l = levels-1; l > level; --l) {
658: Mat R;
659: Vec s;
661: PCMGGetSmoother(pc, l-1, &cksp);
662: KSPGetDM(cksp, &cdm);
663: DMGetGlobalVector(cdm, &cu);
664: PCMGGetRestriction(pc, l, &R);
665: PCMGGetRScale(pc, l, &s);
666: MatRestrict(R, fu, cu);
667: VecPointwiseMult(cu, cu, s);
668: if (l < levels-1) {DMRestoreGlobalVector(fdm, &fu);}
669: fdm = cdm;
670: fu = cu;
671: }
672: if (levels-1 > level) {
673: VecAXPY(du, 1.0, cu);
674: DMRestoreGlobalVector(cdm, &cu);
675: }
676: }
677: /* Calculate error */
678: DMGetGlobalVector(dm, &r);
679: DMProjectFunction(dm, 0.0, user->exactFuncs, NULL, INSERT_ALL_VALUES, r);
680: VecAXPY(r,-1.0,du);
681: PetscObjectSetName((PetscObject) r, "solution error");
682: /* View error */
683: PetscSNPrintf(buf, 256, "ex12-%D.h5", level);
684: PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);
685: VecView(r, viewer);
686: /* Cleanup */
687: PetscViewerDestroy(&viewer);
688: DMRestoreGlobalVector(dm, &r);
689: return(0);
690: }
694: /*@C
695: SNESMonitorError - Outputs the error at each iteration of an iterative solver.
697: Collective on SNES
699: Input Parameters:
700: + snes - the SNES
701: . its - iteration number
702: . rnorm - 2-norm of residual
703: - ctx - user context
705: Level: intermediate
707: .keywords: SNES, nonlinear, default, monitor, norm
708: .seealso: SNESMonitorDefault(), SNESMonitorSet(), SNESMonitorSolution()
709: @*/
710: static PetscErrorCode SNESMonitorError(SNES snes, PetscInt its, PetscReal rnorm, void *ctx)
711: {
712: AppCtx *user = (AppCtx *) ctx;
713: DM dm;
714: Vec u, r;
715: PetscInt level;
716: PetscBool hasLevel;
717: PetscViewer viewer;
718: char buf[256];
722: SNESGetDM(snes, &dm);
723: /* Calculate error */
724: SNESGetSolution(snes, &u);
725: DMGetGlobalVector(dm, &r);
726: PetscObjectSetName((PetscObject) r, "solution error");
727: DMProjectFunction(dm, 0.0, user->exactFuncs, NULL, INSERT_ALL_VALUES, r);
728: VecAXPY(r, -1.0, u);
729: /* View error */
730: PetscObjectComposedDataGetInt((PetscObject) snes, PetscMGLevelId, level, hasLevel);
731: PetscSNPrintf(buf, 256, "ex12-%D.h5", level);
732: PetscViewerHDF5Open(PETSC_COMM_WORLD, buf, FILE_MODE_APPEND, &viewer);
733: VecView(r, viewer);
734: /* Cleanup */
735: PetscViewerDestroy(&viewer);
736: DMRestoreGlobalVector(dm, &r);
737: return(0);
738: }
742: int main(int argc, char **argv)
743: {
744: DM dm; /* Problem specification */
745: SNES snes; /* nonlinear solver */
746: Vec u,r; /* solution, residual vectors */
747: Mat A,J; /* Jacobian matrix */
748: MatNullSpace nullSpace; /* May be necessary for Neumann conditions */
749: AppCtx user; /* user-defined work context */
750: JacActionCtx userJ; /* context for Jacobian MF action */
751: PetscInt its; /* iterations for convergence */
752: PetscReal error = 0.0; /* L_2 error in the solution */
753: PetscBool isFAS;
756: PetscInitialize(&argc, &argv, NULL, help);
757: ProcessOptions(PETSC_COMM_WORLD, &user);
758: SNESCreate(PETSC_COMM_WORLD, &snes);
759: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
760: SNESSetDM(snes, dm);
761: DMSetApplicationContext(dm, &user);
763: PetscMalloc(1 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs);
764: SetupDiscretization(dm, &user);
766: DMCreateGlobalVector(dm, &u);
767: PetscObjectSetName((PetscObject) u, "potential");
768: VecDuplicate(u, &r);
770: DMSetMatType(dm,MATAIJ);
771: DMCreateMatrix(dm, &J);
772: if (user.jacobianMF) {
773: PetscInt M, m, N, n;
775: MatGetSize(J, &M, &N);
776: MatGetLocalSize(J, &m, &n);
777: MatCreate(PETSC_COMM_WORLD, &A);
778: MatSetSizes(A, m, n, M, N);
779: MatSetType(A, MATSHELL);
780: MatSetUp(A);
781: #if 0
782: MatShellSetOperation(A, MATOP_MULT, (void (*)(void))FormJacobianAction);
783: #endif
785: userJ.dm = dm;
786: userJ.J = J;
787: userJ.user = &user;
789: DMCreateLocalVector(dm, &userJ.u);
790: DMProjectFunctionLocal(dm, 0.0, user.exactFuncs, NULL, INSERT_BC_VALUES, userJ.u);
791: MatShellSetContext(A, &userJ);
792: } else {
793: A = J;
794: }
795: if (user.bcType == NEUMANN) {
796: MatNullSpaceCreate(PetscObjectComm((PetscObject) dm), PETSC_TRUE, 0, NULL, &nullSpace);
797: MatSetNullSpace(A, nullSpace);
798: }
800: DMPlexSetSNESLocalFEM(dm,&user,&user,&user);
801: SNESSetJacobian(snes, A, J, NULL, NULL);
803: SNESSetFromOptions(snes);
805: DMProjectFunction(dm, 0.0, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
806: if (user.restart) {
807: #if defined(PETSC_HAVE_HDF5)
808: PetscViewer viewer;
810: PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
811: PetscViewerSetType(viewer, PETSCVIEWERHDF5);
812: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
813: PetscViewerFileSetName(viewer, user.filename);
814: PetscViewerHDF5PushGroup(viewer, "/fields");
815: VecLoad(u, viewer);
816: PetscViewerHDF5PopGroup(viewer);
817: PetscViewerDestroy(&viewer);
818: #endif
819: }
820: if (user.showInitial) {
821: Vec lv;
822: DMGetLocalVector(dm, &lv);
823: DMGlobalToLocalBegin(dm, u, INSERT_VALUES, lv);
824: DMGlobalToLocalEnd(dm, u, INSERT_VALUES, lv);
825: DMPrintLocalVec(dm, "Local function", 1.0e-10, lv);
826: DMRestoreLocalVector(dm, &lv);
827: }
828: if (user.viewHierarchy) {
829: SNES lsnes;
830: KSP ksp;
831: PC pc;
832: PetscInt numLevels, l;
833: PetscBool isMG;
835: PetscObjectTypeCompare((PetscObject) snes, SNESFAS, &isFAS);
836: if (isFAS) {
837: SNESFASGetLevels(snes, &numLevels);
838: for (l = 0; l < numLevels; ++l) {
839: SNESFASGetCycleSNES(snes, l, &lsnes);
840: SNESMonitorSet(lsnes, SNESMonitorError, &user, NULL);
841: }
842: } else {
843: SNESGetKSP(snes, &ksp);
844: KSPGetPC(ksp, &pc);
845: PetscObjectTypeCompare((PetscObject) pc, PCMG, &isMG);
846: if (isMG) {
847: user.pcmg = pc;
848: PCMGGetLevels(pc, &numLevels);
849: for (l = 0; l < numLevels; ++l) {
850: PCMGGetSmootherDown(pc, l, &ksp);
851: KSPMonitorSet(ksp, KSPMonitorError, &user, NULL);
852: }
853: }
854: }
855: }
856: if (user.runType == RUN_FULL || user.runType == RUN_EXACT) {
857: PetscErrorCode (*initialGuess[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {zero};
859: if (user.runType == RUN_FULL) {
860: DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
861: }
862: if (user.debug) {
863: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
864: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
865: }
866: SNESSolve(snes, NULL, u);
867: SNESGetIterationNumber(snes, &its);
868: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);
869: DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);
870: if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
871: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
872: if (user.showSolution) {
873: PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
874: VecChop(u, 3.0e-9);
875: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
876: }
877: } else if (user.runType == RUN_PERF) {
878: PetscReal res = 0.0;
880: SNESComputeFunction(snes, u, r);
881: PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
882: VecChop(r, 1.0e-10);
883: VecNorm(r, NORM_2, &res);
884: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
885: } else {
886: PetscReal res = 0.0;
888: /* Check discretization error */
889: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
890: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
891: DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);
892: if (error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
893: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);}
894: /* Check residual */
895: SNESComputeFunction(snes, u, r);
896: PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
897: VecChop(r, 1.0e-10);
898: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
899: VecNorm(r, NORM_2, &res);
900: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
901: /* Check Jacobian */
902: {
903: Vec b;
905: SNESComputeJacobian(snes, u, A, A);
906: VecDuplicate(u, &b);
907: VecSet(r, 0.0);
908: SNESComputeFunction(snes, r, b);
909: MatMult(A, u, r);
910: VecAXPY(r, 1.0, b);
911: VecDestroy(&b);
912: PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
913: VecChop(r, 1.0e-10);
914: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
915: VecNorm(r, NORM_2, &res);
916: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
917: }
918: }
919: VecViewFromOptions(u, NULL, "-vec_view");
921: if (user.bcType == NEUMANN) {MatNullSpaceDestroy(&nullSpace);}
922: if (user.jacobianMF) {VecDestroy(&userJ.u);}
923: if (A != J) {MatDestroy(&A);}
924: MatDestroy(&J);
925: VecDestroy(&u);
926: VecDestroy(&r);
927: SNESDestroy(&snes);
928: DMDestroy(&dm);
929: PetscFree(user.exactFuncs);
930: PetscFinalize();
931: return 0;
932: }