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 const PetscInt *points; 281 282 PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "spectral_plane_%" PetscInt_FMT, p)); 283 PetscCall(DMGetLabel(dm, name, &label)); 284 PetscCall(DMLabelGetStratumIS(label, 1, &stratum)); 285 PetscCall(ISGetLocalSize(stratum, &n)); 286 PetscCall(ISGetIndices(stratum, &points)); 287 PetscCall(PetscMalloc2(n, &ray, n, &svals)); 288 for (i = 0; i < n; ++i) { 289 PetscCall(PetscSectionGetOffset(coordSection, points[i], &off)); 290 PetscCall(PetscSectionGetOffset(section, points[i], &offu)); 291 ray[i] = PetscRealPart(coords[off + ((planeDir[p] + 1) % 2)]); 292 svals[i] = array[offu]; 293 } 294 /* Gather the ray data to proc 0 */ 295 if (size > 1) { 296 PetscMPIInt *cnt, *displs, p; 297 298 PetscCall(PetscCalloc2(size, &cnt, size, &displs)); 299 PetscCallMPI(MPI_Gather(&n, 1, MPIU_INT, cnt, 1, MPIU_INT, 0, comm)); 300 for (p = 1; p < size; ++p) displs[p] = displs[p - 1] + cnt[p - 1]; 301 N = displs[size - 1] + cnt[size - 1]; 302 PetscCall(PetscMalloc2(N, &gray, N, &gsvals)); 303 PetscCallMPI(MPI_Gatherv(ray, n, MPIU_REAL, gray, cnt, displs, MPIU_REAL, 0, comm)); 304 PetscCallMPI(MPI_Gatherv(svals, n, MPIU_SCALAR, gsvals, cnt, displs, MPIU_SCALAR, 0, comm)); 305 PetscCall(PetscFree2(cnt, displs)); 306 } else { 307 N = n; 308 gray = ray; 309 gsvals = svals; 310 } 311 if (rank == 0) { 312 /* Sort point along ray */ 313 PetscCall(PetscMalloc2(N, &perm, N, &nperm)); 314 for (i = 0; i < N; ++i) perm[i] = i; 315 PetscCall(PetscSortRealWithPermutation(N, gray, perm)); 316 /* Count duplicates and squish mapping */ 317 nperm[0] = perm[0]; 318 for (i = 1, j = 1; i < N; ++i) { 319 if (PetscAbsReal(gray[perm[i]] - gray[perm[i - 1]]) > PETSC_SMALL) nperm[j++] = perm[i]; 320 } 321 /* Create FFT structs */ 322 PetscCall(MatCreateFFT(PETSC_COMM_SELF, 1, &j, MATFFTW, &F)); 323 PetscCall(MatCreateVecs(F, &x, &y)); 324 PetscCall(PetscObjectSetName((PetscObject)y, name)); 325 PetscCall(VecGetArray(x, &rvals)); 326 for (i = 0, j = 0; i < N; ++i) { 327 if (i > 0 && PetscAbsReal(gray[perm[i]] - gray[perm[i - 1]]) < PETSC_SMALL) continue; 328 rvals[j] = gsvals[nperm[j]]; 329 ++j; 330 } 331 PetscCall(PetscFree2(perm, nperm)); 332 if (size > 1) PetscCall(PetscFree2(gray, gsvals)); 333 PetscCall(VecRestoreArray(x, &rvals)); 334 /* Do FFT along the ray */ 335 PetscCall(MatMult(F, x, y)); 336 /* Chop FFT */ 337 PetscCall(VecFilter(y, PETSC_SMALL)); 338 PetscCall(VecViewFromOptions(x, NULL, "-real_view")); 339 PetscCall(VecViewFromOptions(y, NULL, "-fft_view")); 340 PetscCall(VecDestroy(&x)); 341 PetscCall(VecDestroy(&y)); 342 PetscCall(MatDestroy(&F)); 343 } 344 PetscCall(ISRestoreIndices(stratum, &points)); 345 PetscCall(ISDestroy(&stratum)); 346 PetscCall(PetscFree2(ray, svals)); 347 } 348 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 349 PetscCall(VecRestoreArrayRead(uloc, &array)); 350 PetscCall(DMRestoreLocalVector(dm, &uloc)); 351 PetscFunctionReturn(PETSC_SUCCESS); 352 } 353 354 static PetscErrorCode ComputeAdjoint(Vec u, AppCtx *user) 355 { 356 PetscFunctionBegin; 357 if (!user->adjoint) PetscFunctionReturn(PETSC_SUCCESS); 358 DM dm, dmAdj; 359 SNES snesAdj; 360 Vec uAdj; 361 362 PetscCall(VecGetDM(u, &dm)); 363 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snesAdj)); 364 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)snesAdj, "adjoint_")); 365 PetscCall(DMClone(dm, &dmAdj)); 366 PetscCall(SNESSetDM(snesAdj, dmAdj)); 367 PetscCall(SetupDiscretization(dmAdj, "adjoint", SetupAdjointProblem, user)); 368 PetscCall(DMCreateGlobalVector(dmAdj, &uAdj)); 369 PetscCall(VecSet(uAdj, 0.0)); 370 PetscCall(PetscObjectSetName((PetscObject)uAdj, "adjoint")); 371 PetscCall(DMPlexSetSNESLocalFEM(dmAdj, user, user, user)); 372 PetscCall(SNESSetFromOptions(snesAdj)); 373 PetscCall(SNESSolve(snesAdj, NULL, uAdj)); 374 PetscCall(SNESGetSolution(snesAdj, &uAdj)); 375 PetscCall(VecViewFromOptions(uAdj, NULL, "-adjoint_view")); 376 /* Error representation */ 377 { 378 DM dmErr, dmErrAux, dms[2]; 379 Vec errorEst, errorL2, uErr, uErrLoc, uAdjLoc, uAdjProj; 380 IS *subis; 381 PetscReal errorEstTot, errorL2Norm, errorL2Tot; 382 PetscInt N, i; 383 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {user->homogeneous ? trig_homogeneous_u : trig_inhomogeneous_u}; 384 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}; 385 void *ctxs[1] = {0}; 386 387 ctxs[0] = user; 388 PetscCall(DMClone(dm, &dmErr)); 389 PetscCall(SetupDiscretization(dmErr, "error", SetupErrorProblem, user)); 390 PetscCall(DMGetGlobalVector(dmErr, &errorEst)); 391 PetscCall(DMGetGlobalVector(dmErr, &errorL2)); 392 /* Compute auxiliary data (solution and projection of adjoint solution) */ 393 PetscCall(DMGetLocalVector(dmAdj, &uAdjLoc)); 394 PetscCall(DMGlobalToLocalBegin(dmAdj, uAdj, INSERT_VALUES, uAdjLoc)); 395 PetscCall(DMGlobalToLocalEnd(dmAdj, uAdj, INSERT_VALUES, uAdjLoc)); 396 PetscCall(DMGetGlobalVector(dm, &uAdjProj)); 397 PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, uAdjLoc)); 398 PetscCall(DMProjectField(dm, 0.0, u, identity, INSERT_VALUES, uAdjProj)); 399 PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL)); 400 PetscCall(DMRestoreLocalVector(dmAdj, &uAdjLoc)); 401 /* Attach auxiliary data */ 402 dms[0] = dm; 403 dms[1] = dm; 404 PetscCall(DMCreateSuperDM(dms, 2, &subis, &dmErrAux)); 405 if (0) { 406 PetscSection sec; 407 408 PetscCall(DMGetLocalSection(dms[0], &sec)); 409 PetscCall(PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD)); 410 PetscCall(DMGetLocalSection(dms[1], &sec)); 411 PetscCall(PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD)); 412 PetscCall(DMGetLocalSection(dmErrAux, &sec)); 413 PetscCall(PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD)); 414 } 415 PetscCall(DMViewFromOptions(dmErrAux, NULL, "-dm_err_view")); 416 PetscCall(ISViewFromOptions(subis[0], NULL, "-super_is_view")); 417 PetscCall(ISViewFromOptions(subis[1], NULL, "-super_is_view")); 418 PetscCall(DMGetGlobalVector(dmErrAux, &uErr)); 419 PetscCall(VecViewFromOptions(u, NULL, "-map_vec_view")); 420 PetscCall(VecViewFromOptions(uAdjProj, NULL, "-map_vec_view")); 421 PetscCall(VecViewFromOptions(uErr, NULL, "-map_vec_view")); 422 PetscCall(VecISCopy(uErr, subis[0], SCATTER_FORWARD, u)); 423 PetscCall(VecISCopy(uErr, subis[1], SCATTER_FORWARD, uAdjProj)); 424 PetscCall(DMRestoreGlobalVector(dm, &uAdjProj)); 425 for (i = 0; i < 2; ++i) PetscCall(ISDestroy(&subis[i])); 426 PetscCall(PetscFree(subis)); 427 PetscCall(DMGetLocalVector(dmErrAux, &uErrLoc)); 428 PetscCall(DMGlobalToLocalBegin(dm, uErr, INSERT_VALUES, uErrLoc)); 429 PetscCall(DMGlobalToLocalEnd(dm, uErr, INSERT_VALUES, uErrLoc)); 430 PetscCall(DMRestoreGlobalVector(dmErrAux, &uErr)); 431 PetscCall(DMSetAuxiliaryVec(dmAdj, NULL, 0, 0, uErrLoc)); 432 /* Compute cellwise error estimate */ 433 PetscCall(VecSet(errorEst, 0.0)); 434 PetscCall(DMPlexComputeCellwiseIntegralFEM(dmAdj, uAdj, errorEst, user)); 435 PetscCall(DMSetAuxiliaryVec(dmAdj, NULL, 0, 0, NULL)); 436 PetscCall(DMRestoreLocalVector(dmErrAux, &uErrLoc)); 437 PetscCall(DMDestroy(&dmErrAux)); 438 /* Plot cellwise error vector */ 439 PetscCall(VecViewFromOptions(errorEst, NULL, "-error_view")); 440 /* Compute ratio of estimate (sum over cells) with actual L_2 error */ 441 PetscCall(DMComputeL2Diff(dm, 0.0, funcs, ctxs, u, &errorL2Norm)); 442 PetscCall(DMPlexComputeL2DiffVec(dm, 0.0, funcs, ctxs, u, errorL2)); 443 PetscCall(VecViewFromOptions(errorL2, NULL, "-l2_error_view")); 444 PetscCall(VecNorm(errorL2, NORM_INFINITY, &errorL2Tot)); 445 PetscCall(VecNorm(errorEst, NORM_INFINITY, &errorEstTot)); 446 PetscCall(VecGetSize(errorEst, &N)); 447 PetscCall(VecPointwiseDivide(errorEst, errorEst, errorL2)); 448 PetscCall(PetscObjectSetName((PetscObject)errorEst, "Error ratio")); 449 PetscCall(VecViewFromOptions(errorEst, NULL, "-error_ratio_view")); 450 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)))); 451 PetscCall(DMRestoreGlobalVector(dmErr, &errorEst)); 452 PetscCall(DMRestoreGlobalVector(dmErr, &errorL2)); 453 PetscCall(DMDestroy(&dmErr)); 454 } 455 PetscCall(DMDestroy(&dmAdj)); 456 PetscCall(VecDestroy(&uAdj)); 457 PetscCall(SNESDestroy(&snesAdj)); 458 PetscFunctionReturn(PETSC_SUCCESS); 459 } 460 461 static PetscErrorCode ErrorView(Vec u, AppCtx *user) 462 { 463 PetscErrorCode (*sol)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *); 464 void *ctx; 465 DM dm; 466 PetscDS ds; 467 PetscReal error; 468 PetscInt N; 469 470 PetscFunctionBegin; 471 if (!user->viewError) PetscFunctionReturn(PETSC_SUCCESS); 472 PetscCall(VecGetDM(u, &dm)); 473 PetscCall(DMGetDS(dm, &ds)); 474 PetscCall(PetscDSGetExactSolution(ds, 0, &sol, &ctx)); 475 PetscCall(VecGetSize(u, &N)); 476 PetscCall(DMComputeL2Diff(dm, 0.0, &sol, &ctx, u, &error)); 477 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "N: %" PetscInt_FMT " L2 error: %g\n", N, (double)error)); 478 PetscFunctionReturn(PETSC_SUCCESS); 479 } 480 481 int main(int argc, char **argv) 482 { 483 DM dm; /* Problem specification */ 484 SNES snes; /* Nonlinear solver */ 485 Vec u; /* Solutions */ 486 AppCtx user; /* User-defined work context */ 487 PetscInt planeDir[2] = {0, 1}; 488 PetscReal planeCoord[2] = {0., 1.}; 489 490 PetscFunctionBeginUser; 491 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 492 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 493 /* Primal system */ 494 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 495 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 496 PetscCall(SNESSetDM(snes, dm)); 497 PetscCall(SetupDiscretization(dm, "potential", SetupPrimalProblem, &user)); 498 PetscCall(DMCreateGlobalVector(dm, &u)); 499 PetscCall(VecSet(u, 0.0)); 500 PetscCall(PetscObjectSetName((PetscObject)u, "potential")); 501 PetscCall(DMPlexSetSNESLocalFEM(dm, &user, &user, &user)); 502 PetscCall(SNESSetFromOptions(snes)); 503 PetscCall(SNESSolve(snes, NULL, u)); 504 PetscCall(SNESGetSolution(snes, &u)); 505 PetscCall(VecViewFromOptions(u, NULL, "-potential_view")); 506 PetscCall(ErrorView(u, &user)); 507 PetscCall(ComputeSpectral(u, 2, planeDir, planeCoord, &user)); 508 PetscCall(ComputeAdjoint(u, &user)); 509 /* Cleanup */ 510 PetscCall(VecDestroy(&u)); 511 PetscCall(SNESDestroy(&snes)); 512 PetscCall(DMDestroy(&dm)); 513 PetscCall(PetscFinalize()); 514 return 0; 515 } 516 517 /*TEST 518 519 test: 520 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9 521 suffix: 2d_p1_conv 522 requires: triangle 523 args: -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 524 test: 525 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9 526 suffix: 2d_p2_conv 527 requires: triangle 528 args: -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 529 test: 530 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9 531 suffix: 2d_p3_conv 532 requires: triangle 533 args: -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2 534 test: 535 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9 536 suffix: 2d_q1_conv 537 args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 538 test: 539 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9 540 suffix: 2d_q2_conv 541 args: -dm_plex_simplex 0 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 542 test: 543 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9 544 suffix: 2d_q3_conv 545 args: -dm_plex_simplex 0 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2 546 test: 547 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9 548 suffix: 2d_q1_ceed_conv 549 requires: libceed 550 args: -dm_plex_use_ceed -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 551 test: 552 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9 553 suffix: 2d_q2_ceed_conv 554 requires: libceed 555 args: -dm_plex_use_ceed -dm_plex_simplex 0 -potential_petscspace_degree 2 -cdm_default_quadrature_order 2 \ 556 -snes_convergence_estimate -convest_num_refine 2 557 test: 558 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9 559 suffix: 2d_q3_ceed_conv 560 requires: libceed 561 args: -dm_plex_use_ceed -dm_plex_simplex 0 -potential_petscspace_degree 3 -cdm_default_quadrature_order 3 \ 562 -snes_convergence_estimate -convest_num_refine 2 563 test: 564 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9 565 suffix: 2d_q1_shear_conv 566 args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2 567 test: 568 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9 569 suffix: 2d_q2_shear_conv 570 args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2 571 test: 572 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9 573 suffix: 2d_q3_shear_conv 574 args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2 575 test: 576 # Using -convest_num_refine 3 we get L_2 convergence rate: 1.7 577 suffix: 3d_p1_conv 578 requires: ctetgen 579 args: -dm_plex_dim 3 -dm_refine 1 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 580 test: 581 # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 2.8 582 suffix: 3d_p2_conv 583 requires: ctetgen 584 args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 1 585 test: 586 # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 4.0 587 suffix: 3d_p3_conv 588 requires: ctetgen 589 args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 1 590 test: 591 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.8 592 suffix: 3d_q1_conv 593 args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 1 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1 594 test: 595 # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.8 596 suffix: 3d_q2_conv 597 args: -dm_plex_dim 3 -dm_plex_simplex 0 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 1 598 test: 599 # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 3.8 600 suffix: 3d_q3_conv 601 args: -dm_plex_dim 3 -dm_plex_simplex 0 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 1 602 test: 603 suffix: 2d_p1_fas_full 604 requires: triangle 605 args: -potential_petscspace_degree 1 -dm_refine_hierarchy 5 \ 606 -snes_max_it 1 -snes_type fas -snes_fas_levels 5 -snes_fas_type full -snes_fas_full_total \ 607 -fas_coarse_snes_monitor -fas_coarse_snes_max_it 1 -fas_coarse_ksp_atol 1.e-13 \ 608 -fas_levels_snes_monitor -fas_levels_snes_max_it 1 -fas_levels_snes_type newtonls \ 609 -fas_levels_pc_type none -fas_levels_ksp_max_it 2 -fas_levels_ksp_converged_maxits -fas_levels_ksp_type chebyshev \ 610 -fas_levels_esteig_ksp_type cg -fas_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -fas_levels_esteig_ksp_max_it 10 611 test: 612 suffix: 2d_p1_fas_full_homogeneous 613 requires: triangle 614 args: -homogeneous -potential_petscspace_degree 1 -dm_refine_hierarchy 5 \ 615 -snes_max_it 1 -snes_type fas -snes_fas_levels 5 -snes_fas_type full \ 616 -fas_coarse_snes_monitor -fas_coarse_snes_max_it 1 -fas_coarse_ksp_atol 1.e-13 \ 617 -fas_levels_snes_monitor -fas_levels_snes_max_it 1 -fas_levels_snes_type newtonls \ 618 -fas_levels_pc_type none -fas_levels_ksp_max_it 2 -fas_levels_ksp_converged_maxits -fas_levels_ksp_type chebyshev \ 619 -fas_levels_esteig_ksp_type cg -fas_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -fas_levels_esteig_ksp_max_it 10 620 621 test: 622 suffix: 2d_p1_scalable 623 requires: triangle 624 args: -potential_petscspace_degree 1 -dm_refine 3 \ 625 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned \ 626 -pc_type gamg -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 \ 627 -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \ 628 -pc_gamg_coarse_eq_limit 1000 \ 629 -pc_gamg_threshold 0.05 \ 630 -pc_gamg_threshold_scale .0 \ 631 -mg_levels_ksp_type chebyshev \ 632 -mg_levels_ksp_max_it 1 \ 633 -mg_levels_pc_type jacobi \ 634 -matptap_via scalable 635 test: 636 suffix: 2d_p1_gmg_vcycle 637 requires: triangle 638 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 639 -ksp_rtol 5e-10 -pc_type mg \ 640 -mg_levels_ksp_max_it 1 \ 641 -mg_levels_esteig_ksp_type cg \ 642 -mg_levels_esteig_ksp_max_it 10 \ 643 -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \ 644 -mg_levels_pc_type jacobi 645 test: 646 suffix: 2d_p1_gmg_fcycle 647 requires: triangle 648 args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 649 -ksp_rtol 5e-10 -pc_type mg -pc_mg_type full \ 650 -mg_levels_ksp_max_it 2 \ 651 -mg_levels_esteig_ksp_type cg \ 652 -mg_levels_esteig_ksp_max_it 10 \ 653 -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \ 654 -mg_levels_pc_type jacobi 655 test: 656 suffix: 2d_p1_gmg_vcycle_adapt 657 requires: triangle 658 args: -petscpartitioner_type simple -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \ 659 -ksp_rtol 5e-10 -pc_type mg -pc_mg_galerkin -pc_mg_adapt_interp_coarse_space harmonic -pc_mg_adapt_interp_n 8 \ 660 -mg_levels_ksp_max_it 1 \ 661 -mg_levels_esteig_ksp_type cg \ 662 -mg_levels_esteig_ksp_max_it 10 \ 663 -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \ 664 -mg_levels_pc_type jacobi 665 test: 666 suffix: 2d_p1_spectral_0 667 requires: triangle fftw !complex 668 args: -dm_plex_box_faces 1,1 -potential_petscspace_degree 1 -dm_refine 6 -spectral -fft_view 669 test: 670 suffix: 2d_p1_spectral_1 671 requires: triangle fftw !complex 672 nsize: 2 673 args: -dm_plex_box_faces 4,4 -potential_petscspace_degree 1 -spectral -fft_view 674 test: 675 suffix: 2d_p1_adj_0 676 requires: triangle 677 args: -potential_petscspace_degree 1 -dm_refine 1 -adjoint -adjoint_petscspace_degree 1 -error_petscspace_degree 0 678 test: 679 nsize: 2 680 requires: kokkos_kernels 681 suffix: kokkos 682 args: -dm_plex_dim 3 -dm_plex_box_faces 2,3,6 -petscpartitioner_type simple -dm_plex_simplex 0 -potential_petscspace_degree 1 \ 683 -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 \ 684 -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 \ 685 -ksp_monitor -snes_monitor -dm_view -dm_mat_type aijkokkos -dm_vec_type kokkos 686 687 TEST*/ 688