1 static char help[] = "Time dependent Navier-Stokes problem in 2d and 3d with finite elements.\n\ 2 We solve the Navier-Stokes in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\ 4 This example supports discretized auxiliary fields (Re) as well as\n\ 5 multilevel nonlinear solvers.\n\ 6 Contributed by: Julian Andrej <juan@tf.uni-kiel.de>\n\n\n"; 7 8 #include <petscdmplex.h> 9 #include <petscsnes.h> 10 #include <petscts.h> 11 #include <petscds.h> 12 13 /* 14 Navier-Stokes equation: 15 16 du/dt + u . grad u - \Delta u - grad p = f 17 div u = 0 18 */ 19 20 typedef struct { 21 PetscInt mms; 22 } AppCtx; 23 24 #define REYN 400.0 25 26 /* MMS1 27 28 u = t + x^2 + y^2; 29 v = t + 2*x^2 - 2*x*y; 30 p = x + y - 1; 31 32 f_x = -2*t*(x + y) + 2*x*y^2 - 4*x^2*y - 2*x^3 + 4.0/Re - 1.0 33 f_y = -2*t*x + 2*y^3 - 4*x*y^2 - 2*x^2*y + 4.0/Re - 1.0 34 35 so that 36 37 u_t + u \cdot \nabla u - 1/Re \Delta u + \nabla p + f = <1, 1> + <t (2x + 2y) + 2x^3 + 4x^2y - 2xy^2, t 2x + 2x^2y + 4xy^2 - 2y^3> - 1/Re <4, 4> + <1, 1> 38 + <-t (2x + 2y) + 2xy^2 - 4x^2y - 2x^3 + 4/Re - 1, -2xt + 2y^3 - 4xy^2 - 2x^2y + 4/Re - 1> = 0 39 \nabla \cdot u = 2x - 2x = 0 40 41 where 42 43 <u, v> . <<u_x, v_x>, <u_y, v_y>> = <u u_x + v u_y, u v_x + v v_y> 44 */ 45 PetscErrorCode mms1_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) { 46 u[0] = time + x[0] * x[0] + x[1] * x[1]; 47 u[1] = time + 2.0 * x[0] * x[0] - 2.0 * x[0] * x[1]; 48 return 0; 49 } 50 51 PetscErrorCode mms1_p_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx) { 52 *p = x[0] + x[1] - 1.0; 53 return 0; 54 } 55 56 /* MMS 2*/ 57 58 static PetscErrorCode mms2_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) { 59 u[0] = PetscSinReal(time + x[0]) * PetscSinReal(time + x[1]); 60 u[1] = PetscCosReal(time + x[0]) * PetscCosReal(time + x[1]); 61 return 0; 62 } 63 64 static PetscErrorCode mms2_p_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx) { 65 *p = PetscSinReal(time + x[0] - x[1]); 66 return 0; 67 } 68 69 static void f0_mms1_u(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[]) { 70 const PetscReal Re = REYN; 71 const PetscInt Ncomp = dim; 72 PetscInt c, d; 73 74 for (c = 0; c < Ncomp; ++c) { 75 for (d = 0; d < dim; ++d) f0[c] += u[d] * u_x[c * dim + d]; 76 } 77 f0[0] += u_t[0]; 78 f0[1] += u_t[1]; 79 80 f0[0] += -2.0 * t * (x[0] + x[1]) + 2.0 * x[0] * x[1] * x[1] - 4.0 * x[0] * x[0] * x[1] - 2.0 * x[0] * x[0] * x[0] + 4.0 / Re - 1.0; 81 f0[1] += -2.0 * t * x[0] + 2.0 * x[1] * x[1] * x[1] - 4.0 * x[0] * x[1] * x[1] - 2.0 * x[0] * x[0] * x[1] + 4.0 / Re - 1.0; 82 } 83 84 static void f0_mms2_u(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[]) { 85 const PetscReal Re = REYN; 86 const PetscInt Ncomp = dim; 87 PetscInt c, d; 88 89 for (c = 0; c < Ncomp; ++c) { 90 for (d = 0; d < dim; ++d) f0[c] += u[d] * u_x[c * dim + d]; 91 } 92 f0[0] += u_t[0]; 93 f0[1] += u_t[1]; 94 95 f0[0] -= (Re * ((1.0L / 2.0L) * PetscSinReal(2 * t + 2 * x[0]) + PetscSinReal(2 * t + x[0] + x[1]) + PetscCosReal(t + x[0] - x[1])) + 2.0 * PetscSinReal(t + x[0]) * PetscSinReal(t + x[1])) / Re; 96 f0[1] -= (-Re * ((1.0L / 2.0L) * PetscSinReal(2 * t + 2 * x[1]) + PetscSinReal(2 * t + x[0] + x[1]) + PetscCosReal(t + x[0] - x[1])) + 2.0 * PetscCosReal(t + x[0]) * PetscCosReal(t + x[1])) / Re; 97 } 98 99 static void f1_u(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[]) { 100 const PetscReal Re = REYN; 101 const PetscInt Ncomp = dim; 102 PetscInt comp, d; 103 104 for (comp = 0; comp < Ncomp; ++comp) { 105 for (d = 0; d < dim; ++d) f1[comp * dim + d] = 1.0 / Re * u_x[comp * dim + d]; 106 f1[comp * dim + comp] -= u[Ncomp]; 107 } 108 } 109 110 static void f0_p(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[]) { 111 PetscInt d; 112 for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d * dim + d]; 113 } 114 115 static void f1_p(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[]) { 116 PetscInt d; 117 for (d = 0; d < dim; ++d) f1[d] = 0.0; 118 } 119 120 /* 121 (psi_i, u_j grad_j u_i) ==> (\psi_i, \phi_j grad_j u_i) 122 */ 123 static void g0_uu(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[]) { 124 PetscInt NcI = dim, NcJ = dim; 125 PetscInt fc, gc; 126 PetscInt d; 127 128 for (d = 0; d < dim; ++d) g0[d * dim + d] = u_tShift; 129 130 for (fc = 0; fc < NcI; ++fc) { 131 for (gc = 0; gc < NcJ; ++gc) g0[fc * NcJ + gc] += u_x[fc * NcJ + gc]; 132 } 133 } 134 135 /* 136 (psi_i, u_j grad_j u_i) ==> (\psi_i, \u_j grad_j \phi_i) 137 */ 138 static void g1_uu(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[]) { 139 PetscInt NcI = dim; 140 PetscInt NcJ = dim; 141 PetscInt fc, gc, dg; 142 for (fc = 0; fc < NcI; ++fc) { 143 for (gc = 0; gc < NcJ; ++gc) { 144 for (dg = 0; dg < dim; ++dg) { 145 /* kronecker delta */ 146 if (fc == gc) g1[(fc * NcJ + gc) * dim + dg] += u[dg]; 147 } 148 } 149 } 150 } 151 152 /* < q, \nabla\cdot u > 153 NcompI = 1, NcompJ = dim */ 154 static void g1_pu(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[]) { 155 PetscInt d; 156 for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */ 157 } 158 159 /* -< \nabla\cdot v, p > 160 NcompI = dim, NcompJ = 1 */ 161 static void g2_up(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[]) { 162 PetscInt d; 163 for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */ 164 } 165 166 /* < \nabla v, \nabla u + {\nabla u}^T > 167 This just gives \nabla u, give the perdiagonal for the transpose */ 168 static void g3_uu(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[]) { 169 const PetscReal Re = REYN; 170 const PetscInt Ncomp = dim; 171 PetscInt compI, d; 172 173 for (compI = 0; compI < Ncomp; ++compI) { 174 for (d = 0; d < dim; ++d) g3[((compI * Ncomp + compI) * dim + d) * dim + d] = 1.0 / Re; 175 } 176 } 177 178 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) { 179 PetscFunctionBeginUser; 180 options->mms = 1; 181 182 PetscOptionsBegin(comm, "", "Navier-Stokes Equation Options", "DMPLEX"); 183 PetscCall(PetscOptionsInt("-mms", "The manufactured solution to use", "ex46.c", options->mms, &options->mms, NULL)); 184 PetscOptionsEnd(); 185 PetscFunctionReturn(0); 186 } 187 188 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx) { 189 PetscFunctionBeginUser; 190 PetscCall(DMCreate(comm, dm)); 191 PetscCall(DMSetType(*dm, DMPLEX)); 192 PetscCall(DMSetFromOptions(*dm)); 193 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 194 PetscFunctionReturn(0); 195 } 196 197 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx) { 198 PetscDS ds; 199 DMLabel label; 200 const PetscInt id = 1; 201 PetscInt dim; 202 203 PetscFunctionBeginUser; 204 PetscCall(DMGetDimension(dm, &dim)); 205 PetscCall(DMGetDS(dm, &ds)); 206 PetscCall(DMGetLabel(dm, "marker", &label)); 207 switch (dim) { 208 case 2: 209 switch (ctx->mms) { 210 case 1: 211 PetscCall(PetscDSSetResidual(ds, 0, f0_mms1_u, f1_u)); 212 PetscCall(PetscDSSetResidual(ds, 1, f0_p, f1_p)); 213 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, g1_uu, NULL, g3_uu)); 214 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL)); 215 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL)); 216 PetscCall(PetscDSSetExactSolution(ds, 0, mms1_u_2d, ctx)); 217 PetscCall(PetscDSSetExactSolution(ds, 1, mms1_p_2d, ctx)); 218 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))mms1_u_2d, NULL, ctx, NULL)); 219 break; 220 case 2: 221 PetscCall(PetscDSSetResidual(ds, 0, f0_mms2_u, f1_u)); 222 PetscCall(PetscDSSetResidual(ds, 1, f0_p, f1_p)); 223 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_uu, g1_uu, NULL, g3_uu)); 224 PetscCall(PetscDSSetJacobian(ds, 0, 1, NULL, NULL, g2_up, NULL)); 225 PetscCall(PetscDSSetJacobian(ds, 1, 0, NULL, g1_pu, NULL, NULL)); 226 PetscCall(PetscDSSetExactSolution(ds, 0, mms2_u_2d, ctx)); 227 PetscCall(PetscDSSetExactSolution(ds, 1, mms2_p_2d, ctx)); 228 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))mms2_u_2d, NULL, ctx, NULL)); 229 break; 230 default: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid MMS %" PetscInt_FMT, ctx->mms); 231 } 232 break; 233 default: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %" PetscInt_FMT, dim); 234 } 235 PetscFunctionReturn(0); 236 } 237 238 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx) { 239 MPI_Comm comm; 240 DM cdm = dm; 241 PetscFE fe[2]; 242 PetscInt dim; 243 PetscBool simplex; 244 245 PetscFunctionBeginUser; 246 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 247 PetscCall(DMGetDimension(dm, &dim)); 248 PetscCall(DMPlexIsSimplex(dm, &simplex)); 249 PetscCall(PetscFECreateDefault(comm, dim, dim, simplex, "vel_", PETSC_DEFAULT, &fe[0])); 250 PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity")); 251 PetscCall(PetscFECreateDefault(comm, dim, 1, simplex, "pres_", PETSC_DEFAULT, &fe[1])); 252 PetscCall(PetscFECopyQuadrature(fe[0], fe[1])); 253 PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure")); 254 /* Set discretization and boundary conditions for each mesh */ 255 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe[0])); 256 PetscCall(DMSetField(dm, 1, NULL, (PetscObject)fe[1])); 257 PetscCall(DMCreateDS(dm)); 258 PetscCall(SetupProblem(dm, ctx)); 259 while (cdm) { 260 PetscObject pressure; 261 MatNullSpace nsp; 262 263 PetscCall(DMGetField(cdm, 1, NULL, &pressure)); 264 PetscCall(MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nsp)); 265 PetscCall(PetscObjectCompose(pressure, "nullspace", (PetscObject)nsp)); 266 PetscCall(MatNullSpaceDestroy(&nsp)); 267 268 PetscCall(DMCopyDisc(dm, cdm)); 269 PetscCall(DMGetCoarseDM(cdm, &cdm)); 270 } 271 PetscCall(PetscFEDestroy(&fe[0])); 272 PetscCall(PetscFEDestroy(&fe[1])); 273 PetscFunctionReturn(0); 274 } 275 276 static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx) { 277 PetscSimplePointFunc funcs[2]; 278 void *ctxs[2]; 279 DM dm; 280 PetscDS ds; 281 PetscReal ferrors[2]; 282 283 PetscFunctionBeginUser; 284 PetscCall(TSGetDM(ts, &dm)); 285 PetscCall(DMGetDS(dm, &ds)); 286 PetscCall(PetscDSGetExactSolution(ds, 0, &funcs[0], &ctxs[0])); 287 PetscCall(PetscDSGetExactSolution(ds, 1, &funcs[1], &ctxs[1])); 288 PetscCall(DMComputeL2FieldDiff(dm, crtime, funcs, ctxs, u, ferrors)); 289 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Timestep: %04d time = %-8.4g \t L_2 Error: [%2.3g, %2.3g]\n", (int)step, (double)crtime, (double)ferrors[0], (double)ferrors[1])); 290 PetscFunctionReturn(0); 291 } 292 293 int main(int argc, char **argv) { 294 AppCtx ctx; 295 DM dm; 296 TS ts; 297 Vec u, r; 298 299 PetscFunctionBeginUser; 300 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 301 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx)); 302 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm, &ctx)); 303 PetscCall(DMSetApplicationContext(dm, &ctx)); 304 PetscCall(SetupDiscretization(dm, &ctx)); 305 PetscCall(DMPlexCreateClosureIndex(dm, NULL)); 306 307 PetscCall(DMCreateGlobalVector(dm, &u)); 308 PetscCall(VecDuplicate(u, &r)); 309 310 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 311 PetscCall(TSMonitorSet(ts, MonitorError, &ctx, NULL)); 312 PetscCall(TSSetDM(ts, dm)); 313 PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx)); 314 PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx)); 315 PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx)); 316 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 317 PetscCall(TSSetFromOptions(ts)); 318 PetscCall(DMTSCheckFromOptions(ts, u)); 319 320 { 321 PetscSimplePointFunc funcs[2]; 322 void *ctxs[2]; 323 PetscDS ds; 324 325 PetscCall(DMGetDS(dm, &ds)); 326 PetscCall(PetscDSGetExactSolution(ds, 0, &funcs[0], &ctxs[0])); 327 PetscCall(PetscDSGetExactSolution(ds, 1, &funcs[1], &ctxs[1])); 328 PetscCall(DMProjectFunction(dm, 0.0, funcs, ctxs, INSERT_ALL_VALUES, u)); 329 } 330 PetscCall(TSSolve(ts, u)); 331 PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view")); 332 333 PetscCall(VecDestroy(&u)); 334 PetscCall(VecDestroy(&r)); 335 PetscCall(TSDestroy(&ts)); 336 PetscCall(DMDestroy(&dm)); 337 PetscCall(PetscFinalize()); 338 return 0; 339 } 340 341 /*TEST 342 343 # Full solves 344 test: 345 suffix: 2d_p2p1_r1 346 requires: !single triangle 347 filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g" 348 args: -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \ 349 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -ts_monitor -dmts_check \ 350 -snes_monitor_short -snes_converged_reason \ 351 -ksp_monitor_short -ksp_converged_reason \ 352 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full \ 353 -fieldsplit_velocity_pc_type lu \ 354 -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi 355 356 test: 357 suffix: 2d_q2q1_r1 358 requires: !single 359 filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g" -e "s~ 0\]~ 0.0\]~g" 360 args: -dm_plex_simplex 0 -dm_refine 1 -vel_petscspace_degree 2 -pres_petscspace_degree 1 \ 361 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -ts_monitor -dmts_check \ 362 -snes_monitor_short -snes_converged_reason \ 363 -ksp_monitor_short -ksp_converged_reason \ 364 -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full \ 365 -fieldsplit_velocity_pc_type lu \ 366 -fieldsplit_pressure_ksp_rtol 1.0e-10 -fieldsplit_pressure_pc_type jacobi 367 368 TEST*/ 369