Actual source code: ex31.c
petsc-3.4.5 2014-06-29
1: static char help[] = "Stokes Problem with Temperature 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: TODO for Mantle Convection:
7: - Variable viscosity
8: - Free-slip boundary condition on upper surface
9: - Parse Citcom input
11: The isoviscous Stokes problem, which we discretize using the finite
12: element method on an unstructured mesh. The weak form equations are
14: < \nabla v, \nabla u + {\nabla u}^T > - < \nabla\cdot v, p > + < v, f > = 0
15: < q, \nabla\cdot v > = 0
16: < \nabla t, \nabla T> = q_T
18: Boundary Conditions:
20: -bc_type dirichlet Dirichlet conditions on the entire boundary, coming from the exact solution functions
22: -dm_plex_separate_marker
23: -bc_type freeslip Dirichlet conditions on the normal velocity at each boundary
25: Discretization:
27: We use a Python script to generate a tabulation of the finite element basis
28: functions at quadrature points, which we put in a C header file. The generic
29: command would be:
31: bin/pythonscripts/PetscGenerateFEMQuadrature.py dim order dim 1 laplacian dim order 1 1 gradient src/snes/examples/tutorials/ex62.h
33: We can currently generate an arbitrary order Lagrange element. The underlying
34: FIAT code is capable of handling more exotic elements, but these have not been
35: tested with this code.
37: Field Data:
39: Sieve data is organized by point, and the closure operation just stacks up the
40: data from each sieve point in the closure. Thus, for a P_2-P_1 Stokes element, we
41: have
43: cl{e} = {f e_0 e_1 e_2 v_0 v_1 v_2}
44: 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}]
46: The problem here is that we would like to loop over each field separately for
47: integration. Therefore, the closure visitor in DMPlexVecGetClosure() reorders
48: the data so that each field is contiguous
50: 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}]
52: Likewise, DMPlexVecSetClosure() takes data partitioned by field, and correctly
53: puts it into the Sieve ordering.
54: */
56: #include <petscdmplex.h>
57: #include <petscsnes.h>
59: /*------------------------------------------------------------------------------
60: This code can be generated using 'bin/pythonscripts/PetscGenerateFEMQuadrature.py dim order dim 1 laplacian dim order 1 1 gradient dim order 1 1 identity src/snes/examples/tutorials/ex31.h'
61: -----------------------------------------------------------------------------*/
62: #include "ex31.h"
64: typedef enum {DIRICHLET, FREE_SLIP} BCType;
65: typedef enum {RUN_FULL, RUN_TEST} RunType;
66: typedef enum {FORCING_CONSTANT, FORCING_LINEAR, FORCING_CUBIC} ForcingType;
68: typedef struct {
69: DM dm; /* REQUIRED in order to use SNES evaluation functions */
70: PetscFEM fem; /* REQUIRED to use DMPlexComputeResidualFEM() */
71: PetscInt debug; /* The debugging level */
72: PetscMPIInt rank; /* The process rank */
73: PetscMPIInt numProcs; /* The number of processes */
74: RunType runType; /* Whether to run tests, or solve the full problem */
75: PetscBool jacobianMF; /* Whether to calculate the Jacobian action on the fly */
76: PetscLogEvent createMeshEvent;
77: PetscBool showInitial, showSolution;
78: /* Domain and mesh definition */
79: PetscInt dim; /* The topological mesh dimension */
80: PetscBool interpolate; /* Generate intermediate mesh elements */
81: PetscReal refinementLimit; /* The largest allowable cell volume */
82: char partitioner[2048]; /* The graph partitioner */
83: /* GPU partitioning */
84: PetscInt numBatches; /* The number of cell batches per kernel */
85: PetscInt numBlocks; /* The number of concurrent blocks per kernel */
86: /* Element quadrature */
87: PetscQuadrature q[NUM_FIELDS];
88: /* Problem definition */
89: void (*f0Funcs[NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]); /* f0_u(x,y,z), and f0_p(x,y,z) */
90: void (*f1Funcs[NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]); /* f1_u(x,y,z), and f1_p(x,y,z) */
91: void (*g0Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g0[]); /* g0_uu(x,y,z), g0_up(x,y,z), g0_pu(x,y,z), and g0_pp(x,y,z) */
92: void (*g1Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g1[]); /* g1_uu(x,y,z), g1_up(x,y,z), g1_pu(x,y,z), and g1_pp(x,y,z) */
93: void (*g2Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g2[]); /* g2_uu(x,y,z), g2_up(x,y,z), g2_pu(x,y,z), and g2_pp(x,y,z) */
94: void (*g3Funcs[NUM_FIELDS*NUM_FIELDS])(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[]); /* g3_uu(x,y,z), g3_up(x,y,z), g3_pu(x,y,z), and g3_pp(x,y,z) */
95: void (*exactFuncs[NUM_BASIS_COMPONENTS_TOTAL])(const PetscReal x[], PetscScalar *u); /* The exact solution function u(x,y,z), v(x,y,z), p(x,y,z), and T(x,y,z) */
96: BCType bcType; /* The type of boundary conditions */
97: ForcingType forcingType; /* The type of rhs */
98: } AppCtx;
100: void zero(const PetscReal coords[], PetscScalar *u)
101: {
102: *u = 0.0;
103: }
105: /*
106: In 2D, for constant forcing,
108: f_x = f_y = 3
110: we use the exact solution,
112: u = x^2 + y^2
113: v = 2 x^2 - 2xy
114: p = x + y - 1
115: T = x + y
117: so that
119: -\Delta u + \nabla p + f = <-4, -4> + <1, 1> + <3, 3> = 0
120: \nabla \cdot u = 2x - 2x = 0
121: -\Delta T + q_T = 0
122: */
123: void quadratic_u_2d(const PetscReal x[], PetscScalar *u)
124: {
125: *u = x[0]*x[0] + x[1]*x[1];
126: };
128: void quadratic_v_2d(const PetscReal x[], PetscScalar *v)
129: {
130: *v = 2.0*x[0]*x[0] - 2.0*x[0]*x[1];
131: };
133: void linear_p_2d(const PetscReal x[], PetscScalar *p)
134: {
135: *p = x[0] + x[1] - 1.0;
136: };
138: void linear_T_2d(const PetscReal x[], PetscScalar *T)
139: {
140: *T = x[0] + x[1];
141: };
143: /*
144: In 2D, for linear forcing,
146: f_x = 3 - 8y
147: f_y = -5 + 8x
149: we use the exact solution,
151: u = 2 x (x-1) (1 - 2 y)
152: v = -2 y (y-1) (1 - 2 x)
153: p = x + y - 1
154: T = x + y
156: so that
158: -\Delta u + \nabla p + f = <-4+8y, 4-8x> + <1, 1> + <3-8y, 8x-5> = 0
159: \nabla \cdot u = (4x-2) (1-2y) - (4y-2) (1-2x) = 0
160: -\Delta T + q_T = 0
161: */
162: void cubic_u_2d(const PetscReal x[], PetscScalar *u)
163: {
164: *u = 2.0*x[0]*(x[0]-1.0)*(1.0 - 2.0*x[1]);
165: };
167: void cubic_v_2d(const PetscReal x[], PetscScalar *v)
168: {
169: *v = -2.0*x[1]*(x[1]-1.0)*(1.0 - 2.0*x[0]);
170: };
172: /*
173: Let \sigma = (\nabla u + \nabla u^T) = < \sigma_{ij} >, where \sigma_{ij} = \sigma_{ji}
174: Then at the top and bottom (t = <1,0>),
175: <\sigma_{00}, \sigma_{01}> = 0 so \sigma_{00} = A(x,y) y(1-y) \sigma_{01} = B(x,y) y(1-y)
176: Using the left and right (t = <0,1>),
177: <\sigma_{10}, \sigma_{11}> = 0 so \sigma_{10} = C(x,y) x(1-x) \sigma_{11} = D(x,y) x(1-x)
178: Which means
179: \sigma_{00} = A(x,y) y(1-y) = 2 u_x
180: \sigma_{01} = E(x,y) x(1-x) y(1-y) = u_y + v_x
181: \sigma_{11} = D(x,y) x(1-x) = 2 v_y
182: Also we have
183: u(x=0,1) = 0 ==> u = A'(x,y) x(1-x)
184: v(y=0,1) = 0 ==> v = D'(x,y) y(1-y)
185: Thus we need
186: \int x - x^2 = x^2/2 - x^3/3 + C ==> 3 x^2 - 2 x^3 + 1 = 0 at x=0,1
187: so that
188: u = (3 x^2 - 2 x^3 + 1) y(1-y)
189: v = -(3 y^2 - 2 y^3 + 1) x(1-x)
190: u_x = 6 x(1-x) y(1-y)
191: v_y = -6 x(1-x) y(1-y)
192: u_xx = 6 (1-2x) y(1-y)
193: v_yy = -6 (1-2y) x(1-x)
195: In 2D, for cubic forcing,
197: f_x = -1 + 6 (1-2x) y(1-y)
198: f_y = -1 - 6 (1-2y) x(1-x)
200: we use the exact solution,
202: u = (3 x^2 - 2 x^3 + 1) y(1-y)
203: v = -(3 y^2 - 2 y^3 + 1) x(1-x)
204: p = x + y - 1
205: T = x + y
207: so that
209: -\Delta u + \nabla p + f = <-6 (1-2x) y(1-y), 6 (1-2y) x(1-x)> + <1, 1> + <-1 + 6 (1-2x) y(1-y), -1 - 6 (1-2y) x(1-x)> = 0
210: \nabla \cdot u = 6 x(1-x) y(1-y) -6 (1-2y) x(1-x) = 0
211: -\Delta T + q_T = 0
212: */
213: void quintic_u_2d(const PetscReal x[], PetscScalar *u)
214: {
215: *u = (3.0*x[0]*x[0] - 2.0*x[0]*x[0]*x[0] + 1.0)*x[1]*(1.0-x[1]);
216: };
218: void quintic_v_2d(const PetscReal x[], PetscScalar *v)
219: {
220: *v = -(3.0*x[1]*x[1] - 2.0*x[1]*x[1]*x[1] + 1.0)*x[0]*(1.0-x[0]);
221: };
223: void f0_u_constant(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
224: {
225: const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
226: PetscInt comp;
228: for (comp = 0; comp < Ncomp; ++comp) f0[comp] = 3.0;
229: }
231: void f0_u_linear_2d(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
232: {
233: f0[0] = 3.0 - 8.0*x[1];
234: f0[1] = -5.0 + 8.0*x[0];
235: }
237: void f0_u_cubic_2d(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
238: {
239: f0[0] = -1.0 + 6.0*(1.0 - 2.0*x[0])*x[1]*(1.0 - x[1]);
240: f0[1] = -1.0 - 6.0*(1.0 - 2.0*x[1])*x[0]*(1.0 - x[0]);
241: }
243: /* 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}
244: u[Ncomp] = {p} */
245: void f1_u(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[])
246: {
247: const PetscInt dim = SPATIAL_DIM_0;
248: const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
249: PetscInt comp, d;
251: for (comp = 0; comp < Ncomp; ++comp) {
252: for (d = 0; d < dim; ++d) {
253: /* f1[comp*dim+d] = 0.5*(gradU[comp*dim+d] + gradU[d*dim+comp]); */
254: f1[comp*dim+d] = gradU[comp*dim+d];
255: }
256: f1[comp*dim+comp] -= u[Ncomp];
257: }
258: }
260: /* 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} */
261: void f0_p(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
262: {
263: const PetscInt dim = SPATIAL_DIM_0;
264: PetscInt d;
266: f0[0] = 0.0;
267: for (d = 0; d < dim; ++d) f0[0] += gradU[d*dim+d];
268: }
270: void f1_p(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[])
271: {
272: const PetscInt dim = SPATIAL_DIM_0;
273: PetscInt d;
275: for (d = 0; d < dim; ++d) f1[d] = 0.0;
276: }
278: void f0_T(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[])
279: {
280: f0[0] = 0.0;
281: }
283: void f1_T(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[])
284: {
285: const PetscInt dim = SPATIAL_DIM_2;
286: const PetscInt off = SPATIAL_DIM_0*NUM_BASIS_COMPONENTS_0+SPATIAL_DIM_1*NUM_BASIS_COMPONENTS_1;
287: PetscInt d;
289: for (d = 0; d < dim; ++d) f1[d] = gradU[off+d];
290: }
292: /* < v_t, I t > */
293: void g0_TT(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g0[])
294: {
295: g0[0] = 1.0;
296: }
298: /* < q, \nabla\cdot v >
299: NcompI = 1, NcompJ = dim */
300: void g1_pu(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g1[])
301: {
302: const PetscInt dim = SPATIAL_DIM_0;
303: PetscInt d;
305: for (d = 0; d < dim; ++d) g1[d*dim+d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
306: }
308: /* -< \nabla\cdot v, p >
309: NcompI = dim, NcompJ = 1 */
310: void g2_up(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g2[])
311: {
312: const PetscInt dim = SPATIAL_DIM_0;
313: PetscInt d;
315: for (d = 0; d < dim; ++d) g2[d*dim+d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */
316: }
318: /* < \nabla v, \nabla u + {\nabla u}^T >
319: This just gives \nabla u, give the perdiagonal for the transpose */
320: void g3_uu(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[])
321: {
322: const PetscInt dim = SPATIAL_DIM_0;
323: const PetscInt Ncomp = NUM_BASIS_COMPONENTS_0;
324: PetscInt compI, d;
326: for (compI = 0; compI < Ncomp; ++compI) {
327: for (d = 0; d < dim; ++d) {
328: g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
329: }
330: }
331: }
333: /* < \nabla t, \nabla T + {\nabla u}^T >
334: This just gives \nabla T, give the perdiagonal for the transpose */
335: void g3_TT(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar g3[])
336: {
337: const PetscInt dim = SPATIAL_DIM_2;
338: const PetscInt Ncomp = NUM_BASIS_COMPONENTS_2;
339: PetscInt compI, d;
341: for (compI = 0; compI < Ncomp; ++compI) {
342: for (d = 0; d < dim; ++d) {
343: g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0;
344: }
345: }
346: }
348: /*
349: In 3D we use exact solution:
351: u = x^2 + y^2
352: v = y^2 + z^2
353: w = x^2 + y^2 - 2(x+y)z
354: p = x + y + z - 3/2
355: f_x = f_y = f_z = 3
357: so that
359: -\Delta u + \nabla p + f = <-4, -4, -4> + <1, 1, 1> + <3, 3, 3> = 0
360: \nabla \cdot u = 2x + 2y - 2(x + y) = 0
361: */
362: void quadratic_u_3d(const PetscReal x[], PetscScalar *u)
363: {
364: *u = x[0]*x[0] + x[1]*x[1];
365: };
367: void quadratic_v_3d(const PetscReal x[], PetscScalar *v)
368: {
369: *v = x[1]*x[1] + x[2]*x[2];
370: };
372: void quadratic_w_3d(const PetscReal x[], PetscScalar *w)
373: {
374: *w = x[0]*x[0] + x[1]*x[1] - 2.0*(x[0] + x[1])*x[2];
375: };
377: void linear_p_3d(const PetscReal x[], PetscScalar *p)
378: {
379: *p = x[0] + x[1] + x[2] - 1.5;
380: };
384: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
385: {
386: const char *bcTypes[2] = {"dirichlet", "freeslip"};
387: const char *forcingTypes[3] = {"constant", "linear", "cubic"};
388: const char *runTypes[2] = {"full", "test"};
389: PetscInt bc, forcing, run;
393: options->debug = 0;
394: options->runType = RUN_FULL;
395: options->dim = 2;
396: options->interpolate = PETSC_FALSE;
397: options->refinementLimit = 0.0;
398: options->bcType = DIRICHLET;
399: options->forcingType = FORCING_CONSTANT;
400: options->numBatches = 1;
401: options->numBlocks = 1;
402: options->jacobianMF = PETSC_FALSE;
403: options->showInitial = PETSC_FALSE;
404: options->showSolution = PETSC_TRUE;
406: options->fem.quad = (PetscQuadrature*) &options->q;
407: options->fem.f0Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->f0Funcs;
408: options->fem.f1Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->f1Funcs;
409: options->fem.g0Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g0Funcs;
410: options->fem.g1Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g1Funcs;
411: options->fem.g2Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g2Funcs;
412: options->fem.g3Funcs = (void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[])) &options->g3Funcs;
414: MPI_Comm_size(comm, &options->numProcs);
415: MPI_Comm_rank(comm, &options->rank);
416: PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
417: PetscOptionsInt("-debug", "The debugging level", "ex31.c", options->debug, &options->debug, NULL);
418: run = options->runType;
419: PetscOptionsEList("-run_type", "The run type", "ex31.c", runTypes, 2, runTypes[options->runType], &run, NULL);
421: options->runType = (RunType) run;
423: PetscOptionsInt("-dim", "The topological mesh dimension", "ex31.c", options->dim, &options->dim, NULL);
424: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex31.c", options->interpolate, &options->interpolate, NULL);
425: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex31.c", options->refinementLimit, &options->refinementLimit, NULL);
426: PetscStrcpy(options->partitioner, "chaco");
427: PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, NULL);
428: bc = options->bcType;
429: PetscOptionsEList("-bc_type","Type of boundary condition","ex31.c",bcTypes,2,bcTypes[options->bcType],&bc,NULL);
431: options->bcType = (BCType) bc;
432: forcing = options->forcingType;
434: PetscOptionsEList("-forcing_type","Type of forcing function","ex31.c",forcingTypes,3,forcingTypes[options->forcingType],&forcing,NULL);
436: options->forcingType = (ForcingType) forcing;
438: PetscOptionsInt("-gpu_batches", "The number of cell batches per kernel", "ex31.c", options->numBatches, &options->numBatches, NULL);
439: PetscOptionsInt("-gpu_blocks", "The number of concurrent blocks per kernel", "ex31.c", options->numBlocks, &options->numBlocks, NULL);
440: PetscOptionsBool("-jacobian_mf", "Calculate the action of the Jacobian on the fly", "ex31.c", options->jacobianMF, &options->jacobianMF, NULL);
441: PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex31.c", options->showInitial, &options->showInitial, NULL);
442: PetscOptionsBool("-show_solution", "Output the solution for verification", "ex31.c", options->showSolution, &options->showSolution, NULL);
443: PetscOptionsEnd();
445: PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
446: return(0);
447: };
451: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
452: {
453: Vec lv;
454: PetscInt p;
455: PetscMPIInt rank, numProcs;
459: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
460: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
461: DMGetLocalVector(dm, &lv);
462: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
463: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
464: PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
465: for (p = 0; p < numProcs; ++p) {
466: if (p == rank) {VecView(lv, PETSC_VIEWER_STDOUT_SELF);}
467: PetscBarrier((PetscObject) dm);
468: }
469: DMRestoreLocalVector(dm, &lv);
470: return(0);
471: }
475: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
476: {
477: PetscInt dim = user->dim;
478: PetscBool interpolate = user->interpolate;
479: PetscReal refinementLimit = user->refinementLimit;
480: const char *partitioner = user->partitioner;
484: PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
485: DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
486: {
487: DM refinedMesh = NULL;
488: DM distributedMesh = NULL;
490: /* Refine mesh using a volume constraint */
491: DMPlexSetRefinementLimit(*dm, refinementLimit);
492: DMRefine(*dm, comm, &refinedMesh);
493: if (refinedMesh) {
494: DMDestroy(dm);
495: *dm = refinedMesh;
496: }
497: /* Distribute mesh over processes */
498: DMPlexDistribute(*dm, partitioner, 0, &distributedMesh);
499: if (distributedMesh) {
500: DMDestroy(dm);
501: *dm = distributedMesh;
502: }
503: }
504: DMSetFromOptions(*dm);
505: PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
506: user->dm = *dm;
507: return(0);
508: }
512: PetscErrorCode PointOnBoundary_2D(const PetscScalar coords[], PetscBool onBd[])
513: {
514: const PetscInt corner = 0, bottom = 1, right = 2, top = 3, left = 4;
515: const PetscReal eps = 1.0e-10;
518: onBd[bottom] = PetscAbsScalar(coords[1]) < eps ? PETSC_TRUE : PETSC_FALSE;
519: onBd[right] = PetscAbsScalar(coords[0] - 1.0) < eps ? PETSC_TRUE : PETSC_FALSE;
520: onBd[top] = PetscAbsScalar(coords[1] - 1.0) < eps ? PETSC_TRUE : PETSC_FALSE;
521: onBd[left] = PetscAbsScalar(coords[0]) < eps ? PETSC_TRUE : PETSC_FALSE;
522: onBd[corner] = onBd[bottom] + onBd[right] + onBd[top] + onBd[left] > 1 ? PETSC_TRUE : PETSC_FALSE;
523: return(0);
524: }
528: PetscErrorCode CreateBoundaryPointIS_Square(DM dm, PetscInt *numBoundaries, PetscInt **numBoundaryConstraints, IS **boundaryPoints, IS **constraintIndices)
529: {
530: MPI_Comm comm;
531: PetscSection coordSection;
532: Vec coordinates;
533: PetscScalar *coords;
534: PetscInt vStart, vEnd;
535: IS bcPoints;
536: const PetscInt *points;
537: const PetscInt corner = 0, bottom = 1, right = 2, top = 3, left = 4;
538: PetscInt numBoundaryPoints[5] = {0, 0, 0, 0, 0}, bd, numPoints, p;
539: PetscInt *bdPoints[5], *idx;
543: PetscObjectGetComm((PetscObject)dm,&comm);
544: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
545: /* boundary 0: corners
546: boundary 1: bottom
547: boundary 2: right
548: boundary 3: top
549: boundary 4: left
550: */
551: *numBoundaries = 5;
553: PetscMalloc(*numBoundaries * sizeof(PetscInt), numBoundaryConstraints);
554: PetscMalloc(*numBoundaries * sizeof(IS), boundaryPoints);
555: PetscMalloc(*numBoundaries * sizeof(IS), constraintIndices);
557: /* Set number of constraints for each boundary */
558: (*numBoundaryConstraints)[corner] = 2;
559: (*numBoundaryConstraints)[bottom] = 1;
560: (*numBoundaryConstraints)[right] = 1;
561: (*numBoundaryConstraints)[top] = 1;
562: (*numBoundaryConstraints)[left] = 1;
563: /* Set local constraint indices for each boundary */
564: PetscMalloc((*numBoundaryConstraints)[corner] * sizeof(PetscInt), &idx);
565: idx[0] = 0; idx[1] = 1;
566: ISCreateGeneral(comm, (*numBoundaryConstraints)[corner], idx, PETSC_OWN_POINTER, &(*constraintIndices)[corner]);
567: PetscMalloc((*numBoundaryConstraints)[bottom] * sizeof(PetscInt), &idx);
568: idx[0] = 1;
569: ISCreateGeneral(comm, (*numBoundaryConstraints)[bottom], idx, PETSC_OWN_POINTER, &(*constraintIndices)[bottom]);
570: PetscMalloc((*numBoundaryConstraints)[right] * sizeof(PetscInt), &idx);
571: idx[0] = 0;
572: ISCreateGeneral(comm, (*numBoundaryConstraints)[right], idx, PETSC_OWN_POINTER, &(*constraintIndices)[right]);
573: PetscMalloc((*numBoundaryConstraints)[top] * sizeof(PetscInt), &idx);
574: idx[0] = 1;
575: ISCreateGeneral(comm, (*numBoundaryConstraints)[top], idx, PETSC_OWN_POINTER, &(*constraintIndices)[top]);
576: PetscMalloc((*numBoundaryConstraints)[left] * sizeof(PetscInt), &idx);
577: idx[0] = 0;
578: ISCreateGeneral(comm, (*numBoundaryConstraints)[left], idx, PETSC_OWN_POINTER, &(*constraintIndices)[left]);
580: /* Count points on each boundary */
581: DMPlexGetCoordinateSection(dm, &coordSection);
582: DMGetCoordinatesLocal(dm, &coordinates);
583: VecGetArray(coordinates, &coords);
584: DMPlexGetStratumIS(dm, "marker", 1, &bcPoints);
585: ISGetLocalSize(bcPoints, &numPoints);
586: ISGetIndices(bcPoints, &points);
587: for (p = 0; p < numPoints; ++p) {
588: PetscBool onBd[5];
589: PetscInt off, bd;
591: if ((points[p] >= vStart) && (points[p] < vEnd)) {
592: PetscSectionGetOffset(coordSection, points[p], &off);
593: PointOnBoundary_2D(&coords[off], onBd);
594: } else {
595: PetscInt *closure = NULL;
596: PetscInt closureSize, q, r;
598: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
599: /* Compress out non-vertices */
600: for (q = 0, r = 0; q < closureSize*2; q += 2) {
601: if ((closure[q] >= vStart) && (closure[q] < vEnd)) {
602: closure[r] = closure[q];
603: ++r;
604: }
605: }
606: closureSize = r;
607: for (q = 0; q < closureSize; ++q) {
608: PetscSectionGetOffset(coordSection, closure[q], &off);
609: PointOnBoundary_2D(&coords[off], onBd);
610: if (!onBd[corner]) break;
611: }
612: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
613: if (q == closureSize) SETERRQ1(comm, PETSC_ERR_PLIB, "Cannot handle face %d which has every vertex on a corner", points[p]);
614: }
616: for (bd = 0; bd < 5; ++bd) {
617: if (onBd[bd]) {
618: ++numBoundaryPoints[bd];
619: break;
620: }
621: }
622: }
623: /* Set points on each boundary */
624: for (bd = 0; bd < 5; ++bd) {
625: PetscMalloc(numBoundaryPoints[bd] * sizeof(PetscInt), &bdPoints[bd]);
626: numBoundaryPoints[bd] = 0;
627: }
628: for (p = 0; p < numPoints; ++p) {
629: PetscBool onBd[5];
630: PetscInt off, bd;
632: if ((points[p] >= vStart) && (points[p] < vEnd)) {
633: PetscSectionGetOffset(coordSection, points[p], &off);
634: PointOnBoundary_2D(&coords[off], onBd);
635: } else {
636: PetscInt *closure = NULL;
637: PetscInt closureSize, q, r;
639: DMPlexGetTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
640: /* Compress out non-vertices */
641: for (q = 0, r = 0; q < closureSize*2; q += 2) {
642: if ((closure[q] >= vStart) && (closure[q] < vEnd)) {
643: closure[r] = closure[q];
644: ++r;
645: }
646: }
647: closureSize = r;
648: for (q = 0; q < closureSize; ++q) {
649: PetscSectionGetOffset(coordSection, closure[q], &off);
650: PointOnBoundary_2D(&coords[off], onBd);
651: if (!onBd[corner]) break;
652: }
653: DMPlexRestoreTransitiveClosure(dm, points[p], PETSC_TRUE, &closureSize, &closure);
654: if (q == closureSize) SETERRQ1(comm, PETSC_ERR_PLIB, "Cannot handle face %d which has every vertex on a corner", points[p]);
655: }
657: for (bd = 0; bd < 5; ++bd) {
658: if (onBd[bd]) {
659: bdPoints[bd][numBoundaryPoints[bd]++] = points[p];
660: break;
661: }
662: }
663: }
664: VecRestoreArray(coordinates, &coords);
665: ISRestoreIndices(bcPoints, &points);
666: ISDestroy(&bcPoints);
667: for (bd = 0; bd < 5; ++bd) {
668: ISCreateGeneral(comm, numBoundaryPoints[bd], bdPoints[bd], PETSC_OWN_POINTER, &(*boundaryPoints)[bd]);
669: }
670: return(0);
671: }
675: PetscErrorCode CreateBoundaryPointIS_Cube(DM dm, PetscInt *numBoundaries, PetscInt **numBoundaryConstraints, IS **boundaryPoints, IS **constraintIndices)
676: {
678: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Just lazy");
679: return(0);
680: }
684: /* This will only work for the square/cube, but I think the interface is robust */
685: PetscErrorCode CreateBoundaryPointIS(DM dm, PetscInt *numBoundaries, PetscInt **numBoundaryConstraints, IS **boundaryPoints, IS **constraintIndices)
686: {
687: PetscInt dim;
691: DMPlexGetDimension(dm, &dim);
692: switch (dim) {
693: case 2:
694: CreateBoundaryPointIS_Square(dm, numBoundaries, numBoundaryConstraints, boundaryPoints, constraintIndices);
695: break;
696: case 3:
697: CreateBoundaryPointIS_Cube(dm, numBoundaries, numBoundaryConstraints, boundaryPoints, constraintIndices);
698: break;
699: default:
700: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "No boundary creatin routine for dimension %d", dim);
701: }
702: return(0);
703: }
707: PetscErrorCode SetupQuadrature(AppCtx *user)
708: {
710: user->fem.quad[0].numQuadPoints = NUM_QUADRATURE_POINTS_0;
711: user->fem.quad[0].quadPoints = points_0;
712: user->fem.quad[0].quadWeights = weights_0;
713: user->fem.quad[0].numBasisFuncs = NUM_BASIS_FUNCTIONS_0;
714: user->fem.quad[0].numComponents = NUM_BASIS_COMPONENTS_0;
715: user->fem.quad[0].basis = Basis_0;
716: user->fem.quad[0].basisDer = BasisDerivatives_0;
717: user->fem.quad[1].numQuadPoints = NUM_QUADRATURE_POINTS_1;
718: user->fem.quad[1].quadPoints = points_1;
719: user->fem.quad[1].quadWeights = weights_1;
720: user->fem.quad[1].numBasisFuncs = NUM_BASIS_FUNCTIONS_1;
721: user->fem.quad[1].numComponents = NUM_BASIS_COMPONENTS_1;
722: user->fem.quad[1].basis = Basis_1;
723: user->fem.quad[1].basisDer = BasisDerivatives_1;
724: user->fem.quad[2].numQuadPoints = NUM_QUADRATURE_POINTS_2;
725: user->fem.quad[2].quadPoints = points_2;
726: user->fem.quad[2].quadWeights = weights_2;
727: user->fem.quad[2].numBasisFuncs = NUM_BASIS_FUNCTIONS_2;
728: user->fem.quad[2].numComponents = NUM_BASIS_COMPONENTS_2;
729: user->fem.quad[2].basis = Basis_2;
730: user->fem.quad[2].basisDer = BasisDerivatives_2;
731: return(0);
732: }
736: /*
737: There is a problem here with uninterpolated meshes. The index in numDof[] is not dimension in this case,
738: but sieve depth.
739: */
740: PetscErrorCode SetupSection(DM dm, AppCtx *user)
741: {
742: PetscSection section;
743: const PetscInt numFields = NUM_FIELDS;
744: PetscInt dim = user->dim;
745: PetscInt numBC = 0;
746: PetscInt numComp[NUM_FIELDS] = {NUM_BASIS_COMPONENTS_0, NUM_BASIS_COMPONENTS_1, NUM_BASIS_COMPONENTS_2};
747: PetscInt bcFields[2] = {0, 2};
748: IS bcPoints[2] = {NULL, NULL};
749: PetscInt numDof[NUM_FIELDS*(SPATIAL_DIM_0+1)];
750: PetscInt f, d;
751: PetscBool view;
755: if (dim != SPATIAL_DIM_0) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Spatial dimension %d should be %d", dim, SPATIAL_DIM_0);
756: if (dim != SPATIAL_DIM_1) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Spatial dimension %d should be %d", dim, SPATIAL_DIM_1);
757: for (d = 0; d <= dim; ++d) {
758: numDof[0*(dim+1)+d] = numDof_0[d];
759: numDof[1*(dim+1)+d] = numDof_1[d];
760: numDof[2*(dim+1)+d] = numDof_2[d];
761: }
762: for (f = 0; f < numFields; ++f) {
763: for (d = 1; d < dim; ++d) {
764: if ((numDof[f*(dim+1)+d] > 0) && !user->interpolate) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces.");
765: }
766: }
767: if (user->bcType == FREE_SLIP) {
768: PetscInt numBoundaries, b;
769: PetscInt *numBoundaryConstraints;
770: IS *boundaryPoints, *constraintIndices;
772: DMPlexCreateSectionInitial(dm, dim, numFields, numComp, numDof, §ion);
773: /* Velocity conditions */
774: CreateBoundaryPointIS(dm, &numBoundaries, &numBoundaryConstraints, &boundaryPoints, &constraintIndices);
775: for (b = 0; b < numBoundaries; ++b) {
776: DMPlexCreateSectionBCDof(dm, 1, &bcFields[0], &boundaryPoints[b], numBoundaryConstraints[b], section);
777: }
778: /* Temperature conditions */
779: DMPlexGetStratumIS(dm, "marker", 1, &bcPoints[0]);
780: DMPlexCreateSectionBCDof(dm, 1, &bcFields[1], &bcPoints[0], PETSC_DETERMINE, section);
781: PetscSectionSetUp(section);
782: for (b = 0; b < numBoundaries; ++b) {
783: DMPlexCreateSectionBCIndicesField(dm, bcFields[0], boundaryPoints[b], constraintIndices[b], section);
784: }
785: DMPlexCreateSectionBCIndicesField(dm, bcFields[1], bcPoints[0], NULL, section);
786: DMPlexCreateSectionBCIndices(dm, section);
787: } else {
788: if (user->bcType == DIRICHLET) {
789: numBC = 2;
790: DMPlexGetStratumIS(dm, "marker", 1, &bcPoints[0]);
791: bcPoints[1] = bcPoints[0];
792: PetscObjectReference((PetscObject) bcPoints[1]);
793: }
794: DMPlexCreateSection(dm, dim, numFields, numComp, numDof, numBC, bcFields, bcPoints, §ion);
795: }
796: PetscSectionSetFieldName(section, 0, "velocity");
797: PetscSectionSetFieldName(section, 1, "pressure");
798: PetscSectionSetFieldName(section, 2, "temperature");
799: DMSetDefaultSection(dm, section);
800: ISDestroy(&bcPoints[0]);
801: ISDestroy(&bcPoints[1]);
802: PetscOptionsHasName(((PetscObject) dm)->prefix, "-section_view", &view);
803: if ((user->bcType == FREE_SLIP) && view) {
804: PetscSection s, gs;
806: DMGetDefaultSection(dm, &s);
807: PetscSectionView(s, PETSC_VIEWER_STDOUT_WORLD);
808: DMGetDefaultGlobalSection(dm, &gs);
809: PetscSectionView(gs, PETSC_VIEWER_STDOUT_WORLD);
810: }
811: return(0);
812: }
816: PetscErrorCode SetupExactSolution(DM dm, AppCtx *user)
817: {
818: PetscFEM *fem = &user->fem;
822: switch (user->forcingType) {
823: case FORCING_CONSTANT:
824: if (user->bcType == FREE_SLIP) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Constant forcing is incompatible with freeslip boundary conditions");
825: fem->f0Funcs[0] = f0_u_constant;
826: break;
827: case FORCING_LINEAR:
828: switch (user->bcType) {
829: case DIRICHLET:
830: case FREE_SLIP:
831: switch (user->dim) {
832: case 2:
833: fem->f0Funcs[0] = f0_u_linear_2d;
834: break;
835: default:
836: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
837: }
838: break;
839: default:
840: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary condition type %d", user->bcType);
841: }
842: break;
843: case FORCING_CUBIC:
844: switch (user->bcType) {
845: case DIRICHLET:
846: case FREE_SLIP:
847: switch (user->dim) {
848: case 2:
849: fem->f0Funcs[0] = f0_u_cubic_2d;
850: break;
851: default:
852: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
853: }
854: break;
855: default:
856: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary condition type %d", user->bcType);
857: }
858: break;
859: }
860: fem->f0Funcs[1] = f0_p;
861: fem->f0Funcs[2] = f0_T;
862: fem->f1Funcs[0] = f1_u;
863: fem->f1Funcs[1] = f1_p;
864: fem->f1Funcs[2] = f1_T;
865: fem->g0Funcs[0] = NULL;
866: fem->g0Funcs[1] = NULL;
867: fem->g0Funcs[2] = NULL;
868: fem->g0Funcs[3] = NULL;
869: fem->g0Funcs[4] = NULL;
870: fem->g0Funcs[5] = NULL;
871: fem->g0Funcs[6] = NULL;
872: fem->g0Funcs[7] = NULL;
873: fem->g0Funcs[8] = NULL;
874: fem->g1Funcs[0] = NULL;
875: fem->g1Funcs[1] = NULL;
876: fem->g1Funcs[2] = NULL;
877: fem->g1Funcs[3] = g1_pu; /* < q, \nabla\cdot v > */
878: fem->g1Funcs[4] = NULL;
879: fem->g1Funcs[5] = NULL;
880: fem->g1Funcs[6] = NULL;
881: fem->g1Funcs[7] = NULL;
882: fem->g1Funcs[8] = NULL;
883: fem->g2Funcs[0] = NULL;
884: fem->g2Funcs[1] = g2_up; /* < \nabla\cdot v, p > */
885: fem->g2Funcs[2] = NULL;
886: fem->g2Funcs[3] = NULL;
887: fem->g2Funcs[4] = NULL;
888: fem->g2Funcs[5] = NULL;
889: fem->g2Funcs[6] = NULL;
890: fem->g2Funcs[7] = NULL;
891: fem->g2Funcs[8] = NULL;
892: fem->g3Funcs[0] = g3_uu; /* < \nabla v, \nabla u + {\nabla u}^T > */
893: fem->g3Funcs[1] = NULL;
894: fem->g3Funcs[2] = NULL;
895: fem->g3Funcs[3] = NULL;
896: fem->g3Funcs[4] = NULL;
897: fem->g3Funcs[5] = NULL;
898: fem->g3Funcs[6] = NULL;
899: fem->g3Funcs[7] = NULL;
900: fem->g3Funcs[8] = g3_TT; /* < \nabla t, \nabla T + {\nabla T}^T > */
901: switch (user->forcingType) {
902: case FORCING_CONSTANT:
903: switch (user->bcType) {
904: case DIRICHLET:
905: switch (user->dim) {
906: case 2:
907: user->exactFuncs[0] = quadratic_u_2d;
908: user->exactFuncs[1] = quadratic_v_2d;
909: user->exactFuncs[2] = linear_p_2d;
910: user->exactFuncs[3] = linear_T_2d;
911: break;
912: case 3:
913: user->exactFuncs[0] = quadratic_u_3d;
914: user->exactFuncs[1] = quadratic_v_3d;
915: user->exactFuncs[2] = quadratic_w_3d;
916: user->exactFuncs[3] = linear_p_3d;
917: user->exactFuncs[4] = linear_T_2d;
918: break;
919: default:
920: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
921: }
922: break;
923: default:
924: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary condition type %d", user->bcType);
925: }
926: break;
927: case FORCING_LINEAR:
928: switch (user->bcType) {
929: case DIRICHLET:
930: switch (user->dim) {
931: case 2:
932: user->exactFuncs[0] = cubic_u_2d;
933: user->exactFuncs[1] = cubic_v_2d;
934: user->exactFuncs[2] = linear_p_2d;
935: user->exactFuncs[3] = linear_T_2d;
936: break;
937: default:
938: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
939: }
940: break;
941: default:
942: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary condition type %d", user->bcType);
943: }
944: break;
945: case FORCING_CUBIC:
946: switch (user->bcType) {
947: case DIRICHLET:
948: case FREE_SLIP:
949: switch (user->dim) {
950: case 2:
951: user->exactFuncs[0] = quintic_u_2d;
952: user->exactFuncs[1] = quintic_v_2d;
953: user->exactFuncs[2] = linear_p_2d;
954: user->exactFuncs[3] = linear_T_2d;
955: break;
956: default:
957: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
958: }
959: break;
960: default:
961: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid boundary condition type %d", user->bcType);
962: }
963: break;
964: default:
965: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid forcing type %d", user->forcingType);
966: }
967: DMPlexSetFEMIntegration(dm, FEMIntegrateResidualBatch, FEMIntegrateJacobianActionBatch, FEMIntegrateJacobianBatch);
968: return(0);
969: }
973: /*
974: . field - The field whose diagonal block (of the Jacobian) has this null space
975: */
976: PetscErrorCode CreateNullSpaces(DM dm, PetscInt field, MatNullSpace *nullSpace)
977: {
978: AppCtx *user;
979: Vec nullVec, localNullVec;
980: PetscSection section;
981: PetscScalar *a;
982: PetscInt pressure = field;
983: PetscInt pStart, pEnd, p;
987: DMGetApplicationContext(dm, (void**) &user);
988: DMGetGlobalVector(dm, &nullVec);
989: DMGetLocalVector(dm, &localNullVec);
990: VecSet(nullVec, 0.0);
991: /* Put a constant in for all pressures */
992: DMGetDefaultSection(dm, §ion);
993: PetscSectionGetChart(section, &pStart, &pEnd);
994: VecGetArray(localNullVec, &a);
995: for (p = pStart; p < pEnd; ++p) {
996: PetscInt fDim, off, d;
998: PetscSectionGetFieldDof(section, p, pressure, &fDim);
999: PetscSectionGetFieldOffset(section, p, pressure, &off);
1000: for (d = 0; d < fDim; ++d) a[off+d] = 1.0;
1001: }
1002: VecRestoreArray(localNullVec, &a);
1003: DMLocalToGlobalBegin(dm, localNullVec, INSERT_VALUES, nullVec);
1004: DMLocalToGlobalEnd(dm, localNullVec, INSERT_VALUES, nullVec);
1005: DMRestoreLocalVector(dm, &localNullVec);
1006: VecNormalize(nullVec, NULL);
1007: if (user->debug) {
1008: PetscPrintf(PetscObjectComm((PetscObject)dm), "Pressure Null Space\n");
1009: VecView(nullVec, PETSC_VIEWER_STDOUT_WORLD);
1010: }
1011: MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &nullVec, nullSpace);
1012: DMRestoreGlobalVector(dm, &nullVec);
1013: return(0);
1014: }
1018: /*
1019: FormJacobianAction - Form the global Jacobian action Y = JX from the global input X
1021: Input Parameters:
1022: + mat - The Jacobian shell matrix
1023: - X - Global input vector
1025: Output Parameter:
1026: . Y - Local output vector
1028: Note:
1029: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1030: like a GPU, or vectorize on a multicore machine.
1032: .seealso: FormJacobianActionLocal()
1033: */
1034: PetscErrorCode FormJacobianAction(Mat J, Vec X, Vec Y)
1035: {
1036: JacActionCtx *ctx;
1037: DM dm;
1038: Vec dummy, localX, localY;
1039: PetscInt N, n;
1046: MatShellGetContext(J, &ctx);
1047: dm = ctx->dm;
1049: /* determine whether X = localX */
1050: DMGetLocalVector(dm, &dummy);
1051: DMGetLocalVector(dm, &localX);
1052: DMGetLocalVector(dm, &localY);
1053: /* TODO: THIS dummy restore is necessary here so that the first available local vector has boundary conditions in it
1054: I think the right thing to do is have the user put BC into a local vector and give it to us
1055: */
1056: DMRestoreLocalVector(dm, &dummy);
1057: VecGetSize(X, &N);
1058: VecGetSize(localX, &n);
1060: if (n != N) { /* X != localX */
1061: VecSet(localX, 0.0);
1062: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
1063: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
1064: } else {
1065: DMRestoreLocalVector(dm, &localX);
1066: localX = X;
1067: }
1068: DMPlexComputeJacobianActionFEM(dm, J, localX, localY, ctx->user);
1069: if (n != N) {
1070: DMRestoreLocalVector(dm, &localX);
1071: }
1072: VecSet(Y, 0.0);
1073: DMLocalToGlobalBegin(dm, localY, ADD_VALUES, Y);
1074: DMLocalToGlobalEnd(dm, localY, ADD_VALUES, Y);
1075: DMRestoreLocalVector(dm, &localY);
1076: if (0) {
1077: Vec r;
1078: PetscReal norm;
1080: VecDuplicate(X, &r);
1081: MatMult(ctx->J, X, r);
1082: VecAXPY(r, -1.0, Y);
1083: VecNorm(r, NORM_2, &norm);
1084: if (norm > 1.0e-8) {
1085: PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Input:\n");
1086: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
1087: PetscPrintf(PETSC_COMM_WORLD, "Jacobian Action Result:\n");
1088: VecView(Y, PETSC_VIEWER_STDOUT_WORLD);
1089: PetscPrintf(PETSC_COMM_WORLD, "Difference:\n");
1090: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
1091: SETERRQ1(PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONG, "The difference with assembled multiply is too large %g", norm);
1092: }
1093: VecDestroy(&r);
1094: }
1095: return(0);
1096: }
1100: int main(int argc, char **argv)
1101: {
1102: MPI_Comm comm;
1103: SNES snes; /* nonlinear solver */
1104: Vec u,r; /* solution, residual vectors */
1105: Mat A,J; /* Jacobian matrix */
1106: MatNullSpace nullSpace = 0; /* May be necessary for pressure */
1107: AppCtx user; /* user-defined work context */
1108: JacActionCtx userJ; /* context for Jacobian MF action */
1109: PetscInt its; /* iterations for convergence */
1110: PetscReal error = 0.0; /* L_2 error in the solution */
1111: const PetscInt numComponents = NUM_BASIS_COMPONENTS_TOTAL;
1114: PetscInitialize(&argc, &argv, NULL, help);
1115: comm = PETSC_COMM_WORLD;
1116: ProcessOptions(comm, &user);
1117: SNESCreate(comm, &snes);
1118: CreateMesh(comm, &user, &user.dm);
1119: SNESSetDM(snes, user.dm);
1120: DMSetApplicationContext(user.dm, &user);
1122: SetupExactSolution(user.dm, &user);
1123: SetupQuadrature(&user);
1124: SetupSection(user.dm, &user);
1126: DMCreateGlobalVector(user.dm, &u);
1127: PetscObjectSetName((PetscObject) u, "solution");
1128: VecDuplicate(u, &r);
1130: DMCreateMatrix(user.dm, MATAIJ, &J);
1131: if (user.jacobianMF) {
1132: PetscInt M, m, N, n;
1134: MatGetSize(J, &M, &N);
1135: MatGetLocalSize(J, &m, &n);
1136: MatCreate(comm, &A);
1137: MatSetSizes(A, m, n, M, N);
1138: MatSetType(A, MATSHELL);
1139: MatSetUp(A);
1140: MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) FormJacobianAction);
1142: userJ.dm = user.dm;
1143: userJ.J = J;
1144: userJ.user = &user;
1146: DMCreateLocalVector(user.dm, &userJ.u);
1147: MatShellSetContext(A, &userJ);
1148: } else {
1149: A = J;
1150: }
1151: DMSetNullSpaceConstructor(user.dm, 1, CreateNullSpaces);
1153: DMSNESSetFunctionLocal(user.dm, (PetscErrorCode (*)(DM,Vec,Vec,void*))DMPlexComputeResidualFEM,&user);
1154: DMSNESSetJacobianLocal(user.dm, (PetscErrorCode (*)(DM,Vec,Mat,Mat,MatStructure*,void*))DMPlexComputeJacobianFEM,&user);
1156: SNESSetFromOptions(snes);
1158: {
1159: KSP ksp; PC pc; Vec crd_vec;
1160: const PetscScalar *v;
1161: PetscInt i,k,j,mlocal;
1162: PetscReal *coords;
1164: SNESGetKSP(snes, &ksp);
1165: KSPGetPC(ksp, &pc);
1166: DMGetCoordinatesLocal(user.dm, &crd_vec);
1167: VecGetLocalSize(crd_vec,&mlocal);
1168: PetscMalloc(SPATIAL_DIM_0*mlocal*sizeof(*coords),&coords);
1169: VecGetArrayRead(crd_vec,&v);
1170: for (k=j=0; j<mlocal; j++) {
1171: for (i=0; i<SPATIAL_DIM_0; i++,k++) {
1172: coords[k] = PetscRealPart(v[k]);
1173: }
1174: }
1175: VecRestoreArrayRead(crd_vec,&v);
1176: PCSetCoordinates(pc, SPATIAL_DIM_0, mlocal, coords);
1177: PetscFree(coords);
1178: }
1180: DMPlexProjectFunction(user.dm, numComponents, user.exactFuncs, INSERT_ALL_VALUES, u);
1181: if (user.showInitial) {DMVecViewLocal(user.dm, u, PETSC_VIEWER_STDOUT_SELF);}
1182: if (user.runType == RUN_FULL) {
1183: PetscScalar (*initialGuess[numComponents])(const PetscReal x[]);
1184: PetscInt c;
1186: for (c = 0; c < numComponents; ++c) initialGuess[c] = zero;
1187: DMPlexProjectFunction(user.dm, numComponents, initialGuess, INSERT_VALUES, u);
1188: if (user.showInitial) {DMVecViewLocal(user.dm, u, PETSC_VIEWER_STDOUT_SELF);}
1189: if (user.debug) {
1190: PetscPrintf(comm, "Initial guess\n");
1191: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
1192: }
1193: SNESSolve(snes, NULL, u);
1194: SNESGetIterationNumber(snes, &its);
1195: PetscPrintf(comm, "Number of SNES iterations = %D\n", its);
1196: DMPlexComputeL2Diff(user.dm, user.q, user.exactFuncs, u, &error);
1197: PetscPrintf(comm, "L_2 Error: %.3g\n", error);
1198: if (user.showSolution) {
1199: PetscPrintf(comm, "Solution\n");
1200: VecChop(u, 3.0e-9);
1201: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
1202: }
1203: } else {
1204: PetscReal res = 0.0;
1206: /* Check discretization error */
1207: PetscPrintf(comm, "Initial guess\n");
1208: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
1209: DMPlexComputeL2Diff(user.dm, user.q, user.exactFuncs, u, &error);
1210: PetscPrintf(comm, "L_2 Error: %g\n", error);
1211: /* Check residual */
1212: SNESComputeFunction(snes, u, r);
1213: PetscPrintf(comm, "Initial Residual\n");
1214: VecChop(r, 1.0e-10);
1215: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
1216: VecNorm(r, NORM_2, &res);
1217: PetscPrintf(comm, "L_2 Residual: %g\n", res);
1218: /* Check Jacobian */
1219: {
1220: Vec b;
1221: MatStructure flag;
1222: MatNullSpace nullSpace2;
1223: PetscBool isNull;
1225: CreateNullSpaces(user.dm, 1, &nullSpace2);
1226: MatNullSpaceTest(nullSpace2, J, &isNull);
1227: if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
1228: MatNullSpaceDestroy(&nullSpace2);
1230: SNESComputeJacobian(snes, u, &A, &A, &flag);
1231: VecDuplicate(u, &b);
1232: VecSet(r, 0.0);
1233: SNESComputeFunction(snes, r, b);
1234: MatMult(A, u, r);
1235: VecAXPY(r, 1.0, b);
1236: VecDestroy(&b);
1237: PetscPrintf(comm, "Au - b = Au + F(0)\n");
1238: VecChop(r, 1.0e-10);
1239: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
1240: VecNorm(r, NORM_2, &res);
1241: PetscPrintf(comm, "Linear L_2 Residual: %g\n", res);
1242: }
1243: }
1245: if (user.runType == RUN_FULL) {
1246: PetscContainer c;
1247: PetscSection section;
1248: Vec sol;
1249: PetscViewer viewer;
1250: const char *name;
1252: PetscViewerCreate(comm, &viewer);
1253: PetscViewerSetType(viewer, PETSCVIEWERVTK);
1254: PetscViewerFileSetName(viewer, "ex31_sol.vtk");
1255: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
1256: DMGetLocalVector(user.dm, &sol);
1257: PetscObjectGetName((PetscObject) u, &name);
1258: PetscObjectSetName((PetscObject) sol, name);
1259: DMGlobalToLocalBegin(user.dm, u, INSERT_VALUES, sol);
1260: DMGlobalToLocalEnd(user.dm, u, INSERT_VALUES, sol);
1261: DMGetDefaultSection(user.dm, §ion);
1262: PetscObjectReference((PetscObject) user.dm); /* Needed because viewer destroys the DM */
1263: PetscViewerVTKAddField(viewer, (PetscObject) user.dm, DMPlexVTKWriteAll, PETSC_VTK_POINT_FIELD, (PetscObject) sol);
1264: PetscObjectReference((PetscObject) sol); /* Needed because viewer destroys the Vec */
1265: PetscContainerCreate(comm, &c);
1266: PetscContainerSetPointer(c, section);
1267: PetscObjectCompose((PetscObject) sol, "section", (PetscObject) c);
1268: PetscContainerDestroy(&c);
1269: DMRestoreLocalVector(user.dm, &sol);
1270: PetscViewerDestroy(&viewer);
1271: }
1273: MatNullSpaceDestroy(&nullSpace);
1274: if (user.jacobianMF) {
1275: VecDestroy(&userJ.u);
1276: }
1277: if (A != J) {
1278: MatDestroy(&A);
1279: }
1280: MatDestroy(&J);
1281: VecDestroy(&u);
1282: VecDestroy(&r);
1283: SNESDestroy(&snes);
1284: DMDestroy(&user.dm);
1285: PetscFinalize();
1286: return 0;
1287: }