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