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, PetscCtx ctx)
48: {
49: for (PetscInt d = 0; d < Nc; ++d) u[d] = 0.0;
50: return PETSC_SUCCESS;
51: }
53: static PetscErrorCode constant(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
54: {
55: for (PetscInt d = 0; d < Nc; ++d) u[d] = 1.0;
56: return PETSC_SUCCESS;
57: }
59: /*
60: CASE: quadratic
61: In 2D we use exact solution:
63: u = x^2 + y^2
64: v = 2x^2 - 2xy
65: p = x + y - 1
66: T = x + y
67: f = <2x^3 + 4x^2y - 2xy^2 -4\nu + 1, 4xy^2 + 2x^2y - 2y^3 -4\nu + 1>
68: Q = 3x^2 + y^2 - 2xy
70: so that
72: (1) \nabla \cdot u = 2x - 2x = 0
74: (2) u \cdot \nabla u - \nu \Delta u + \nabla p - f
75: = <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
77: (3) u \cdot \nabla T - \alpha \Delta T - Q = 3x^2 + y^2 - 2xy - \alpha*0 - 3x^2 - y^2 + 2xy = 0
78: */
80: static PetscErrorCode quadratic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
81: {
82: u[0] = X[0] * X[0] + X[1] * X[1];
83: u[1] = 2.0 * X[0] * X[0] - 2.0 * X[0] * X[1];
84: return PETSC_SUCCESS;
85: }
87: static PetscErrorCode linear_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
88: {
89: p[0] = X[0] + X[1] - 1.0;
90: return PETSC_SUCCESS;
91: }
93: static PetscErrorCode linear_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
94: {
95: T[0] = X[0] + X[1];
96: return PETSC_SUCCESS;
97: }
99: 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[])
100: {
101: PetscInt Nc = dim;
102: const PetscReal nu = PetscRealPart(constants[0]);
104: for (PetscInt c = 0; c < Nc; ++c) {
105: for (PetscInt d = 0; d < dim; ++d) f0[c] += u[d] * u_x[c * dim + d];
106: }
107: 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);
108: 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);
109: }
111: 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[])
112: {
113: f0[0] = 0.0;
114: for (PetscInt d = 0; d < dim; ++d) f0[0] += u[uOff[0] + d] * u_x[uOff_x[2] + d];
115: f0[0] -= (3 * X[0] * X[0] + X[1] * X[1] - 2 * X[0] * X[1]);
116: }
118: /*
119: CASE: cubic
120: In 2D we use exact solution:
122: u = x^3 + y^3
123: v = 2x^3 - 3x^2y
124: p = 3/2 x^2 + 3/2 y^2 - 1
125: T = 1/2 x^2 + 1/2 y^2
126: f = <3x^5 + 6x^3y^2 - 6x^2y^3 - \nu(6x + 6y), 6x^2y^3 + 3x^4y - 6xy^4 - \nu(12x - 6y) + 3y>
127: Q = x^4 + xy^3 + 2x^3y - 3x^2y^2 - 2
129: so that
131: \nabla \cdot u = 3x^2 - 3x^2 = 0
133: u \cdot \nabla u - \nu \Delta u + \nabla p - f
134: = <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
136: 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
137: */
139: static PetscErrorCode cubic_u(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *u, PetscCtx ctx)
140: {
141: u[0] = X[0] * X[0] * X[0] + X[1] * X[1] * X[1];
142: u[1] = 2.0 * X[0] * X[0] * X[0] - 3.0 * X[0] * X[0] * X[1];
143: return PETSC_SUCCESS;
144: }
146: static PetscErrorCode quadratic_p(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *p, PetscCtx ctx)
147: {
148: p[0] = 3.0 * X[0] * X[0] / 2.0 + 3.0 * X[1] * X[1] / 2.0 - 1.0;
149: return PETSC_SUCCESS;
150: }
152: static PetscErrorCode quadratic_T(PetscInt Dim, PetscReal time, const PetscReal X[], PetscInt Nf, PetscScalar *T, PetscCtx ctx)
153: {
154: T[0] = X[0] * X[0] / 2.0 + X[1] * X[1] / 2.0;
155: return PETSC_SUCCESS;
156: }
158: 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[])
159: {
160: PetscInt Nc = dim;
161: const PetscReal nu = PetscRealPart(constants[0]);
163: for (PetscInt c = 0; c < Nc; ++c) {
164: for (PetscInt d = 0; d < dim; ++d) f0[c] += u[d] * u_x[c * dim + d];
165: }
166: 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]);
167: 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]);
168: }
170: 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[])
171: {
172: const PetscReal alpha = PetscRealPart(constants[1]);
174: f0[0] = 0.0;
175: for (PetscInt d = 0; d < dim; ++d) f0[0] += u[uOff[0] + d] * u_x[uOff_x[2] + d];
176: 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);
177: }
179: 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[])
180: {
181: f0[0] = 0.0;
182: for (PetscInt d = 0; d < dim; ++d) f0[0] += u_x[d * dim + d];
183: }
185: 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[])
186: {
187: const PetscReal nu = PetscRealPart(constants[0]);
188: const PetscInt Nc = dim;
190: for (PetscInt c = 0; c < Nc; ++c) {
191: for (PetscInt d = 0; d < dim; ++d) {
192: f1[c * dim + d] = nu * (u_x[c * dim + d] + u_x[d * dim + c]);
193: //f1[c*dim+d] = nu*u_x[c*dim+d];
194: }
195: f1[c * dim + c] -= u[uOff[1]];
196: }
197: }
199: 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[])
200: {
201: const PetscReal alpha = PetscRealPart(constants[1]);
203: for (PetscInt d = 0; d < dim; ++d) f1[d] = alpha * u_x[uOff_x[2] + d];
204: }
206: /* Jacobians */
207: 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[])
208: {
209: for (PetscInt d = 0; d < dim; ++d) g1[d * dim + d] = 1.0;
210: }
212: 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[])
213: {
214: const PetscInt Nc = dim;
216: for (PetscInt c = 0; c < Nc; ++c) {
217: for (PetscInt d = 0; d < dim; ++d) g0[c * Nc + d] = u_x[c * Nc + d];
218: }
219: }
221: 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[])
222: {
223: PetscInt NcI = dim;
224: PetscInt NcJ = dim;
225: PetscInt c, d, e;
227: for (c = 0; c < NcI; ++c) {
228: for (d = 0; d < NcJ; ++d) {
229: for (e = 0; e < dim; ++e) {
230: if (c == d) g1[(c * NcJ + d) * dim + e] = u[e];
231: }
232: }
233: }
234: }
236: 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[])
237: {
238: for (PetscInt d = 0; d < dim; ++d) g2[d * dim + d] = -1.0;
239: }
241: 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[])
242: {
243: const PetscReal nu = PetscRealPart(constants[0]);
244: const PetscInt Nc = dim;
246: for (PetscInt c = 0; c < Nc; ++c) {
247: for (PetscInt d = 0; d < dim; ++d) {
248: g3[((c * Nc + c) * dim + d) * dim + d] += nu; // gradU
249: g3[((c * Nc + d) * dim + d) * dim + c] += nu; // gradU transpose
250: }
251: }
252: }
254: 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[])
255: {
256: for (PetscInt d = 0; d < dim; ++d) g0[d] = u_x[uOff_x[2] + d];
257: }
259: 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[])
260: {
261: for (PetscInt d = 0; d < dim; ++d) g1[d] = u[uOff[0] + d];
262: }
264: 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[])
265: {
266: const PetscReal alpha = PetscRealPart(constants[1]);
268: for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = alpha;
269: }
271: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
272: {
273: PetscInt sol;
275: PetscFunctionBeginUser;
276: options->solType = SOL_QUADRATIC;
277: options->showError = PETSC_FALSE;
278: PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
279: sol = options->solType;
280: PetscCall(PetscOptionsEList("-sol_type", "The solution type", "ex62.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL));
281: options->solType = (SolType)sol;
282: PetscCall(PetscOptionsBool("-show_error", "Output the error for verification", "ex62.c", options->showError, &options->showError, NULL));
283: PetscOptionsEnd();
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode SetupParameters(AppCtx *user)
288: {
289: PetscBag bag;
290: Parameter *p;
292: PetscFunctionBeginUser;
293: /* setup PETSc parameter bag */
294: PetscCall(PetscBagGetData(user->bag, &p));
295: PetscCall(PetscBagSetName(user->bag, "par", "Poiseuille flow parameters"));
296: bag = user->bag;
297: PetscCall(PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity"));
298: PetscCall(PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity"));
299: PetscCall(PetscBagRegisterReal(bag, &p->theta, 0.0, "theta", "Angle of pipe wall to x-axis"));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
304: {
305: PetscFunctionBeginUser;
306: PetscCall(DMCreate(comm, dm));
307: PetscCall(DMSetType(*dm, DMPLEX));
308: PetscCall(DMSetFromOptions(*dm));
309: {
310: Parameter *param;
311: Vec coordinates;
312: PetscScalar *coords;
313: PetscReal theta;
314: PetscInt cdim, N, bs, i;
316: PetscCall(DMGetCoordinateDim(*dm, &cdim));
317: PetscCall(DMGetCoordinates(*dm, &coordinates));
318: PetscCall(VecGetLocalSize(coordinates, &N));
319: PetscCall(VecGetBlockSize(coordinates, &bs));
320: PetscCheck(bs == cdim, comm, PETSC_ERR_ARG_WRONG, "Invalid coordinate blocksize %" PetscInt_FMT " != embedding dimension %" PetscInt_FMT, bs, cdim);
321: PetscCall(VecGetArray(coordinates, &coords));
322: PetscCall(PetscBagGetData(user->bag, ¶m));
323: theta = param->theta;
324: for (i = 0; i < N; i += cdim) {
325: PetscScalar x = coords[i + 0];
326: PetscScalar y = coords[i + 1];
328: coords[i + 0] = PetscCosReal(theta) * x - PetscSinReal(theta) * y;
329: coords[i + 1] = PetscSinReal(theta) * x + PetscCosReal(theta) * y;
330: }
331: PetscCall(VecRestoreArray(coordinates, &coords));
332: PetscCall(DMSetCoordinates(*dm, coordinates));
333: }
334: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: static PetscErrorCode SetupProblem(DM dm, AppCtx *user)
339: {
340: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
341: PetscDS prob;
342: DMLabel label;
343: Parameter *ctx;
344: PetscInt id;
346: PetscFunctionBeginUser;
347: PetscCall(DMGetLabel(dm, "marker", &label));
348: PetscCall(DMGetDS(dm, &prob));
349: switch (user->solType) {
350: case SOL_QUADRATIC:
351: PetscCall(PetscDSSetResidual(prob, 0, f0_quadratic_v, f1_v));
352: PetscCall(PetscDSSetResidual(prob, 2, f0_quadratic_w, f1_w));
354: exactFuncs[0] = quadratic_u;
355: exactFuncs[1] = linear_p;
356: exactFuncs[2] = linear_T;
357: break;
358: case SOL_CUBIC:
359: PetscCall(PetscDSSetResidual(prob, 0, f0_cubic_v, f1_v));
360: PetscCall(PetscDSSetResidual(prob, 2, f0_cubic_w, f1_w));
362: exactFuncs[0] = cubic_u;
363: exactFuncs[1] = quadratic_p;
364: exactFuncs[2] = quadratic_T;
365: break;
366: default:
367: SETERRQ(PetscObjectComm((PetscObject)prob), PETSC_ERR_ARG_WRONG, "Unsupported solution type: %s (%d)", solTypes[PetscMin(user->solType, NUM_SOL_TYPES)], user->solType);
368: }
370: PetscCall(PetscDSSetResidual(prob, 1, f0_q, NULL));
372: PetscCall(PetscDSSetJacobian(prob, 0, 0, g0_vu, g1_vu, NULL, g3_vu));
373: PetscCall(PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_vp, NULL));
374: PetscCall(PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL));
375: PetscCall(PetscDSSetJacobian(prob, 2, 0, g0_wu, NULL, NULL, NULL));
376: PetscCall(PetscDSSetJacobian(prob, 2, 2, NULL, g1_wT, NULL, g3_wT));
377: /* Setup constants */
378: {
379: Parameter *param;
380: PetscScalar constants[3];
382: PetscCall(PetscBagGetData(user->bag, ¶m));
384: constants[0] = param->nu;
385: constants[1] = param->alpha;
386: constants[2] = param->theta;
387: PetscCall(PetscDSSetConstants(prob, 3, constants));
388: }
389: /* Setup Boundary Conditions */
390: PetscCall(PetscBagGetData(user->bag, &ctx));
391: id = 3;
392: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall velocity", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], NULL, ctx, NULL));
393: id = 1;
394: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall velocity", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], NULL, ctx, NULL));
395: id = 2;
396: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall velocity", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], NULL, ctx, NULL));
397: id = 4;
398: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall velocity", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)exactFuncs[0], NULL, ctx, NULL));
399: id = 3;
400: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "top wall temp", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exactFuncs[2], NULL, ctx, NULL));
401: id = 1;
402: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "bottom wall temp", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exactFuncs[2], NULL, ctx, NULL));
403: id = 2;
404: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "right wall temp", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exactFuncs[2], NULL, ctx, NULL));
405: id = 4;
406: PetscCall(PetscDSAddBoundary(prob, DM_BC_ESSENTIAL, "left wall temp", label, 1, &id, 2, 0, NULL, (PetscVoidFn *)exactFuncs[2], NULL, ctx, NULL));
408: /*setup exact solution.*/
409: PetscCall(PetscDSSetExactSolution(prob, 0, exactFuncs[0], ctx));
410: PetscCall(PetscDSSetExactSolution(prob, 1, exactFuncs[1], ctx));
411: PetscCall(PetscDSSetExactSolution(prob, 2, exactFuncs[2], ctx));
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
416: {
417: DM cdm = dm;
418: PetscFE fe[3];
419: Parameter *param;
420: MPI_Comm comm;
421: PetscInt dim;
422: PetscBool simplex;
424: PetscFunctionBeginUser;
425: PetscCall(DMGetDimension(dm, &dim));
426: PetscCall(DMPlexIsSimplex(dm, &simplex));
427: /* Create finite element */
428: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
429: PetscCall(PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]));
430: PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity"));
432: PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]));
433: PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
434: PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure"));
436: PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]));
437: PetscCall(PetscFECopyQuadrature(fe[0], fe[2]));
438: PetscCall(PetscObjectSetName((PetscObject)fe[2], "temperature"));
440: /* Set discretization and boundary conditions for each mesh */
441: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
442: PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1]));
443: PetscCall(DMSetField(dm, 2, NULL, (PetscObject)fe[2]));
444: PetscCall(DMCreateDS(dm));
445: PetscCall(SetupProblem(dm, user));
446: PetscCall(PetscBagGetData(user->bag, ¶m));
447: while (cdm) {
448: PetscCall(DMCopyDisc(dm, cdm));
449: PetscCall(DMPlexCreateBasisRotation(cdm, param->theta, 0.0, 0.0));
450: PetscCall(DMGetCoarseDM(cdm, &cdm));
451: }
452: PetscCall(PetscFEDestroy(&fe[0]));
453: PetscCall(PetscFEDestroy(&fe[1]));
454: PetscCall(PetscFEDestroy(&fe[2]));
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: static PetscErrorCode CreatePressureNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace)
459: {
460: Vec vec;
461: PetscErrorCode (*funcs[3])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zero, zero, zero};
463: PetscFunctionBeginUser;
464: PetscCheck(ofield == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Nullspace must be for pressure field at index 1, not %" PetscInt_FMT, ofield);
465: funcs[nfield] = constant;
466: PetscCall(DMCreateGlobalVector(dm, &vec));
467: PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
468: PetscCall(VecNormalize(vec, NULL));
469: PetscCall(PetscObjectSetName((PetscObject)vec, "Pressure Null Space"));
470: PetscCall(VecViewFromOptions(vec, NULL, "-pressure_nullspace_view"));
471: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
472: PetscCall(VecDestroy(&vec));
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: int main(int argc, char **argv)
477: {
478: SNES snes; /* nonlinear solver */
479: DM dm; /* problem definition */
480: Vec u, r; /* solution, residual vectors */
481: AppCtx user; /* user-defined work context */
483: PetscFunctionBeginUser;
484: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
485: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
486: PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag));
487: PetscCall(SetupParameters(&user));
488: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
489: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
490: PetscCall(SNESSetDM(snes, dm));
491: PetscCall(DMSetApplicationContext(dm, &user));
492: /* Setup problem */
493: PetscCall(SetupDiscretization(dm, &user));
494: PetscCall(DMPlexCreateClosureIndex(dm, NULL));
496: PetscCall(DMCreateGlobalVector(dm, &u));
497: PetscCall(PetscObjectSetName((PetscObject)u, "Solution"));
498: PetscCall(VecDuplicate(u, &r));
500: PetscCall(DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace));
501: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
503: PetscCall(SNESSetFromOptions(snes));
504: {
505: PetscDS ds;
506: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
507: PetscCtx ctxs[3];
509: PetscCall(DMGetDS(dm, &ds));
510: PetscCall(PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]));
511: PetscCall(PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]));
512: PetscCall(PetscDSGetExactSolution(ds, 2, &exactFuncs[2], &ctxs[2]));
513: PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, u));
514: PetscCall(PetscObjectSetName((PetscObject)u, "Exact Solution"));
515: PetscCall(VecViewFromOptions(u, NULL, "-exact_vec_view"));
516: }
517: PetscCall(DMSNESCheckFromOptions(snes, u));
518: PetscCall(VecSet(u, 0.0));
519: PetscCall(SNESSolve(snes, NULL, u));
521: if (user.showError) {
522: PetscDS ds;
523: Vec r;
524: PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
525: PetscCtx ctxs[3];
527: PetscCall(DMGetDS(dm, &ds));
528: PetscCall(PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]));
529: PetscCall(PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]));
530: PetscCall(PetscDSGetExactSolution(ds, 2, &exactFuncs[2], &ctxs[2]));
531: PetscCall(DMGetGlobalVector(dm, &r));
532: PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, r));
533: PetscCall(VecAXPY(r, -1.0, u));
534: PetscCall(PetscObjectSetName((PetscObject)r, "Solution Error"));
535: PetscCall(VecViewFromOptions(r, NULL, "-error_vec_view"));
536: PetscCall(DMRestoreGlobalVector(dm, &r));
537: }
538: PetscCall(PetscObjectSetName((PetscObject)u, "Numerical Solution"));
539: PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));
541: PetscCall(VecDestroy(&u));
542: PetscCall(VecDestroy(&r));
543: PetscCall(DMDestroy(&dm));
544: PetscCall(SNESDestroy(&snes));
545: PetscCall(PetscBagDestroy(&user.bag));
546: PetscCall(PetscFinalize());
547: return 0;
548: }
550: /*TEST
552: test:
553: suffix: 2d_tri_p2_p1_p1
554: requires: triangle !single
555: args: -dm_plex_separate_marker -sol_type quadratic -dm_refine 0 \
556: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
557: -dmsnes_check .001 -snes_error_if_not_converged \
558: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
559: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
560: -fieldsplit_0_pc_type lu \
561: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
563: test:
564: # Using -dm_refine 2 -convest_num_refine 3 gives L_2 convergence rate: [2.9, 2.3, 1.9]
565: suffix: 2d_tri_p2_p1_p1_conv
566: requires: triangle !single
567: args: -dm_plex_separate_marker -sol_type cubic -dm_refine 0 \
568: -vel_petscspace_degree 2 -pres_petscspace_degree 1 -temp_petscspace_degree 1 \
569: -snes_error_if_not_converged -snes_convergence_test correct_pressure -snes_convergence_estimate -convest_num_refine 1 \
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: suffix: 2d_tri_p3_p2_p2
577: requires: triangle !single
578: args: -dm_plex_separate_marker -sol_type cubic -dm_refine 0 \
579: -vel_petscspace_degree 3 -pres_petscspace_degree 2 -temp_petscspace_degree 2 \
580: -dmsnes_check .001 -snes_error_if_not_converged \
581: -ksp_type fgmres -ksp_gmres_restart 10 -ksp_rtol 1.0e-9 -ksp_error_if_not_converged \
582: -pc_type fieldsplit -pc_fieldsplit_0_fields 0,2 -pc_fieldsplit_1_fields 1 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full \
583: -fieldsplit_0_pc_type lu \
584: -fieldsplit_pressure_ksp_rtol 1e-10 -fieldsplit_pressure_pc_type jacobi
586: TEST*/