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