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