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