1 static char help[] = "Poisson Problem in 2d and 3d with finite elements.\n\ 2 We solve the Poisson problem in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\ 4 This example supports automatic convergence estimation\n\ 5 and eventually adaptivity.\n\n\n"; 6 7 #include <petscdmplex.h> 8 #include <petscsnes.h> 9 #include <petscds.h> 10 #include <petscconvest.h> 11 12 typedef struct { 13 /* Domain and mesh definition */ 14 PetscBool spectral; /* Look at the spectrum along planes in the solution */ 15 PetscBool shear; /* Shear the domain */ 16 PetscBool adjoint; /* Solve the adjoint problem */ 17 PetscBool homogeneous; /* Use homogeneous boudnary conditions */ 18 PetscBool viewError; /* Output the solution error */ 19 } AppCtx; 20 21 static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 22 { 23 *u = 0.0; 24 return 0; 25 } 26 27 static PetscErrorCode trig_inhomogeneous_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 28 { 29 PetscInt d; 30 *u = 0.0; 31 for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0 * PETSC_PI * x[d]); 32 return 0; 33 } 34 35 static PetscErrorCode trig_homogeneous_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 36 { 37 PetscInt d; 38 *u = 1.0; 39 for (d = 0; d < dim; ++d) *u *= PetscSinReal(2.0 * PETSC_PI * x[d]); 40 return 0; 41 } 42 43 /* Compute integral of (residual of solution)*(adjoint solution - projection of adjoint solution) */ 44 static void obj_error_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 obj[]) 45 { 46 obj[0] = a[aOff[0]] * (u[0] - a[aOff[1]]); 47 } 48 49 static void f0_trig_inhomogeneous_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[]) 50 { 51 PetscInt d; 52 for (d = 0; d < dim; ++d) f0[0] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]); 53 } 54 55 static void f0_trig_homogeneous_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[]) 56 { 57 PetscInt d; 58 for (d = 0; d < dim; ++d) { 59 PetscScalar v = 1.; 60 for (PetscInt e = 0; e < dim; e++) { 61 if (e == d) { 62 v *= -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]); 63 } else { 64 v *= PetscSinReal(2.0 * PETSC_PI * x[d]); 65 } 66 } 67 f0[0] += v; 68 } 69 } 70 71 static void f0_unity_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[]) 72 { 73 f0[0] = 1.0; 74 } 75 76 static void f0_identityaux_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[]) 77 { 78 f0[0] = a[0]; 79 } 80 81 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[]) 82 { 83 PetscInt d; 84 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 85 } 86 87 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[]) 88 { 89 PetscInt d; 90 for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0; 91 } 92 93 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 94 { 95 PetscFunctionBeginUser; 96 options->shear = PETSC_FALSE; 97 options->spectral = PETSC_FALSE; 98 options->adjoint = PETSC_FALSE; 99 options->homogeneous = PETSC_FALSE; 100 options->viewError = PETSC_FALSE; 101 102 PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX"); 103 PetscCall(PetscOptionsBool("-shear", "Shear the domain", "ex13.c", options->shear, &options->shear, NULL)); 104 PetscCall(PetscOptionsBool("-spectral", "Look at the spectrum along planes of the solution", "ex13.c", options->spectral, &options->spectral, NULL)); 105 PetscCall(PetscOptionsBool("-adjoint", "Solve the adjoint problem", "ex13.c", options->adjoint, &options->adjoint, NULL)); 106 PetscCall(PetscOptionsBool("-homogeneous", "Use homogeneous boundary conditions", "ex13.c", options->homogeneous, &options->homogeneous, NULL)); 107 PetscCall(PetscOptionsBool("-error_view", "Output the solution error", "ex13.c", options->viewError, &options->viewError, NULL)); 108 PetscOptionsEnd(); 109 PetscFunctionReturn(0); 110 } 111 112 static PetscErrorCode CreateSpectralPlanes(DM dm, PetscInt numPlanes, const PetscInt planeDir[], const PetscReal planeCoord[], AppCtx *user) 113 { 114 PetscSection coordSection; 115 Vec coordinates; 116 const PetscScalar *coords; 117 PetscInt dim, p, vStart, vEnd, v; 118 119 PetscFunctionBeginUser; 120 PetscCall(DMGetCoordinateDim(dm, &dim)); 121 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 122 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 123 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 124 PetscCall(VecGetArrayRead(coordinates, &coords)); 125 for (p = 0; p < numPlanes; ++p) { 126 DMLabel label; 127 char name[PETSC_MAX_PATH_LEN]; 128 129 PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "spectral_plane_%" PetscInt_FMT, p)); 130 PetscCall(DMCreateLabel(dm, name)); 131 PetscCall(DMGetLabel(dm, name, &label)); 132 PetscCall(DMLabelAddStratum(label, 1)); 133 for (v = vStart; v < vEnd; ++v) { 134 PetscInt off; 135 136 PetscCall(PetscSectionGetOffset(coordSection, v, &off)); 137 if (PetscAbsReal(planeCoord[p] - PetscRealPart(coords[off + planeDir[p]])) < PETSC_SMALL) PetscCall(DMLabelSetValue(label, v, 1)); 138 } 139 } 140 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 141 PetscFunctionReturn(0); 142 } 143 144 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 145 { 146 PetscFunctionBeginUser; 147 PetscCall(DMCreate(comm, dm)); 148 PetscCall(DMSetType(*dm, DMPLEX)); 149 PetscCall(DMSetFromOptions(*dm)); 150 if (user->shear) PetscCall(DMPlexShearGeometry(*dm, DM_X, NULL)); 151 PetscCall(DMSetApplicationContext(*dm, user)); 152 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 153 if (user->spectral) { 154 PetscInt planeDir[2] = {0, 1}; 155 PetscReal planeCoord[2] = {0., 1.}; 156 157 PetscCall(CreateSpectralPlanes(*dm, 2, planeDir, planeCoord, user)); 158 } 159 PetscFunctionReturn(0); 160 } 161 162 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 163 { 164 PetscDS ds; 165 DMLabel label; 166 const PetscInt id = 1; 167 PetscPointFunc f0 = user->homogeneous ? f0_trig_homogeneous_u : f0_trig_inhomogeneous_u; 168 PetscErrorCode (*ex)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = user->homogeneous ? trig_homogeneous_u : trig_inhomogeneous_u; 169 170 PetscFunctionBeginUser; 171 PetscCall(DMGetDS(dm, &ds)); 172 PetscCall(PetscDSSetResidual(ds, 0, f0, f1_u)); 173 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 174 PetscCall(PetscDSSetExactSolution(ds, 0, ex, user)); 175 PetscCall(DMGetLabel(dm, "marker", &label)); 176 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))ex, NULL, user, NULL)); 177 PetscFunctionReturn(0); 178 } 179 180 static PetscErrorCode SetupAdjointProblem(DM dm, AppCtx *user) 181 { 182 PetscDS ds; 183 DMLabel label; 184 const PetscInt id = 1; 185 186 PetscFunctionBeginUser; 187 PetscCall(DMGetDS(dm, &ds)); 188 PetscCall(PetscDSSetResidual(ds, 0, f0_unity_u, f1_u)); 189 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu)); 190 PetscCall(PetscDSSetObjective(ds, 0, obj_error_u)); 191 PetscCall(DMGetLabel(dm, "marker", &label)); 192 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))zero, NULL, user, NULL)); 193 PetscFunctionReturn(0); 194 } 195 196 static PetscErrorCode SetupErrorProblem(DM dm, AppCtx *user) 197 { 198 PetscDS prob; 199 200 PetscFunctionBeginUser; 201 PetscCall(DMGetDS(dm, &prob)); 202 PetscFunctionReturn(0); 203 } 204 205 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 206 { 207 DM cdm = dm; 208 PetscFE fe; 209 DMPolytopeType ct; 210 PetscBool simplex; 211 PetscInt dim, cStart; 212 char prefix[PETSC_MAX_PATH_LEN]; 213 214 PetscFunctionBeginUser; 215 PetscCall(DMGetDimension(dm, &dim)); 216 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 217 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 218 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE; 219 /* Create finite element */ 220 PetscCall(PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name)); 221 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe)); 222 PetscCall(PetscObjectSetName((PetscObject)fe, name)); 223 /* Set discretization and boundary conditions for each mesh */ 224 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 225 PetscCall(DMCreateDS(dm)); 226 PetscCall((*setup)(dm, user)); 227 while (cdm) { 228 PetscCall(DMCopyDisc(dm, cdm)); 229 /* TODO: Check whether the boundary of coarse meshes is marked */ 230 PetscCall(DMGetCoarseDM(cdm, &cdm)); 231 } 232 PetscCall(PetscFEDestroy(&fe)); 233 PetscFunctionReturn(0); 234 } 235 236 static PetscErrorCode ComputeSpectral(DM dm, Vec u, PetscInt numPlanes, const PetscInt planeDir[], const PetscReal planeCoord[], AppCtx *user) 237 { 238 MPI_Comm comm; 239 PetscSection coordSection, section; 240 Vec coordinates, uloc; 241 const PetscScalar *coords, *array; 242 PetscInt p; 243 PetscMPIInt size, rank; 244 245 PetscFunctionBeginUser; 246 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 247 PetscCallMPI(MPI_Comm_size(comm, &size)); 248 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 249 PetscCall(DMGetLocalVector(dm, &uloc)); 250 PetscCall(DMGlobalToLocalBegin(dm, u, INSERT_VALUES, uloc)); 251 PetscCall(DMGlobalToLocalEnd(dm, u, INSERT_VALUES, uloc)); 252 PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, uloc, 0.0, NULL, NULL, NULL)); 253 PetscCall(VecViewFromOptions(uloc, NULL, "-sol_view")); 254 PetscCall(DMGetLocalSection(dm, §ion)); 255 PetscCall(VecGetArrayRead(uloc, &array)); 256 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 257 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 258 PetscCall(VecGetArrayRead(coordinates, &coords)); 259 for (p = 0; p < numPlanes; ++p) { 260 DMLabel label; 261 char name[PETSC_MAX_PATH_LEN]; 262 Mat F; 263 Vec x, y; 264 IS stratum; 265 PetscReal *ray, *gray; 266 PetscScalar *rvals, *svals, *gsvals; 267 PetscInt *perm, *nperm; 268 PetscInt n, N, i, j, off, offu; 269 const PetscInt *points; 270 271 PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "spectral_plane_%" PetscInt_FMT, p)); 272 PetscCall(DMGetLabel(dm, name, &label)); 273 PetscCall(DMLabelGetStratumIS(label, 1, &stratum)); 274 PetscCall(ISGetLocalSize(stratum, &n)); 275 PetscCall(ISGetIndices(stratum, &points)); 276 PetscCall(PetscMalloc2(n, &ray, n, &svals)); 277 for (i = 0; i < n; ++i) { 278 PetscCall(PetscSectionGetOffset(coordSection, points[i], &off)); 279 PetscCall(PetscSectionGetOffset(section, points[i], &offu)); 280 ray[i] = PetscRealPart(coords[off + ((planeDir[p] + 1) % 2)]); 281 svals[i] = array[offu]; 282 } 283 /* Gather the ray data to proc 0 */ 284 if (size > 1) { 285 PetscMPIInt *cnt, *displs, p; 286 287 PetscCall(PetscCalloc2(size, &cnt, size, &displs)); 288 PetscCallMPI(MPI_Gather(&n, 1, MPIU_INT, cnt, 1, MPIU_INT, 0, comm)); 289 for (p = 1; p < size; ++p) displs[p] = displs[p - 1] + cnt[p - 1]; 290 N = displs[size - 1] + cnt[size - 1]; 291 PetscCall(PetscMalloc2(N, &gray, N, &gsvals)); 292 PetscCallMPI(MPI_Gatherv(ray, n, MPIU_REAL, gray, cnt, displs, MPIU_REAL, 0, comm)); 293 PetscCallMPI(MPI_Gatherv(svals, n, MPIU_SCALAR, gsvals, cnt, displs, MPIU_SCALAR, 0, comm)); 294 PetscCall(PetscFree2(cnt, displs)); 295 } else { 296 N = n; 297 gray = ray; 298 gsvals = svals; 299 } 300 if (rank == 0) { 301 /* Sort point along ray */ 302 PetscCall(PetscMalloc2(N, &perm, N, &nperm)); 303 for (i = 0; i < N; ++i) perm[i] = i; 304 PetscCall(PetscSortRealWithPermutation(N, gray, perm)); 305 /* Count duplicates and squish mapping */ 306 nperm[0] = perm[0]; 307 for (i = 1, j = 1; i < N; ++i) { 308 if (PetscAbsReal(gray[perm[i]] - gray[perm[i - 1]]) > PETSC_SMALL) nperm[j++] = perm[i]; 309 } 310 /* Create FFT structs */ 311 PetscCall(MatCreateFFT(PETSC_COMM_SELF, 1, &j, MATFFTW, &F)); 312 PetscCall(MatCreateVecs(F, &x, &y)); 313 PetscCall(PetscObjectSetName((PetscObject)y, name)); 314 PetscCall(VecGetArray(x, &rvals)); 315 for (i = 0, j = 0; i < N; ++i) { 316 if (i > 0 && PetscAbsReal(gray[perm[i]] - gray[perm[i - 1]]) < PETSC_SMALL) continue; 317 rvals[j] = gsvals[nperm[j]]; 318 ++j; 319 } 320 PetscCall(PetscFree2(perm, nperm)); 321 if (size > 1) PetscCall(PetscFree2(gray, gsvals)); 322 PetscCall(VecRestoreArray(x, &rvals)); 323 /* Do FFT along the ray */ 324 PetscCall(MatMult(F, x, y)); 325 /* Chop FFT */ 326 PetscCall(VecChop(y, PETSC_SMALL)); 327 PetscCall(VecViewFromOptions(x, NULL, "-real_view")); 328 PetscCall(VecViewFromOptions(y, NULL, "-fft_view")); 329 PetscCall(VecDestroy(&x)); 330 PetscCall(VecDestroy(&y)); 331 PetscCall(MatDestroy(&F)); 332 } 333 PetscCall(ISRestoreIndices(stratum, &points)); 334 PetscCall(ISDestroy(&stratum)); 335 PetscCall(PetscFree2(ray, svals)); 336 } 337 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 338 PetscCall(VecRestoreArrayRead(uloc, &array)); 339 PetscCall(DMRestoreLocalVector(dm, &uloc)); 340 PetscFunctionReturn(0); 341 } 342 343 int main(int argc, char **argv) 344 { 345 DM dm; /* Problem specification */ 346 SNES snes; /* Nonlinear solver */ 347 Vec u; /* Solutions */ 348 AppCtx user; /* User-defined work context */ 349 350 PetscFunctionBeginUser; 351 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 352 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 353 /* Primal system */ 354 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 355 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 356 PetscCall(SNESSetDM(snes, dm)); 357 PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user)); 358 PetscCall(DMCreateGlobalVector(dm, &u)); 359 PetscCall(VecSet(u, 0.0)); 360 PetscCall(PetscObjectSetName((PetscObject)u, "potential")); 361 PetscCall(DMPlexSetSNESLocalFEM(dm, &user, &user, &user)); 362 PetscCall(SNESSetFromOptions(snes)); 363 PetscCall(SNESSolve(snes, NULL, u)); 364 PetscCall(SNESGetSolution(snes, &u)); 365 PetscCall(VecViewFromOptions(u, NULL, "-potential_view")); 366 if (user.viewError) { 367 PetscErrorCode (*sol)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *); 368 void *ctx; 369 PetscDS ds; 370 PetscReal error; 371 PetscInt N; 372 373 PetscCall(DMGetDS(dm, &ds)); 374 PetscCall(PetscDSGetExactSolution(ds, 0, &sol, &ctx)); 375 PetscCall(VecGetSize(u, &N)); 376 PetscCall(DMComputeL2Diff(dm, 0.0, &sol, &ctx, u, &error)); 377 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "N: %" PetscInt_FMT " L2 error: %g\n", N, (double)error)); 378 } 379 if (user.spectral) { 380 PetscInt planeDir[2] = {0, 1}; 381 PetscReal planeCoord[2] = {0., 1.}; 382 383 PetscCall(ComputeSpectral(dm, u, 2, planeDir, planeCoord, &user)); 384 } 385 /* Adjoint system */ 386 if (user.adjoint) { 387 DM dmAdj; 388 SNES snesAdj; 389 Vec uAdj; 390 391 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snesAdj)); 392 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)snesAdj, "adjoint_")); 393 PetscCall(DMClone(dm, &dmAdj)); 394 PetscCall(SNESSetDM(snesAdj, dmAdj)); 395 PetscCall(SetupDiscretization(dmAdj, "adjoint", SetupAdjointProblem, &user)); 396 PetscCall(DMCreateGlobalVector(dmAdj, &uAdj)); 397 PetscCall(VecSet(uAdj, 0.0)); 398 PetscCall(PetscObjectSetName((PetscObject)uAdj, "adjoint")); 399 PetscCall(DMPlexSetSNESLocalFEM(dmAdj, &user, &user, &user)); 400 PetscCall(SNESSetFromOptions(snesAdj)); 401 PetscCall(SNESSolve(snesAdj, NULL, uAdj)); 402 PetscCall(SNESGetSolution(snesAdj, &uAdj)); 403 PetscCall(VecViewFromOptions(uAdj, NULL, "-adjoint_view")); 404 /* Error representation */ 405 { 406 DM dmErr, dmErrAux, dms[2]; 407 Vec errorEst, errorL2, uErr, uErrLoc, uAdjLoc, uAdjProj; 408 IS *subis; 409 PetscReal errorEstTot, errorL2Norm, errorL2Tot; 410 PetscInt N, i; 411 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {user.homogeneous ? trig_homogeneous_u : trig_inhomogeneous_u}; 412 void (*identity[1])(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]) = {f0_identityaux_u}; 413 void *ctxs[1] = {0}; 414 415 ctxs[0] = &user; 416 PetscCall(DMClone(dm, &dmErr)); 417 PetscCall(SetupDiscretization(dmErr, "error", SetupErrorProblem, &user)); 418 PetscCall(DMGetGlobalVector(dmErr, &errorEst)); 419 PetscCall(DMGetGlobalVector(dmErr, &errorL2)); 420 /* Compute auxiliary data (solution and projection of adjoint solution) */ 421 PetscCall(DMGetLocalVector(dmAdj, &uAdjLoc)); 422 PetscCall(DMGlobalToLocalBegin(dmAdj, uAdj, INSERT_VALUES, uAdjLoc)); 423 PetscCall(DMGlobalToLocalEnd(dmAdj, uAdj, INSERT_VALUES, uAdjLoc)); 424 PetscCall(DMGetGlobalVector(dm, &uAdjProj)); 425 PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, uAdjLoc)); 426 PetscCall(DMProjectField(dm, 0.0, u, identity, INSERT_VALUES, uAdjProj)); 427 PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL)); 428 PetscCall(DMRestoreLocalVector(dmAdj, &uAdjLoc)); 429 /* Attach auxiliary data */ 430 dms[0] = dm; 431 dms[1] = dm; 432 PetscCall(DMCreateSuperDM(dms, 2, &subis, &dmErrAux)); 433 if (0) { 434 PetscSection sec; 435 436 PetscCall(DMGetLocalSection(dms[0], &sec)); 437 PetscCall(PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD)); 438 PetscCall(DMGetLocalSection(dms[1], &sec)); 439 PetscCall(PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD)); 440 PetscCall(DMGetLocalSection(dmErrAux, &sec)); 441 PetscCall(PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD)); 442 } 443 PetscCall(DMViewFromOptions(dmErrAux, NULL, "-dm_err_view")); 444 PetscCall(ISViewFromOptions(subis[0], NULL, "-super_is_view")); 445 PetscCall(ISViewFromOptions(subis[1], NULL, "-super_is_view")); 446 PetscCall(DMGetGlobalVector(dmErrAux, &uErr)); 447 PetscCall(VecViewFromOptions(u, NULL, "-map_vec_view")); 448 PetscCall(VecViewFromOptions(uAdjProj, NULL, "-map_vec_view")); 449 PetscCall(VecViewFromOptions(uErr, NULL, "-map_vec_view")); 450 PetscCall(VecISCopy(uErr, subis[0], SCATTER_FORWARD, u)); 451 PetscCall(VecISCopy(uErr, subis[1], SCATTER_FORWARD, uAdjProj)); 452 PetscCall(DMRestoreGlobalVector(dm, &uAdjProj)); 453 for (i = 0; i < 2; ++i) PetscCall(ISDestroy(&subis[i])); 454 PetscCall(PetscFree(subis)); 455 PetscCall(DMGetLocalVector(dmErrAux, &uErrLoc)); 456 PetscCall(DMGlobalToLocalBegin(dm, uErr, INSERT_VALUES, uErrLoc)); 457 PetscCall(DMGlobalToLocalEnd(dm, uErr, INSERT_VALUES, uErrLoc)); 458 PetscCall(DMRestoreGlobalVector(dmErrAux, &uErr)); 459 PetscCall(DMSetAuxiliaryVec(dmAdj, NULL, 0, 0, uErrLoc)); 460 /* Compute cellwise error estimate */ 461 PetscCall(VecSet(errorEst, 0.0)); 462 PetscCall(DMPlexComputeCellwiseIntegralFEM(dmAdj, uAdj, errorEst, &user)); 463 PetscCall(DMSetAuxiliaryVec(dmAdj, NULL, 0, 0, NULL)); 464 PetscCall(DMRestoreLocalVector(dmErrAux, &uErrLoc)); 465 PetscCall(DMDestroy(&dmErrAux)); 466 /* Plot cellwise error vector */ 467 PetscCall(VecViewFromOptions(errorEst, NULL, "-error_view")); 468 /* Compute ratio of estimate (sum over cells) with actual L_2 error */ 469 PetscCall(DMComputeL2Diff(dm, 0.0, funcs, ctxs, u, &errorL2Norm)); 470 PetscCall(DMPlexComputeL2DiffVec(dm, 0.0, funcs, ctxs, u, errorL2)); 471 PetscCall(VecViewFromOptions(errorL2, NULL, "-l2_error_view")); 472 PetscCall(VecNorm(errorL2, NORM_INFINITY, &errorL2Tot)); 473 PetscCall(VecNorm(errorEst, NORM_INFINITY, &errorEstTot)); 474 PetscCall(VecGetSize(errorEst, &N)); 475 PetscCall(VecPointwiseDivide(errorEst, errorEst, errorL2)); 476 PetscCall(PetscObjectSetName((PetscObject)errorEst, "Error ratio")); 477 PetscCall(VecViewFromOptions(errorEst, NULL, "-error_ratio_view")); 478 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "N: %" PetscInt_FMT " L2 error: %g Error Ratio: %g/%g = %g\n", N, (double)errorL2Norm, (double)errorEstTot, (double)PetscSqrtReal(errorL2Tot), (double)(errorEstTot / PetscSqrtReal(errorL2Tot)))); 479 PetscCall(DMRestoreGlobalVector(dmErr, &errorEst)); 480 PetscCall(DMRestoreGlobalVector(dmErr, &errorL2)); 481 PetscCall(DMDestroy(&dmErr)); 482 } 483 PetscCall(DMDestroy(&dmAdj)); 484 PetscCall(VecDestroy(&uAdj)); 485 PetscCall(SNESDestroy(&snesAdj)); 486 } 487 /* Cleanup */ 488 PetscCall(VecDestroy(&u)); 489 PetscCall(SNESDestroy(&snes)); 490 PetscCall(DMDestroy(&dm)); 491 PetscCall(PetscFinalize()); 492 return 0; 493 } 494 495 /*TEST 496 497 test: 498 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9 499 suffix: 2d_p1_conv 500 requires: triangle 501 args: -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 502 test: 503 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9 504 suffix: 2d_p2_conv 505 requires: triangle 506 args: -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 507 test: 508 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9 509 suffix: 2d_p3_conv 510 requires: triangle 511 args: -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2 512 test: 513 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9 514 suffix: 2d_q1_conv 515 args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 516 test: 517 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9 518 suffix: 2d_q2_conv 519 args: -dm_plex_simplex 0 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 520 test: 521 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9 522 suffix: 2d_q3_conv 523 args: -dm_plex_simplex 0 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2 524 test: 525 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9 526 suffix: 2d_q1_shear_conv 527 args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 528 test: 529 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9 530 suffix: 2d_q2_shear_conv 531 args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 532 test: 533 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9 534 suffix: 2d_q3_shear_conv 535 args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2 536 test: 537 # Using -convest_num_refine 3 we get L_2 convergence rate: 1.7 538 suffix: 3d_p1_conv 539 requires: ctetgen 540 args: -dm_plex_dim 3 -dm_refine 1 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 541 test: 542 # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 2.8 543 suffix: 3d_p2_conv 544 requires: ctetgen 545 args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 1 546 test: 547 # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 4.0 548 suffix: 3d_p3_conv 549 requires: ctetgen 550 args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 1 551 test: 552 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.8 553 suffix: 3d_q1_conv 554 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 1 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 555 test: 556 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.8 557 suffix: 3d_q2_conv 558 args: -dm_plex_dim 3 -dm_plex_simplex 0 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 1 559 test: 560 # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 3.8 561 suffix: 3d_q3_conv 562 args: -dm_plex_dim 3 -dm_plex_simplex 0 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 1 563 test: 564 suffix: 2d_p1_fas_full 565 requires: triangle 566 args: -potential_petscspace_degree 1 -dm_refine_hierarchy 5 \ 567 -snes_max_it 1 -snes_type fas -snes_fas_levels 5 -snes_fas_type full -snes_fas_full_total \ 568 -fas_coarse_snes_monitor -fas_coarse_snes_max_it 1 -fas_coarse_ksp_atol 1.e-13 \ 569 -fas_levels_snes_monitor -fas_levels_snes_max_it 1 -fas_levels_snes_type newtonls \ 570 -fas_levels_pc_type none -fas_levels_ksp_max_it 2 -fas_levels_ksp_converged_maxits -fas_levels_ksp_type chebyshev \ 571 -fas_levels_esteig_ksp_type cg -fas_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -fas_levels_esteig_ksp_max_it 10 572 test: 573 suffix: 2d_p1_fas_full_homogeneous 574 requires: triangle 575 args: -homogeneous -potential_petscspace_degree 1 -dm_refine_hierarchy 5 \ 576 -snes_max_it 1 -snes_type fas -snes_fas_levels 5 -snes_fas_type full \ 577 -fas_coarse_snes_monitor -fas_coarse_snes_max_it 1 -fas_coarse_ksp_atol 1.e-13 \ 578 -fas_levels_snes_monitor -fas_levels_snes_max_it 1 -fas_levels_snes_type newtonls \ 579 -fas_levels_pc_type none -fas_levels_ksp_max_it 2 -fas_levels_ksp_converged_maxits -fas_levels_ksp_type chebyshev \ 580 -fas_levels_esteig_ksp_type cg -fas_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -fas_levels_esteig_ksp_max_it 10 581 582 test: 583 suffix: 2d_p1_scalable 584 requires: triangle 585 args: -potential_petscspace_degree 1 -dm_refine 3 \ 586 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned \ 587 -pc_type gamg -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 \ 588 -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \ 589 -pc_gamg_coarse_eq_limit 1000 \ 590 -pc_gamg_threshold 0.05 \ 591 -pc_gamg_threshold_scale .0 \ 592 -mg_levels_ksp_type chebyshev \ 593 -mg_levels_ksp_max_it 1 \ 594 -mg_levels_pc_type jacobi \ 595 -matptap_via scalable 596 test: 597 suffix: 2d_p1_gmg_vcycle 598 requires: triangle 599 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 600 -ksp_rtol 5e-10 -pc_type mg \ 601 -mg_levels_ksp_max_it 1 \ 602 -mg_levels_esteig_ksp_type cg \ 603 -mg_levels_esteig_ksp_max_it 10 \ 604 -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \ 605 -mg_levels_pc_type jacobi 606 test: 607 suffix: 2d_p1_gmg_fcycle 608 requires: triangle 609 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 610 -ksp_rtol 5e-10 -pc_type mg -pc_mg_type full \ 611 -mg_levels_ksp_max_it 2 \ 612 -mg_levels_esteig_ksp_type cg \ 613 -mg_levels_esteig_ksp_max_it 10 \ 614 -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \ 615 -mg_levels_pc_type jacobi 616 test: 617 suffix: 2d_p1_gmg_vcycle_adapt 618 requires: triangle 619 args: -petscpartitioner_type simple -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 620 -ksp_rtol 5e-10 -pc_type mg -pc_mg_galerkin -pc_mg_adapt_interp_coarse_space harmonic -pc_mg_adapt_interp_n 8 \ 621 -mg_levels_ksp_max_it 1 \ 622 -mg_levels_esteig_ksp_type cg \ 623 -mg_levels_esteig_ksp_max_it 10 \ 624 -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \ 625 -mg_levels_pc_type jacobi 626 test: 627 suffix: 2d_p1_spectral_0 628 requires: triangle fftw !complex 629 args: -dm_plex_box_faces 1,1 -potential_petscspace_degree 1 -dm_refine 6 -spectral -fft_view 630 test: 631 suffix: 2d_p1_spectral_1 632 requires: triangle fftw !complex 633 nsize: 2 634 args: -dm_plex_box_faces 4,4 -potential_petscspace_degree 1 -spectral -fft_view 635 test: 636 suffix: 2d_p1_adj_0 637 requires: triangle 638 args: -potential_petscspace_degree 1 -dm_refine 1 -adjoint -adjoint_petscspace_degree 1 -error_petscspace_degree 0 639 test: 640 nsize: 2 641 requires: kokkos_kernels 642 suffix: kokkos 643 args: -dm_plex_dim 3 -dm_plex_box_faces 2,3,6 -petscpartitioner_type simple -dm_plex_simplex 0 -potential_petscspace_degree 1 \ 644 -dm_refine 0 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -pc_type gamg -pc_gamg_coarse_eq_limit 1000 -pc_gamg_threshold 0.0 \ 645 -pc_gamg_threshold_scale .5 -mg_levels_ksp_type chebyshev -mg_levels_ksp_max_it 2 -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 \ 646 -ksp_monitor -snes_monitor -dm_view -dm_mat_type aijkokkos -dm_vec_type kokkos 647 648 TEST*/ 649