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