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 PETSC_SUCCESS;
 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 PETSC_SUCCESS;
 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 PETSC_SUCCESS;
 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 PETSC_SUCCESS;
 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 PETSC_SUCCESS;
 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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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 PETSC_SUCCESS;
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;

287:   PetscFunctionBeginUser;
288:   options->solType   = SOL_QUADRATIC;
289:   options->showError = PETSC_FALSE;
290:   PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
291:   sol = options->solType;
292:   PetscCall(PetscOptionsEList("-sol_type", "The solution type", "ex62.c", solTypes, NUM_SOL_TYPES, solTypes[options->solType], &sol, NULL));
293:   options->solType = (SolType)sol;
294:   PetscCall(PetscOptionsBool("-show_error", "Output the error for verification", "ex62.c", options->showError, &options->showError, NULL));
295:   PetscOptionsEnd();
296:   PetscFunctionReturn(PETSC_SUCCESS);
297: }

299: static PetscErrorCode SetupParameters(AppCtx *user)
300: {
301:   PetscBag   bag;
302:   Parameter *p;

304:   PetscFunctionBeginUser;
305:   /* setup PETSc parameter bag */
306:   PetscCall(PetscBagGetData(user->bag, (void **)&p));
307:   PetscCall(PetscBagSetName(user->bag, "par", "Poiseuille flow parameters"));
308:   bag = user->bag;
309:   PetscCall(PetscBagRegisterReal(bag, &p->nu, 1.0, "nu", "Kinematic viscosity"));
310:   PetscCall(PetscBagRegisterReal(bag, &p->alpha, 1.0, "alpha", "Thermal diffusivity"));
311:   PetscCall(PetscBagRegisterReal(bag, &p->theta, 0.0, "theta", "Angle of pipe wall to x-axis"));
312:   PetscFunctionReturn(PETSC_SUCCESS);
313: }

315: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
316: {
317:   PetscFunctionBeginUser;
318:   PetscCall(DMCreate(comm, dm));
319:   PetscCall(DMSetType(*dm, DMPLEX));
320:   PetscCall(DMSetFromOptions(*dm));
321:   {
322:     Parameter   *param;
323:     Vec          coordinates;
324:     PetscScalar *coords;
325:     PetscReal    theta;
326:     PetscInt     cdim, N, bs, i;

328:     PetscCall(DMGetCoordinateDim(*dm, &cdim));
329:     PetscCall(DMGetCoordinates(*dm, &coordinates));
330:     PetscCall(VecGetLocalSize(coordinates, &N));
331:     PetscCall(VecGetBlockSize(coordinates, &bs));
332:     PetscCheck(bs == cdim, comm, PETSC_ERR_ARG_WRONG, "Invalid coordinate blocksize %" PetscInt_FMT " != embedding dimension %" PetscInt_FMT, bs, cdim);
333:     PetscCall(VecGetArray(coordinates, &coords));
334:     PetscCall(PetscBagGetData(user->bag, (void **)&param));
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:     PetscCall(VecRestoreArray(coordinates, &coords));
344:     PetscCall(DMSetCoordinates(*dm, coordinates));
345:   }
346:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
347:   PetscFunctionReturn(PETSC_SUCCESS);
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;

358:   PetscFunctionBeginUser;
359:   PetscCall(DMGetLabel(dm, "marker", &label));
360:   PetscCall(DMGetDS(dm, &prob));
361:   switch (user->solType) {
362:   case SOL_QUADRATIC:
363:     PetscCall(PetscDSSetResidual(prob, 0, f0_quadratic_v, f1_v));
364:     PetscCall(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:     PetscCall(PetscDSSetResidual(prob, 0, f0_cubic_v, f1_v));
372:     PetscCall(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:   PetscCall(PetscDSSetResidual(prob, 1, f0_q, NULL));

384:   PetscCall(PetscDSSetJacobian(prob, 0, 0, g0_vu, g1_vu, NULL, g3_vu));
385:   PetscCall(PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_vp, NULL));
386:   PetscCall(PetscDSSetJacobian(prob, 1, 0, NULL, g1_qu, NULL, NULL));
387:   PetscCall(PetscDSSetJacobian(prob, 2, 0, g0_wu, NULL, NULL, NULL));
388:   PetscCall(PetscDSSetJacobian(prob, 2, 2, NULL, g1_wT, NULL, g3_wT));
389:   /* Setup constants */
390:   {
391:     Parameter  *param;
392:     PetscScalar constants[3];

394:     PetscCall(PetscBagGetData(user->bag, (void **)&param));

396:     constants[0] = param->nu;
397:     constants[1] = param->alpha;
398:     constants[2] = param->theta;
399:     PetscCall(PetscDSSetConstants(prob, 3, constants));
400:   }
401:   /* Setup Boundary Conditions */
402:   PetscCall(PetscBagGetData(user->bag, (void **)&ctx));
403:   id = 3;
404:   PetscCall(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:   PetscCall(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:   PetscCall(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:   PetscCall(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:   PetscCall(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:   PetscCall(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:   PetscCall(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:   PetscCall(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:   PetscCall(PetscDSSetExactSolution(prob, 0, exactFuncs[0], ctx));
422:   PetscCall(PetscDSSetExactSolution(prob, 1, exactFuncs[1], ctx));
423:   PetscCall(PetscDSSetExactSolution(prob, 2, exactFuncs[2], ctx));
424:   PetscFunctionReturn(PETSC_SUCCESS);
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;

436:   PetscFunctionBeginUser;
437:   PetscCall(DMGetDimension(dm, &dim));
438:   PetscCall(DMPlexIsSimplex(dm, &simplex));
439:   /* Create finite element */
440:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
441:   PetscCall(PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0]));
442:   PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity"));

444:   PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1]));
445:   PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
446:   PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure"));

448:   PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "temp_", PETSC_DEFAULT, &fe[2]));
449:   PetscCall(PetscFECopyQuadrature(fe[0], fe[2]));
450:   PetscCall(PetscObjectSetName((PetscObject)fe[2], "temperature"));

452:   /* Set discretization and boundary conditions for each mesh */
453:   PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0]));
454:   PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1]));
455:   PetscCall(DMSetField(dm, 2, NULL, (PetscObject)fe[2]));
456:   PetscCall(DMCreateDS(dm));
457:   PetscCall(SetupProblem(dm, user));
458:   PetscCall(PetscBagGetData(user->bag, (void **)&param));
459:   while (cdm) {
460:     PetscCall(DMCopyDisc(dm, cdm));
461:     PetscCall(DMPlexCreateBasisRotation(cdm, param->theta, 0.0, 0.0));
462:     PetscCall(DMGetCoarseDM(cdm, &cdm));
463:   }
464:   PetscCall(PetscFEDestroy(&fe[0]));
465:   PetscCall(PetscFEDestroy(&fe[1]));
466:   PetscCall(PetscFEDestroy(&fe[2]));
467:   PetscFunctionReturn(PETSC_SUCCESS);
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};

475:   PetscFunctionBeginUser;
476:   PetscCheck(ofield == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Nullspace must be for pressure field at index 1, not %" PetscInt_FMT, ofield);
477:   funcs[nfield] = constant;
478:   PetscCall(DMCreateGlobalVector(dm, &vec));
479:   PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
480:   PetscCall(VecNormalize(vec, NULL));
481:   PetscCall(PetscObjectSetName((PetscObject)vec, "Pressure Null Space"));
482:   PetscCall(VecViewFromOptions(vec, NULL, "-pressure_nullspace_view"));
483:   PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
484:   PetscCall(VecDestroy(&vec));
485:   PetscFunctionReturn(PETSC_SUCCESS);
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 */

495:   PetscFunctionBeginUser;
496:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
497:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
498:   PetscCall(PetscBagCreate(PETSC_COMM_WORLD, sizeof(Parameter), &user.bag));
499:   PetscCall(SetupParameters(&user));
500:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
501:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
502:   PetscCall(SNESSetDM(snes, dm));
503:   PetscCall(DMSetApplicationContext(dm, &user));
504:   /* Setup problem */
505:   PetscCall(SetupDiscretization(dm, &user));
506:   PetscCall(DMPlexCreateClosureIndex(dm, NULL));

508:   PetscCall(DMCreateGlobalVector(dm, &u));
509:   PetscCall(PetscObjectSetName((PetscObject)u, "Solution"));
510:   PetscCall(VecDuplicate(u, &r));

512:   PetscCall(DMSetNullSpaceConstructor(dm, 1, CreatePressureNullSpace));
513:   PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));

515:   PetscCall(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:     PetscCall(DMGetDS(dm, &ds));
522:     PetscCall(PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]));
523:     PetscCall(PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]));
524:     PetscCall(PetscDSGetExactSolution(ds, 2, &exactFuncs[2], &ctxs[2]));
525:     PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, u));
526:     PetscCall(PetscObjectSetName((PetscObject)u, "Exact Solution"));
527:     PetscCall(VecViewFromOptions(u, NULL, "-exact_vec_view"));
528:   }
529:   PetscCall(DMSNESCheckFromOptions(snes, u));
530:   PetscCall(VecSet(u, 0.0));
531:   PetscCall(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:     PetscCall(DMGetDS(dm, &ds));
540:     PetscCall(PetscDSGetExactSolution(ds, 0, &exactFuncs[0], &ctxs[0]));
541:     PetscCall(PetscDSGetExactSolution(ds, 1, &exactFuncs[1], &ctxs[1]));
542:     PetscCall(PetscDSGetExactSolution(ds, 2, &exactFuncs[2], &ctxs[2]));
543:     PetscCall(DMGetGlobalVector(dm, &r));
544:     PetscCall(DMProjectFunction(dm, 0.0, exactFuncs, ctxs, INSERT_ALL_VALUES, r));
545:     PetscCall(VecAXPY(r, -1.0, u));
546:     PetscCall(PetscObjectSetName((PetscObject)r, "Solution Error"));
547:     PetscCall(VecViewFromOptions(r, NULL, "-error_vec_view"));
548:     PetscCall(DMRestoreGlobalVector(dm, &r));
549:   }
550:   PetscCall(PetscObjectSetName((PetscObject)u, "Numerical Solution"));
551:   PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view"));

553:   PetscCall(VecDestroy(&u));
554:   PetscCall(VecDestroy(&r));
555:   PetscCall(DMDestroy(&dm));
556:   PetscCall(SNESDestroy(&snes));
557:   PetscCall(PetscBagDestroy(&user.bag));
558:   PetscCall(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*/