Actual source code: ex76.c
1: static char help[] = "Low Mach Flow in 2d and 3d channels with finite elements.\n\
2: We solve the Low Mach flow problem in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\n\n";
5: /*F
6: This Low Mach flow is a steady-state isoviscous Navier-Stokes flow. We discretize using the
7: finite element method on an unstructured mesh. The weak form equations are
9: \begin{align*}
10: < q, \nabla\cdot u > = 0
11: <v, u \cdot \nabla u> + < \nabla v, \nu (\nabla u + {\nabla u}^T) > - < \nabla\cdot v, p > - < v, f > = 0
12: < w, u \cdot \nabla T > - < \nabla w, \alpha \nabla T > - < w, Q > = 0
13: \end{align*}
15: where $\nu$ is the kinematic viscosity and $\alpha$ is thermal diffusivity.
17: For visualization, use
19: -dm_view hdf5:$PWD/sol.h5 -sol_vec_view hdf5:$PWD/sol.h5::append -exact_vec_view hdf5:$PWD/sol.h5::append
20: F*/
22: #include <petscdmplex.h>
23: #include <petscsnes.h>
24: #include <petscds.h>
25: #include <petscbag.h>
27: typedef enum {
28: SOL_QUADRATIC,
29: SOL_CUBIC,
30: NUM_SOL_TYPES
31: } SolType;
32: const char *solTypes[NUM_SOL_TYPES + 1] = {"quadratic", "cubic", "unknown"};
34: typedef struct {
35: PetscReal nu; /* Kinematic viscosity */
36: PetscReal theta; /* Angle of pipe wall to x-axis */
37: PetscReal alpha; /* Thermal diffusivity */
38: PetscReal T_in; /* Inlet temperature*/
39: } Parameter;
41: typedef struct {
42: PetscBool showError;
43: PetscBag bag;
44: SolType solType;
45: } AppCtx;
47: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
48: {
49: PetscInt d;
50: for (d = 0; d < Nc; ++d) u[d] = 0.0;
51: return 0;
52: }
54: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
55: {
56: PetscInt d;
57: for (d = 0; d < Nc; ++d) u[d] = 1.0;
58: return 0;
59: }
61: /*
62: CASE: quadratic
63: In 2D we use exact solution:
65: u = x^2 + y^2
66: v = 2x^2 - 2xy
67: p = x + y - 1
68: T = x + y
69: f = <2x^3 + 4x^2y - 2xy^2 -4\nu + 1, 4xy^2 + 2x^2y - 2y^3 -4\nu + 1>
70: Q = 3x^2 + y^2 - 2xy
72: so that
74: (1) \nabla \cdot u = 2x - 2x = 0
76: (2) u \cdot \nabla u - \nu \Delta u + \nabla p - f
77: = <2x^3 + 4x^2y -2xy^2, 4xy^2 + 2x^2y - 2y^3> -\nu <4, 4> + <1, 1> - <2x^3 + 4x^2y - 2xy^2 -4\nu + 1, 4xy^2 + 2x^2y - 2y^3 - 4\nu + 1> = 0
79: (3) u \cdot \nabla T - \alpha \Delta T - Q = 3x^2 + y^2 - 2xy - \alpha*0 - 3x^2 - y^2 + 2xy = 0
80: */
82: static PetscErrorCode quadratic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
83: {
84: u[0] = X[0] * X[0] + X[1] * X[1];
85: u[1] = 2.0 * X[0] * X[0] - 2.0 * X[0] * X[1];
86: return 0;
87: }
89: static PetscErrorCode linear_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
90: {
91: p[0] = X[0] + X[1] - 1.0;
92: return 0;
93: }
95: static PetscErrorCode linear_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
96: {
97: T[0] = X[0] + X[1];
98: return 0;
99: }
101: static void f0_quadratic_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
102: {
103: PetscInt c, d;
104: PetscInt Nc = dim;
105: const PetscReal nu = PetscRealPart(constants[0]);
107: for (c = 0; c < Nc; ++c) {
108: for (d = 0; d < dim; ++d) f0[c] += u[d] * u_x[c * dim + d];
109: }
110: f0[0] -= (2 * X[0] * X[0] * X[0] + 4 * X[0] * X[0] * X[1] - 2 * X[0] * X[1] * X[1] - 4.0 * nu + 1);
111: f0[1] -= (4 * X[0] * X[1] * X[1] + 2 * X[0] * X[0] * X[1] - 2 * X[1] * X[1] * X[1] - 4.0 * nu + 1);
112: }
114: static void f0_quadratic_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
115: {
116: PetscInt d;
117: for (d = 0, f0[0] = 0; d < dim; ++d) f0[0] += u[uOff[0] + d] * u_x[uOff_x[2] + d];
118: f0[0] -= (3 * X[0] * X[0] + X[1] * X[1] - 2 * X[0] * X[1]);
119: }
121: /*
122: CASE: cubic
123: In 2D we use exact solution:
125: u = x^3 + y^3
126: v = 2x^3 - 3x^2y
127: p = 3/2 x^2 + 3/2 y^2 - 1
128: T = 1/2 x^2 + 1/2 y^2
129: f = <3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y), 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y>
130: Q = x^4 + xy^3 + 2x^3y - 3x^2y^2 - 2
132: so that
134: \nabla \cdot u = 3x^2 - 3x^2 = 0
136: u \cdot \nabla u - \nu \Delta u + \nabla p - f
137: = <3x^5 + 6x^3y^2 - 6x^2y^3, 6x^2y^3 + 3x^4y - 6xy^4> - \nu<6x + 6y, 12x - 6y> + <3x, 3y> - <3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y), 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y> = 0
139: u \cdot \nabla T - \alpha\Delta T - Q = (x^3 + y^3) x + (2x^3 - 3x^2y) y - 2*\alpha - (x^4 + xy^3 + 2x^3y - 3x^2y^2 - 2) = 0
140: */
142: static PetscErrorCode cubic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, void *ctx)
143: {
144: u[0] = X[0] * X[0] * X[0] + X[1] * X[1] * X[1];
145: u[1] = 2.0 * X[0] * X[0] * X[0] - 3.0 * X[0] * X[0] * X[1];
146: return 0;
147: }
149: static PetscErrorCode quadratic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, void *ctx)
150: {
151: p[0] = 3.0 * X[0] * X[0] / 2.0 + 3.0 * X[1] * X[1] / 2.0 - 1.0;
152: return 0;
153: }
155: static PetscErrorCode quadratic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, void *ctx)
156: {
157: T[0] = X[0] * X[0] / 2.0 + X[1] * X[1] / 2.0;
158: return 0;
159: }
161: static void f0_cubic_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
162: {
163: PetscInt c, d;
164: PetscInt Nc = dim;
165: const PetscReal nu = PetscRealPart(constants[0]);
167: for (c = 0; c < Nc; ++c) {
168: for (d = 0; d < dim; ++d) f0[c] += u[d] * u_x[c * dim + d];
169: }
170: f0[0] -= (3 * X[0] * X[0] * X[0] * X[0] * X[0] + 6 * X[0] * X[0] * X[0] * X[1] * X[1] - 6 * X[0] * X[0] * X[1] * X[1] * X[1] - (6 * X[0] + 6 * X[1]) * nu + 3 * X[0]);
171: f0[1] -= (6 * X[0] * X[0] * X[1] * X[1] * X[1] + 3 * X[0] * X[0] * X[0] * X[0] * X[1] - 6 * X[0] * X[1] * X[1] * X[1] * X[1] - (12 * X[0] - 6 * X[1]) * nu + 3 * X[1]);
172: }
174: static void f0_cubic_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
175: {
176: const PetscReal alpha = PetscRealPart(constants[1]);
177: PetscInt d;
179: for (d = 0, f0[0] = 0; d < dim; ++d) f0[0] += u[uOff[0] + d] * u_x[uOff_x[2] + d];
180: f0[0] -= (X[0] * X[0] * X[0] * X[0] + X[0] * X[1] * X[1] * X[1] + 2.0 * X[0] * X[0] * X[0] * X[1] - 3.0 * X[0] * X[0] * X[1] * X[1] - 2.0 * alpha);
181: }
183: static void f0_q(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
184: {
185: PetscInt d;
186: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d * dim + d];
187: }
189: static void f1_v(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
190: {
191: const PetscReal nu = PetscRealPart(constants[0]);
192: const PetscInt Nc = dim;
193: PetscInt c, d;
195: for (c = 0; c < Nc; ++c) {
196: for (d = 0; d < dim; ++d) {
197: f1[c * dim + d] = nu * (u_x[c * dim + d] + u_x[d * dim + c]);
198: //f1[c*dim+d] = nu*u_x[c*dim+d];
199: }
200: f1[c * dim + c] -= u[uOff[1]];
201: }
202: }
204: static void f1_w(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal X[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
205: {
206: const PetscReal alpha = PetscRealPart(constants[1]);
207: PetscInt d;
208: for (d = 0; d < dim; ++d) f1[d] = alpha * u_x[uOff_x[2] + d];
209: }
211: /* Jacobians */
212: static void g1_qu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
213: {
214: PetscInt d;
215: for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
216: }
218: static void g0_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
219: {
220: const PetscInt Nc = dim;
221: PetscInt c, d;
223: for (c = 0; c < Nc; ++c) {
224: for (d = 0; d < dim; ++d) g0[c * Nc + d] = u_x[c * Nc + d];
225: }
226: }
228: static void g1_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
229: {
230: PetscInt NcI = dim;
231: PetscInt NcJ = dim;
232: PetscInt c, d, e;
234: for (c = 0; c < NcI; ++c) {
235: for (d = 0; d < NcJ; ++d) {
236: for (e = 0; e < dim; ++e) {
237: if (c == d) g1[(c * NcJ + d) * dim + e] = u[e];
238: }
239: }
240: }
241: }
243: static void g2_vp(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
244: {
245: PetscInt d;
246: for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0;
247: }
249: static void g3_vu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
250: {
251: const PetscReal nu = PetscRealPart(constants[0]);
252: const PetscInt Nc = dim;
253: PetscInt c, d;
255: for (c = 0; c < Nc; ++c) {
256: for (d = 0; d < dim; ++d) {
257: g3[((c * Nc + c) * dim + d) * dim + d] += nu; // gradU
258: g3[((c * Nc + d) * dim + d) * dim + c] += nu; // gradU transpose
259: }
260: }
261: }
263: static void g0_wu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
264: {
265: PetscInt d;
266: for (d = 0; d < dim; ++d) g0[d] = u_x[uOff_x[2] + d];
267: }
269: static void g1_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
270: {
271: PetscInt d;
272: for (d = 0; d < dim; ++d) g1[d] = u[uOff[0] + d];
273: }
275: static void g3_wT(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
276: {
277: const PetscReal alpha = PetscRealPart(constants[1]);
278: PetscInt d;
280: for (d = 0; d < dim; ++d) g3[d * dim + d] = alpha;
281: }
283: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
284: {
285: PetscInt sol;
288: options->solType = SOL_QUADRATIC;
289: options->showError = PETSC_FALSE;
290: PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
291: sol = options->solType;
292: PetscOptionsEList("-sol_type", "The solution type", "ex62.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL);
293: options->solType = (SolType)sol;
294: PetscOptionsBool("-show_error", "Output the error for verification", "ex62.c", options->showError, &options->showError, NULL);
295: PetscOptionsEnd();
296: return 0;
297: }
299: static PetscErrorCode SetupParameters(AppCtx *user)
300: {
301: PetscBag bag;
302: Parameter *p;
305: /* setup PETSc parameter bag */
306: PetscBagGetData(user->bag, (void **)&p);
307: PetscBagSetName(user->bag, "par", "Poiseuille flow parameters");
308: bag = user->bag;
309: PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity");
310: PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity");
311: PetscBagRegisterReal(bag, &p->theta, 0.0, "theta", "Angle of pipe wall to x-axis");
312: return 0;
313: }
315: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
316: {
318: DMCreate(comm, dm);
319: DMSetType(*dm, DMPLEX);
320: DMSetFromOptions(*dm);
321: {
322: Parameter *param;
323: Vec coordinates;
324: PetscScalar *coords;
325: PetscReal theta;
326: PetscInt cdim, N, bs, i;
328: DMGetCoordinateDim(*dm, &cdim);
329: DMGetCoordinates(*dm, &coordinates);
330: VecGetLocalSize(coordinates, &N);
331: VecGetBlockSize(coordinates, &bs);
333: VecGetArray(coordinates, &coords);
334: PetscBagGetData(user->bag, (void **)¶m);
335: theta = param->theta;
336: for (i = 0; i < N; i += cdim) {
337: PetscScalar x = coords[i + 0];
338: PetscScalar y = coords[i + 1];
340: coords[i + 0] = PetscCosReal(theta) * x - PetscSinReal(theta) * y;
341: coords[i + 1] = PetscSinReal(theta) * x + PetscCosReal(theta) * y;
342: }
343: VecRestoreArray(coordinates, &coords);
344: DMSetCoordinates(*dm, coordinates);
345: }
346: DMViewFromOptions(*dm, NULL, "-dm_view");
347: return 0;
348: }
350: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
351: {
352: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
353: PetscDS prob;
354: DMLabel label;
355: Parameter *ctx;
356: PetscInt id;
359: DMGetLabel(dm, "marker", &label);
360: DMGetDS(dm, &prob);
361: switch (user->solType) {
362: case SOL_QUADRATIC:
363: PetscDSSetResidual(prob, 0, f0_quadratic_v, f1_v);
364: PetscDSSetResidual(prob, 2, f0_quadratic_w, f1_w);
366: exactFuncs[0] = quadratic_u;
367: exactFuncs[1] = linear_p;
368: exactFuncs[2] = linear_T;
369: break;
370: case SOL_CUBIC:
371: PetscDSSetResidual(prob, 0, f0_cubic_v, f1_v);
372: PetscDSSetResidual(prob, 2, f0_cubic_w, f1_w);
374: exactFuncs[0] = cubic_u;
375: exactFuncs[1] = quadratic_p;
376: exactFuncs[2] = quadratic_T;
377: break;
378: default:
379: SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType);
380: }
382: PetscDSSetResidual(prob, 1, f0_q, NULL);
384: PetscDSSetJacobian(prob, 0, 0, g0_vu, g1_vu, NULL, g3_vu);
385: PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_vp, NULL);
386: PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL);
387: PetscDSSetJacobian(prob, 2, 0, g0_wu, NULL, NULL, NULL);
388: PetscDSSetJacobian(prob, 2, 2, NULL, g1_wT, NULL, g3_wT);
389: /* Setup constants */
390: {
391: Parameter *param;
392: PetscScalar constants[3];
394: PetscBagGetData(user->bag, (void **)¶m);
396: constants[0] = param->nu;
397: constants[1] = param->alpha;
398: constants[2] = param->theta;
399: PetscDSSetConstants(prob, 3, constants);
400: }
401: /* Setup Boundary Conditions */
402: PetscBagGetData(user->bag, (void **)&ctx);
403: id = 3;
404: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void))exactFuncs[0], NULL, ctx, NULL);
405: id = 1;
406: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void))exactFuncs[0], NULL, ctx, NULL);
407: id = 2;
408: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void))exactFuncs[0], NULL, ctx, NULL);
409: id = 4;
410: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall velocity", label, 1, &id, 0, 0, NULL, (void (*)(void))exactFuncs[0], NULL, ctx, NULL);
411: id = 3;
412: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void))exactFuncs[2], NULL, ctx, NULL);
413: id = 1;
414: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void))exactFuncs[2], NULL, ctx, NULL);
415: id = 2;
416: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void))exactFuncs[2], NULL, ctx, NULL);
417: id = 4;
418: PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall temp", label, 1, &id, 2, 0, NULL, (void (*)(void))exactFuncs[2], NULL, ctx, NULL);
420: /*setup exact solution.*/
421: PetscDSSetExactSolution(prob, 0, exactFuncs[0], ctx);
422: PetscDSSetExactSolution(prob, 1, exactFuncs[1], ctx);
423: PetscDSSetExactSolution(prob, 2, exactFuncs[2], ctx);
424: return 0;
425: }
427: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
428: {
429: DM cdm = dm;
430: PetscFE fe[3];
431: Parameter *param;
432: MPI_Comm comm;
433: PetscInt dim;
434: PetscBool simplex;
437: DMGetDimension(dm, &dim);
438: DMPlexIsSimplex(dm, &simplex);
439: /* Create finite element */
440: PetscObjectGetComm((PetscObject)dm, &comm);
441: PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]);
442: PetscObjectSetName((PetscObject)fe[0], "velocity");
444: PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]);
445: PetscFECopyQuadrature(fe[0], fe[1]);
446: PetscObjectSetName((PetscObject)fe[1], "pressure");
448: PetscFECreateDefault(comm, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]);
449: PetscFECopyQuadrature(fe[0], fe[2]);
450: PetscObjectSetName((PetscObject)fe[2], "temperature");
452: /* Set discretization and boundary conditions for each mesh */
453: DMSetField(dm, 0, NULL, (PetscObject)fe[0]);
454: DMSetField(dm, 1, NULL, (PetscObject)fe[1]);
455: DMSetField(dm, 2, NULL, (PetscObject)fe[2]);
456: DMCreateDS(dm);
457: SetupProblem(dm, user);
458: PetscBagGetData(user->bag, (void **)¶m);
459: while (cdm) {
460: DMCopyDisc(dm, cdm);
461: DMPlexCreateBasisRotation(cdm, param->theta, 0.0, 0.0);
462: DMGetCoarseDM(cdm, &cdm);
463: }
464: PetscFEDestroy(&fe[0]);
465: PetscFEDestroy(&fe[1]);
466: PetscFEDestroy(&fe[2]);
467: return 0;
468: }
470: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
471: {
472: Vec vec;
473: PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zero, zero, zero};
477: funcs[nfield] = constant;
478: DMCreateGlobalVector(dm, &vec);
479: DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec);
480: VecNormalize(vec, NULL);
481: PetscObjectSetName((PetscObject)vec, "Pressure Null Space");
482: VecViewFromOptions(vec, NULL, "-pressure_nullspace_view");
483: MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace);
484: VecDestroy(&vec);
485: return 0;
486: }
488: int main(int argc, char **argv)
489: {
490: SNES snes; /* nonlinear solver */
491: DM dm; /* problem definition */
492: Vec u, r; /* solution, residual vectors */
493: AppCtx user; /* user-defined work context */
496: PetscInitialize(&argc, &argv, NULL, help);
497: ProcessOptions(PETSC_COMM_WORLD, &user);
498: PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag);
499: SetupParameters(&user);
500: SNESCreate(PETSC_COMM_WORLD, &snes);
501: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
502: SNESSetDM(snes, dm);
503: DMSetApplicationContext(dm, &user);
504: /* Setup problem */
505: SetupDiscretization(dm, &user);
506: DMPlexCreateClosureIndex(dm, NULL);
508: DMCreateGlobalVector(dm, &u);
509: PetscObjectSetName((PetscObject)u, "Solution");
510: VecDuplicate(u, &r);
512: DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace);
513: DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
515: SNESSetFromOptions(snes);
516: {
517: PetscDS ds;
518: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
519: void *ctxs[3];
521: DMGetDS(dm, &ds);
522: PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]);
523: PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]);
524: PetscDSGetExactSolution(ds, 2, &exactFuncs[2], &ctxs[2]);
525: DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, u);
526: PetscObjectSetName((PetscObject)u, "Exact Solution");
527: VecViewFromOptions(u, NULL, "-exact_vec_view");
528: }
529: DMSNESCheckFromOptions(snes, u);
530: VecSet(u, 0.0);
531: SNESSolve(snes, NULL, u);
533: if (user.showError) {
534: PetscDS ds;
535: Vec r;
536: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
537: void *ctxs[3];
539: DMGetDS(dm, &ds);
540: PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]);
541: PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]);
542: PetscDSGetExactSolution(ds, 2, &exactFuncs[2], &ctxs[2]);
543: DMGetGlobalVector(dm, &r);
544: DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, r);
545: VecAXPY(r, -1.0, u);
546: PetscObjectSetName((PetscObject)r, "Solution Error");
547: VecViewFromOptions(r, NULL, "-error_vec_view");
548: DMRestoreGlobalVector(dm, &r);
549: }
550: PetscObjectSetName((PetscObject)u, "Numerical Solution");
551: VecViewFromOptions(u, NULL, "-sol_vec_view");
553: VecDestroy(&u);
554: VecDestroy(&r);
555: DMDestroy(&dm);
556: SNESDestroy(&snes);
557: PetscBagDestroy(&user.bag);
558: PetscFinalize();
559: return 0;
560: }
562: /*TEST
564: test:
565: suffix: 2d_tri_p2_p1_p1
566: requires: triangle !single
567: args: -dm_plex_separate_marker -sol_type quadratic -dm_refine 0 \
568: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
569: -dmsnes_check .001 -snes_error_if_not_converged \
570: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
571: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
572: -fieldsplit_0_pc_type lu \
573: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
575: test:
576: # Using -dm_refine 2 -convest_num_refine 3 gives L_2 convergence rate: [2.9, 2.3, 1.9]
577: suffix: 2d_tri_p2_p1_p1_conv
578: requires: triangle !single
579: args: -dm_plex_separate_marker -sol_type cubic -dm_refine 0 \
580: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
581: -snes_error_if_not_converged -snes_convergence_test correct_pressure -snes_convergence_estimate -convest_num_refine 1 \
582: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
583: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
584: -fieldsplit_0_pc_type lu \
585: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
587: test:
588: suffix: 2d_tri_p3_p2_p2
589: requires: triangle !single
590: args: -dm_plex_separate_marker -sol_type cubic -dm_refine 0 \
591: -vel_petscspace_degree 3 -pres_petscspace_degree 2 -temp_petscspace_degree 2 \
592: -dmsnes_check .001 -snes_error_if_not_converged \
593: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
594: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
595: -fieldsplit_0_pc_type lu \
596: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
598: TEST*/