1 #include <petsc/private/pcpatchimpl.h> /*I "petscpc.h" I*/ 2 #include <petsc/private/kspimpl.h> /* For ksp->setfromoptionscalled */ 3 #include <petsc/private/dmpleximpl.h> /* For DMPlexComputeJacobian_Patch_Internal() */ 4 #include <petscsf.h> 5 #include <petscbt.h> 6 #include <petscds.h> 7 8 PetscLogEvent PC_Patch_CreatePatches, PC_Patch_ComputeOp, PC_Patch_Solve, PC_Patch_Scatter, PC_Patch_Apply, PC_Patch_Prealloc; 9 10 PETSC_STATIC_INLINE PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format) 11 { 12 PetscErrorCode ierr; 13 14 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 15 ierr = PetscObjectView(obj, viewer);CHKERRQ(ierr); 16 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 17 return(0); 18 } 19 20 static PetscErrorCode PCPatchConstruct_Star(void *vpatch, DM dm, PetscInt point, PetscHSetI ht) 21 { 22 PetscInt starSize; 23 PetscInt *star = NULL, si; 24 PetscErrorCode ierr; 25 26 PetscFunctionBegin; 27 PetscHSetIClear(ht); 28 /* To start with, add the point we care about */ 29 ierr = PetscHSetIAdd(ht, point);CHKERRQ(ierr); 30 /* Loop over all the points that this point connects to */ 31 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 32 for (si = 0; si < starSize*2; si += 2) {ierr = PetscHSetIAdd(ht, star[si]);CHKERRQ(ierr);} 33 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 34 PetscFunctionReturn(0); 35 } 36 37 static PetscErrorCode PCPatchConstruct_Vanka(void *vpatch, DM dm, PetscInt point, PetscHSetI ht) 38 { 39 PC_PATCH *patch = (PC_PATCH *) vpatch; 40 PetscInt starSize; 41 PetscInt *star = NULL; 42 PetscBool shouldIgnore = PETSC_FALSE; 43 PetscInt cStart, cEnd, iStart, iEnd, si; 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 ierr = PetscHSetIClear(ht);CHKERRQ(ierr); 48 /* To start with, add the point we care about */ 49 ierr = PetscHSetIAdd(ht, point);CHKERRQ(ierr); 50 /* Should we ignore any points of a certain dimension? */ 51 if (patch->vankadim >= 0) { 52 shouldIgnore = PETSC_TRUE; 53 ierr = DMPlexGetDepthStratum(dm, patch->vankadim, &iStart, &iEnd);CHKERRQ(ierr); 54 } 55 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 56 /* Loop over all the cells that this point connects to */ 57 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 58 for (si = 0; si < starSize*2; si += 2) { 59 const PetscInt cell = star[si]; 60 PetscInt closureSize; 61 PetscInt *closure = NULL, ci; 62 63 if (cell < cStart || cell >= cEnd) continue; 64 /* now loop over all entities in the closure of that cell */ 65 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 66 for (ci = 0; ci < closureSize*2; ci += 2) { 67 const PetscInt newpoint = closure[ci]; 68 69 /* We've been told to ignore entities of this type.*/ 70 if (shouldIgnore && newpoint >= iStart && newpoint < iEnd) continue; 71 ierr = PetscHSetIAdd(ht, newpoint);CHKERRQ(ierr); 72 } 73 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 74 } 75 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 79 static PetscErrorCode PCPatchConstruct_Pardecomp(void *vpatch, DM dm, PetscInt point, PetscHSetI ht) 80 { 81 PC_PATCH *patch = (PC_PATCH *) vpatch; 82 DMLabel ghost = NULL; 83 const PetscInt *leaves; 84 PetscInt nleaves, pStart, pEnd, loc; 85 PetscBool isFiredrake; 86 DM plex; 87 PetscBool flg; 88 PetscInt starSize; 89 PetscInt *star = NULL; 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 PetscHSetIClear(ht); 94 95 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 96 ierr = DMPlexGetChart(plex, &pStart, &pEnd);CHKERRQ(ierr); 97 98 ierr = DMHasLabel(dm, "pyop2_ghost", &isFiredrake);CHKERRQ(ierr); 99 if (isFiredrake) { 100 ierr = DMGetLabel(dm, "pyop2_ghost", &ghost);CHKERRQ(ierr); 101 ierr = DMLabelCreateIndex(ghost, pStart, pEnd);CHKERRQ(ierr); 102 } else { 103 PetscSF sf; 104 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 105 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 106 nleaves = PetscMax(nleaves, 0); 107 } 108 109 for (PetscInt point = pStart; point < pEnd; ++point) { 110 if (ghost) {ierr = DMLabelHasPoint(ghost, point, &flg);CHKERRQ(ierr);} 111 else {ierr = PetscFindInt(point, nleaves, leaves, &loc);CHKERRQ(ierr); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;} 112 /* Not an owned entity, don't make a cell patch. */ 113 if (flg) continue; 114 ierr = PetscHSetIAdd(ht, point);CHKERRQ(ierr); 115 } 116 117 /* Now build the overlap for the patch */ 118 for (PetscInt overlapi = 0; overlapi < patch->pardecomp_overlap; ++overlapi) { 119 PetscInt index = 0; 120 PetscInt *htpoints = NULL; 121 PetscInt htsize; 122 123 ierr = PetscHSetIGetSize(ht, &htsize);CHKERRQ(ierr); 124 ierr = PetscMalloc1(htsize, &htpoints);CHKERRQ(ierr); 125 ierr = PetscHSetIGetElems(ht, &index, htpoints);CHKERRQ(ierr); 126 127 for (PetscInt i = 0; i < htsize; ++i) { 128 point = htpoints[i]; 129 130 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 131 for (PetscInt si = 0; si < starSize*2; si += 2) { 132 const PetscInt starp = star[si]; 133 PetscInt closureSize; 134 PetscInt *closure = NULL, ci; 135 136 /* now loop over all entities in the closure of starp */ 137 ierr = DMPlexGetTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 138 for (PetscInt ci = 0; ci < closureSize*2; ci += 2) { 139 const PetscInt closstarp = closure[ci]; 140 ierr = PetscHSetIAdd(ht, closstarp);CHKERRQ(ierr); 141 } 142 ierr = DMPlexRestoreTransitiveClosure(dm, starp, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 143 } 144 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 145 } 146 ierr = PetscFree(htpoints);CHKERRQ(ierr); 147 } 148 149 PetscFunctionReturn(0); 150 } 151 152 /* The user's already set the patches in patch->userIS. Build the hash tables */ 153 static PetscErrorCode PCPatchConstruct_User(void *vpatch, DM dm, PetscInt point, PetscHSetI ht) 154 { 155 PC_PATCH *patch = (PC_PATCH *) vpatch; 156 IS patchis = patch->userIS[point]; 157 PetscInt n; 158 const PetscInt *patchdata; 159 PetscInt pStart, pEnd, i; 160 PetscErrorCode ierr; 161 162 PetscFunctionBegin; 163 ierr = PetscHSetIClear(ht);CHKERRQ(ierr); 164 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 165 ierr = ISGetLocalSize(patchis, &n);CHKERRQ(ierr); 166 ierr = ISGetIndices(patchis, &patchdata);CHKERRQ(ierr); 167 for (i = 0; i < n; ++i) { 168 const PetscInt ownedpoint = patchdata[i]; 169 170 if (ownedpoint < pStart || ownedpoint >= pEnd) { 171 SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D was not in [%D, %D)", ownedpoint, pStart, pEnd); 172 } 173 ierr = PetscHSetIAdd(ht, ownedpoint);CHKERRQ(ierr); 174 } 175 ierr = ISRestoreIndices(patchis, &patchdata);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178 179 static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs) 180 { 181 PC_PATCH *patch = (PC_PATCH *) pc->data; 182 PetscInt i; 183 PetscErrorCode ierr; 184 185 PetscFunctionBegin; 186 if (n == 1 && bs[0] == 1) { 187 patch->defaultSF = sf[0]; 188 ierr = PetscObjectReference((PetscObject) patch->defaultSF);CHKERRQ(ierr); 189 } else { 190 PetscInt allRoots = 0, allLeaves = 0; 191 PetscInt leafOffset = 0; 192 PetscInt *ilocal = NULL; 193 PetscSFNode *iremote = NULL; 194 PetscInt *remoteOffsets = NULL; 195 PetscInt index = 0; 196 PetscHMapI rankToIndex; 197 PetscInt numRanks = 0; 198 PetscSFNode *remote = NULL; 199 PetscSF rankSF; 200 PetscInt *ranks = NULL; 201 PetscInt *offsets = NULL; 202 MPI_Datatype contig; 203 PetscHSetI ranksUniq; 204 205 /* First figure out how many dofs there are in the concatenated numbering. 206 * allRoots: number of owned global dofs; 207 * allLeaves: number of visible dofs (global + ghosted). 208 */ 209 for (i = 0; i < n; ++i) { 210 PetscInt nroots, nleaves; 211 212 ierr = PetscSFGetGraph(sf[i], &nroots, &nleaves, NULL, NULL);CHKERRQ(ierr); 213 allRoots += nroots * bs[i]; 214 allLeaves += nleaves * bs[i]; 215 } 216 ierr = PetscMalloc1(allLeaves, &ilocal);CHKERRQ(ierr); 217 ierr = PetscMalloc1(allLeaves, &iremote);CHKERRQ(ierr); 218 /* Now build an SF that just contains process connectivity. */ 219 ierr = PetscHSetICreate(&ranksUniq);CHKERRQ(ierr); 220 for (i = 0; i < n; ++i) { 221 const PetscMPIInt *ranks = NULL; 222 PetscInt nranks, j; 223 224 ierr = PetscSFSetUp(sf[i]);CHKERRQ(ierr); 225 ierr = PetscSFGetRanks(sf[i], &nranks, &ranks, NULL, NULL, NULL);CHKERRQ(ierr); 226 /* These are all the ranks who communicate with me. */ 227 for (j = 0; j < nranks; ++j) { 228 ierr = PetscHSetIAdd(ranksUniq, (PetscInt) ranks[j]);CHKERRQ(ierr); 229 } 230 } 231 ierr = PetscHSetIGetSize(ranksUniq, &numRanks);CHKERRQ(ierr); 232 ierr = PetscMalloc1(numRanks, &remote);CHKERRQ(ierr); 233 ierr = PetscMalloc1(numRanks, &ranks);CHKERRQ(ierr); 234 ierr = PetscHSetIGetElems(ranksUniq, &index, ranks);CHKERRQ(ierr); 235 236 ierr = PetscHMapICreate(&rankToIndex);CHKERRQ(ierr); 237 for (i = 0; i < numRanks; ++i) { 238 remote[i].rank = ranks[i]; 239 remote[i].index = 0; 240 ierr = PetscHMapISet(rankToIndex, ranks[i], i);CHKERRQ(ierr); 241 } 242 ierr = PetscFree(ranks);CHKERRQ(ierr); 243 ierr = PetscHSetIDestroy(&ranksUniq);CHKERRQ(ierr); 244 ierr = PetscSFCreate(PetscObjectComm((PetscObject) pc), &rankSF);CHKERRQ(ierr); 245 ierr = PetscSFSetGraph(rankSF, 1, numRanks, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);CHKERRQ(ierr); 246 ierr = PetscSFSetUp(rankSF);CHKERRQ(ierr); 247 /* OK, use it to communicate the root offset on the remote 248 * processes for each subspace. */ 249 ierr = PetscMalloc1(n, &offsets);CHKERRQ(ierr); 250 ierr = PetscMalloc1(n*numRanks, &remoteOffsets);CHKERRQ(ierr); 251 252 offsets[0] = 0; 253 for (i = 1; i < n; ++i) { 254 PetscInt nroots; 255 256 ierr = PetscSFGetGraph(sf[i-1], &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 257 offsets[i] = offsets[i-1] + nroots*bs[i-1]; 258 } 259 /* Offsets are the offsets on the current process of the 260 * global dof numbering for the subspaces. */ 261 ierr = MPI_Type_contiguous(n, MPIU_INT, &contig);CHKERRQ(ierr); 262 ierr = MPI_Type_commit(&contig);CHKERRQ(ierr); 263 264 ierr = PetscSFBcastBegin(rankSF, contig, offsets, remoteOffsets);CHKERRQ(ierr); 265 ierr = PetscSFBcastEnd(rankSF, contig, offsets, remoteOffsets);CHKERRQ(ierr); 266 ierr = MPI_Type_free(&contig);CHKERRQ(ierr); 267 ierr = PetscFree(offsets);CHKERRQ(ierr); 268 ierr = PetscSFDestroy(&rankSF);CHKERRQ(ierr); 269 /* Now remoteOffsets contains the offsets on the remote 270 * processes who communicate with me. So now we can 271 * concatenate the list of SFs into a single one. */ 272 index = 0; 273 for (i = 0; i < n; ++i) { 274 const PetscSFNode *remote = NULL; 275 const PetscInt *local = NULL; 276 PetscInt nroots, nleaves, j; 277 278 ierr = PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote);CHKERRQ(ierr); 279 for (j = 0; j < nleaves; ++j) { 280 PetscInt rank = remote[j].rank; 281 PetscInt idx, rootOffset, k; 282 283 ierr = PetscHMapIGet(rankToIndex, rank, &idx);CHKERRQ(ierr); 284 if (idx == -1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Didn't find rank, huh?"); 285 /* Offset on given rank for ith subspace */ 286 rootOffset = remoteOffsets[n*idx + i]; 287 for (k = 0; k < bs[i]; ++k) { 288 ilocal[index] = (local ? local[j] : j)*bs[i] + k + leafOffset; 289 iremote[index].rank = remote[j].rank; 290 iremote[index].index = remote[j].index*bs[i] + k + rootOffset; 291 ++index; 292 } 293 } 294 leafOffset += nleaves * bs[i]; 295 } 296 ierr = PetscHMapIDestroy(&rankToIndex);CHKERRQ(ierr); 297 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 298 ierr = PetscSFCreate(PetscObjectComm((PetscObject)pc), &patch->defaultSF);CHKERRQ(ierr); 299 ierr = PetscSFSetGraph(patch->defaultSF, allRoots, allLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);CHKERRQ(ierr); 300 } 301 PetscFunctionReturn(0); 302 } 303 304 /* TODO: Docs */ 305 PetscErrorCode PCPatchSetIgnoreDim(PC pc, PetscInt dim) 306 { 307 PC_PATCH *patch = (PC_PATCH *) pc->data; 308 PetscFunctionBegin; 309 patch->ignoredim = dim; 310 PetscFunctionReturn(0); 311 } 312 313 /* TODO: Docs */ 314 PetscErrorCode PCPatchGetIgnoreDim(PC pc, PetscInt *dim) 315 { 316 PC_PATCH *patch = (PC_PATCH *) pc->data; 317 PetscFunctionBegin; 318 *dim = patch->ignoredim; 319 PetscFunctionReturn(0); 320 } 321 322 /* TODO: Docs */ 323 PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg) 324 { 325 PC_PATCH *patch = (PC_PATCH *) pc->data; 326 PetscFunctionBegin; 327 patch->save_operators = flg; 328 PetscFunctionReturn(0); 329 } 330 331 /* TODO: Docs */ 332 PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg) 333 { 334 PC_PATCH *patch = (PC_PATCH *) pc->data; 335 PetscFunctionBegin; 336 *flg = patch->save_operators; 337 PetscFunctionReturn(0); 338 } 339 340 /* TODO: Docs */ 341 PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg) 342 { 343 PC_PATCH *patch = (PC_PATCH *) pc->data; 344 PetscFunctionBegin; 345 patch->partition_of_unity = flg; 346 PetscFunctionReturn(0); 347 } 348 349 /* TODO: Docs */ 350 PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg) 351 { 352 PC_PATCH *patch = (PC_PATCH *) pc->data; 353 PetscFunctionBegin; 354 *flg = patch->partition_of_unity; 355 PetscFunctionReturn(0); 356 } 357 358 /* TODO: Docs */ 359 PetscErrorCode PCPatchSetLocalComposition(PC pc, PCCompositeType type) 360 { 361 PC_PATCH *patch = (PC_PATCH *) pc->data; 362 PetscFunctionBegin; 363 if (type != PC_COMPOSITE_ADDITIVE && type != PC_COMPOSITE_MULTIPLICATIVE) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Only supports additive or multiplicative as the local type"); 364 patch->local_composition_type = type; 365 PetscFunctionReturn(0); 366 } 367 368 /* TODO: Docs */ 369 PetscErrorCode PCPatchGetLocalComposition(PC pc, PCCompositeType *type) 370 { 371 PC_PATCH *patch = (PC_PATCH *) pc->data; 372 PetscFunctionBegin; 373 *type = patch->local_composition_type; 374 PetscFunctionReturn(0); 375 } 376 377 /* TODO: Docs */ 378 PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type) 379 { 380 PC_PATCH *patch = (PC_PATCH *) pc->data; 381 PetscErrorCode ierr; 382 383 PetscFunctionBegin; 384 if (patch->sub_mat_type) {ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr);} 385 ierr = PetscStrallocpy(sub_mat_type, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 386 PetscFunctionReturn(0); 387 } 388 389 /* TODO: Docs */ 390 PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type) 391 { 392 PC_PATCH *patch = (PC_PATCH *) pc->data; 393 PetscFunctionBegin; 394 *sub_mat_type = patch->sub_mat_type; 395 PetscFunctionReturn(0); 396 } 397 398 /* TODO: Docs */ 399 PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering) 400 { 401 PC_PATCH *patch = (PC_PATCH *) pc->data; 402 PetscErrorCode ierr; 403 404 PetscFunctionBegin; 405 patch->cellNumbering = cellNumbering; 406 ierr = PetscObjectReference((PetscObject) cellNumbering);CHKERRQ(ierr); 407 PetscFunctionReturn(0); 408 } 409 410 /* TODO: Docs */ 411 PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering) 412 { 413 PC_PATCH *patch = (PC_PATCH *) pc->data; 414 PetscFunctionBegin; 415 *cellNumbering = patch->cellNumbering; 416 PetscFunctionReturn(0); 417 } 418 419 /* TODO: Docs */ 420 PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx) 421 { 422 PC_PATCH *patch = (PC_PATCH *) pc->data; 423 424 PetscFunctionBegin; 425 patch->ctype = ctype; 426 switch (ctype) { 427 case PC_PATCH_STAR: 428 patch->user_patches = PETSC_FALSE; 429 patch->patchconstructop = PCPatchConstruct_Star; 430 break; 431 case PC_PATCH_VANKA: 432 patch->user_patches = PETSC_FALSE; 433 patch->patchconstructop = PCPatchConstruct_Vanka; 434 break; 435 case PC_PATCH_PARDECOMP: 436 patch->user_patches = PETSC_FALSE; 437 patch->patchconstructop = PCPatchConstruct_Pardecomp; 438 break; 439 case PC_PATCH_USER: 440 case PC_PATCH_PYTHON: 441 patch->user_patches = PETSC_TRUE; 442 patch->patchconstructop = PCPatchConstruct_User; 443 if (func) { 444 patch->userpatchconstructionop = func; 445 patch->userpatchconstructctx = ctx; 446 } 447 break; 448 default: 449 SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype); 450 } 451 PetscFunctionReturn(0); 452 } 453 454 /* TODO: Docs */ 455 PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx) 456 { 457 PC_PATCH *patch = (PC_PATCH *) pc->data; 458 459 PetscFunctionBegin; 460 *ctype = patch->ctype; 461 switch (patch->ctype) { 462 case PC_PATCH_STAR: 463 case PC_PATCH_VANKA: 464 case PC_PATCH_PARDECOMP: 465 break; 466 case PC_PATCH_USER: 467 case PC_PATCH_PYTHON: 468 *func = patch->userpatchconstructionop; 469 *ctx = patch->userpatchconstructctx; 470 break; 471 default: 472 SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype); 473 } 474 PetscFunctionReturn(0); 475 } 476 477 /* TODO: Docs */ 478 PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, 479 const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 480 { 481 PC_PATCH *patch = (PC_PATCH *) pc->data; 482 DM dm; 483 PetscSF *sfs; 484 PetscInt cStart, cEnd, i, j; 485 PetscErrorCode ierr; 486 487 PetscFunctionBegin; 488 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 489 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 490 ierr = PetscMalloc1(nsubspaces, &sfs);CHKERRQ(ierr); 491 ierr = PetscMalloc1(nsubspaces, &patch->dofSection);CHKERRQ(ierr); 492 ierr = PetscMalloc1(nsubspaces, &patch->bs);CHKERRQ(ierr); 493 ierr = PetscMalloc1(nsubspaces, &patch->nodesPerCell);CHKERRQ(ierr); 494 ierr = PetscMalloc1(nsubspaces, &patch->cellNodeMap);CHKERRQ(ierr); 495 ierr = PetscMalloc1(nsubspaces+1, &patch->subspaceOffsets);CHKERRQ(ierr); 496 497 patch->nsubspaces = nsubspaces; 498 patch->totalDofsPerCell = 0; 499 for (i = 0; i < nsubspaces; ++i) { 500 ierr = DMGetDefaultSection(dms[i], &patch->dofSection[i]);CHKERRQ(ierr); 501 ierr = PetscObjectReference((PetscObject) patch->dofSection[i]);CHKERRQ(ierr); 502 ierr = DMGetDefaultSF(dms[i], &sfs[i]);CHKERRQ(ierr); 503 patch->bs[i] = bs[i]; 504 patch->nodesPerCell[i] = nodesPerCell[i]; 505 patch->totalDofsPerCell += nodesPerCell[i]*bs[i]; 506 ierr = PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);CHKERRQ(ierr); 507 for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j]; 508 patch->subspaceOffsets[i] = subspaceOffsets[i]; 509 } 510 ierr = PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs);CHKERRQ(ierr); 511 ierr = PetscFree(sfs);CHKERRQ(ierr); 512 513 patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces]; 514 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);CHKERRQ(ierr); 515 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);CHKERRQ(ierr); 516 PetscFunctionReturn(0); 517 } 518 519 /* TODO: Docs */ 520 PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 521 { 522 PC_PATCH *patch = (PC_PATCH *) pc->data; 523 PetscInt cStart, cEnd, i, j; 524 PetscErrorCode ierr; 525 526 PetscFunctionBegin; 527 patch->combined = PETSC_TRUE; 528 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 529 ierr = DMGetNumFields(dm, &patch->nsubspaces);CHKERRQ(ierr); 530 ierr = PetscCalloc1(patch->nsubspaces, &patch->dofSection);CHKERRQ(ierr); 531 ierr = PetscMalloc1(patch->nsubspaces, &patch->bs);CHKERRQ(ierr); 532 ierr = PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell);CHKERRQ(ierr); 533 ierr = PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap);CHKERRQ(ierr); 534 ierr = PetscCalloc1(patch->nsubspaces+1, &patch->subspaceOffsets);CHKERRQ(ierr); 535 ierr = DMGetDefaultSection(dm, &patch->dofSection[0]);CHKERRQ(ierr); 536 ierr = PetscObjectReference((PetscObject) patch->dofSection[0]);CHKERRQ(ierr); 537 ierr = PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]);CHKERRQ(ierr); 538 patch->totalDofsPerCell = 0; 539 for (i = 0; i < patch->nsubspaces; ++i) { 540 patch->bs[i] = 1; 541 patch->nodesPerCell[i] = nodesPerCell[i]; 542 patch->totalDofsPerCell += nodesPerCell[i]; 543 ierr = PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);CHKERRQ(ierr); 544 for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j]; 545 } 546 ierr = DMGetDefaultSF(dm, &patch->defaultSF);CHKERRQ(ierr); 547 ierr = PetscObjectReference((PetscObject) patch->defaultSF);CHKERRQ(ierr); 548 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);CHKERRQ(ierr); 549 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);CHKERRQ(ierr); 550 PetscFunctionReturn(0); 551 } 552 553 /*@C 554 555 PCPatchSetComputeFunction - Set the callback used to compute patch residuals 556 557 Input Parameters: 558 + pc - The PC 559 . func - The callback 560 - ctx - The user context 561 562 Level: advanced 563 564 Note: 565 The callback has signature: 566 + usercomputef(pc, point, x, f, cellIS, n, u, ctx) 567 + pc - The PC 568 + point - The point 569 + x - The input solution (not used in linear problems) 570 + f - The patch residual vector 571 + cellIS - An array of the cell numbers 572 + n - The size of g2l 573 + g2l - The global to local dof translation table 574 + ctx - The user context 575 and can assume that the matrix entries have been set to zero before the call. 576 577 .seealso: PCPatchSetComputeOperator(), PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo() 578 @*/ 579 PetscErrorCode PCPatchSetComputeFunction(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx) 580 { 581 PC_PATCH *patch = (PC_PATCH *) pc->data; 582 583 PetscFunctionBegin; 584 patch->usercomputef = func; 585 patch->usercomputefctx = ctx; 586 PetscFunctionReturn(0); 587 } 588 589 /*@C 590 591 PCPatchSetComputeOperator - Set the callback used to compute patch matrices 592 593 Input Parameters: 594 + pc - The PC 595 . func - The callback 596 - ctx - The user context 597 598 Level: advanced 599 600 Note: 601 The callback has signature: 602 + usercomputeop(pc, point, x, mat, cellIS, n, u, ctx) 603 + pc - The PC 604 + point - The point 605 + x - The input solution (not used in linear problems) 606 + mat - The patch matrix 607 + cellIS - An array of the cell numbers 608 + n - The size of g2l 609 + g2l - The global to local dof translation table 610 + ctx - The user context 611 and can assume that the matrix entries have been set to zero before the call. 612 613 .seealso: PCPatchGetComputeOperator(), PCPatchSetComputeFunction(), PCPatchSetDiscretisationInfo() 614 @*/ 615 PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx) 616 { 617 PC_PATCH *patch = (PC_PATCH *) pc->data; 618 619 PetscFunctionBegin; 620 patch->usercomputeop = func; 621 patch->usercomputeopctx = ctx; 622 PetscFunctionReturn(0); 623 } 624 625 /* On entry, ht contains the topological entities whose dofs we are responsible for solving for; 626 on exit, cht contains all the topological entities we need to compute their residuals. 627 In full generality this should incorporate knowledge of the sparsity pattern of the matrix; 628 here we assume a standard FE sparsity pattern.*/ 629 /* TODO: Use DMPlexGetAdjacency() */ 630 static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht) 631 { 632 DM dm; 633 PetscHashIter hi; 634 PetscInt point; 635 PetscInt *star = NULL, *closure = NULL; 636 PetscInt ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci; 637 PetscErrorCode ierr; 638 639 PetscFunctionBegin; 640 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 641 ierr = PCPatchGetIgnoreDim(pc, &ignoredim);CHKERRQ(ierr); 642 if (ignoredim >= 0) {ierr = DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd);CHKERRQ(ierr);} 643 ierr = PetscHSetIClear(cht);CHKERRQ(ierr); 644 PetscHashIterBegin(ht, hi); 645 while (!PetscHashIterAtEnd(ht, hi)) { 646 647 PetscHashIterGetKey(ht, hi, point); 648 PetscHashIterNext(ht, hi); 649 650 /* Loop over all the cells that this point connects to */ 651 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 652 for (si = 0; si < starSize*2; si += 2) { 653 const PetscInt ownedpoint = star[si]; 654 /* TODO Check for point in cht before running through closure again */ 655 /* now loop over all entities in the closure of that cell */ 656 ierr = DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 657 for (ci = 0; ci < closureSize*2; ci += 2) { 658 const PetscInt seenpoint = closure[ci]; 659 if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue; 660 ierr = PetscHSetIAdd(cht, seenpoint);CHKERRQ(ierr); 661 } 662 ierr = DMPlexRestoreTransitiveClosure(dm, ownedpoint, PETSC_TRUE, NULL, &closure);CHKERRQ(ierr); 663 } 664 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, NULL, &star);CHKERRQ(ierr); 665 } 666 PetscFunctionReturn(0); 667 } 668 669 static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off) 670 { 671 PetscErrorCode ierr; 672 673 PetscFunctionBegin; 674 if (combined) { 675 if (f < 0) { 676 if (dof) {ierr = PetscSectionGetDof(dofSection[0], p, dof);CHKERRQ(ierr);} 677 if (off) {ierr = PetscSectionGetOffset(dofSection[0], p, off);CHKERRQ(ierr);} 678 } else { 679 if (dof) {ierr = PetscSectionGetFieldDof(dofSection[0], p, f, dof);CHKERRQ(ierr);} 680 if (off) {ierr = PetscSectionGetFieldOffset(dofSection[0], p, f, off);CHKERRQ(ierr);} 681 } 682 } else { 683 if (f < 0) { 684 PC_PATCH *patch = (PC_PATCH *) pc->data; 685 PetscInt fdof, g; 686 687 if (dof) { 688 *dof = 0; 689 for (g = 0; g < patch->nsubspaces; ++g) { 690 ierr = PetscSectionGetDof(dofSection[g], p, &fdof);CHKERRQ(ierr); 691 *dof += fdof; 692 } 693 } 694 if (off) { 695 *off = 0; 696 for (g = 0; g < patch->nsubspaces; ++g) { 697 ierr = PetscSectionGetOffset(dofSection[g], p, &fdof);CHKERRQ(ierr); 698 *off += fdof; 699 } 700 } 701 } else { 702 if (dof) {ierr = PetscSectionGetDof(dofSection[f], p, dof);CHKERRQ(ierr);} 703 if (off) {ierr = PetscSectionGetOffset(dofSection[f], p, off);CHKERRQ(ierr);} 704 } 705 } 706 PetscFunctionReturn(0); 707 } 708 709 /* Given a hash table with a set of topological entities (pts), compute the degrees of 710 freedom in global concatenated numbering on those entities. 711 For Vanka smoothing, this needs to do something special: ignore dofs of the 712 constraint subspace on entities that aren't the base entity we're building the patch 713 around. */ 714 static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI* subspaces_to_exclude) 715 { 716 PC_PATCH *patch = (PC_PATCH *) pc->data; 717 PetscHashIter hi; 718 PetscInt ldof, loff; 719 PetscInt k, p; 720 PetscErrorCode ierr; 721 722 PetscFunctionBegin; 723 ierr = PetscHSetIClear(dofs);CHKERRQ(ierr); 724 for (k = 0; k < patch->nsubspaces; ++k) { 725 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 726 PetscInt bs = patch->bs[k]; 727 PetscInt j, l; 728 729 if (subspaces_to_exclude != NULL) { 730 PetscBool should_exclude_k = PETSC_FALSE; 731 PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k); 732 if (should_exclude_k) { 733 /* only get this subspace dofs at the base entity, not any others */ 734 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff);CHKERRQ(ierr); 735 if (0 == ldof) continue; 736 for (j = loff; j < ldof + loff; ++j) { 737 for (l = 0; l < bs; ++l) { 738 PetscInt dof = bs*j + l + subspaceOffset; 739 ierr = PetscHSetIAdd(dofs, dof);CHKERRQ(ierr); 740 } 741 } 742 continue; /* skip the other dofs of this subspace */ 743 } 744 } 745 746 PetscHashIterBegin(pts, hi); 747 while (!PetscHashIterAtEnd(pts, hi)) { 748 PetscHashIterGetKey(pts, hi, p); 749 PetscHashIterNext(pts, hi); 750 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);CHKERRQ(ierr); 751 if (0 == ldof) continue; 752 for (j = loff; j < ldof + loff; ++j) { 753 for (l = 0; l < bs; ++l) { 754 PetscInt dof = bs*j + l + subspaceOffset; 755 ierr = PetscHSetIAdd(dofs, dof);CHKERRQ(ierr); 756 } 757 } 758 } 759 } 760 PetscFunctionReturn(0); 761 } 762 763 /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */ 764 static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C) 765 { 766 PetscHashIter hi; 767 PetscInt key; 768 PetscBool flg; 769 PetscErrorCode ierr; 770 771 PetscFunctionBegin; 772 ierr = PetscHSetIClear(C);CHKERRQ(ierr); 773 PetscHashIterBegin(B, hi); 774 while (!PetscHashIterAtEnd(B, hi)) { 775 PetscHashIterGetKey(B, hi, key); 776 PetscHashIterNext(B, hi); 777 ierr = PetscHSetIHas(A, key, &flg);CHKERRQ(ierr); 778 if (!flg) {ierr = PetscHSetIAdd(C, key);CHKERRQ(ierr);} 779 } 780 PetscFunctionReturn(0); 781 } 782 783 /* 784 * PCPatchCreateCellPatches - create patches. 785 * 786 * Input Parameters: 787 * + dm - The DMPlex object defining the mesh 788 * 789 * Output Parameters: 790 * + cellCounts - Section with counts of cells around each vertex 791 * . cells - IS of the cell point indices of cells in each patch 792 * . pointCounts - Section with counts of cells around each vertex 793 * - point - IS of the cell point indices of cells in each patch 794 */ 795 static PetscErrorCode PCPatchCreateCellPatches(PC pc) 796 { 797 PC_PATCH *patch = (PC_PATCH *) pc->data; 798 DMLabel ghost = NULL; 799 DM dm, plex; 800 PetscHSetI ht, cht; 801 PetscSection cellCounts, pointCounts; 802 PetscInt *cellsArray, *pointsArray; 803 PetscInt numCells, numPoints; 804 const PetscInt *leaves; 805 PetscInt nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, v; 806 PetscBool isFiredrake; 807 PetscErrorCode ierr; 808 809 PetscFunctionBegin; 810 /* Used to keep track of the cells in the patch. */ 811 ierr = PetscHSetICreate(&ht);CHKERRQ(ierr); 812 ierr = PetscHSetICreate(&cht);CHKERRQ(ierr); 813 814 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 815 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC\n"); 816 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 817 ierr = DMPlexGetChart(plex, &pStart, &pEnd);CHKERRQ(ierr); 818 ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr); 819 820 if (patch->user_patches) { 821 ierr = patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);CHKERRQ(ierr); 822 vStart = 0; vEnd = patch->npatch; 823 } else if (patch->ctype == PC_PATCH_PARDECOMP) { 824 vStart = 0; vEnd = 1; 825 } else if (patch->codim < 0) { 826 if (patch->dim < 0) {ierr = DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);CHKERRQ(ierr);} 827 else {ierr = DMPlexGetDepthStratum(plex, patch->dim, &vStart, &vEnd);CHKERRQ(ierr);} 828 } else {ierr = DMPlexGetHeightStratum(plex, patch->codim, &vStart, &vEnd);CHKERRQ(ierr);} 829 patch->npatch = vEnd - vStart; 830 831 /* These labels mark the owned points. We only create patches around points that this process owns. */ 832 ierr = DMHasLabel(dm, "pyop2_ghost", &isFiredrake);CHKERRQ(ierr); 833 if (isFiredrake) { 834 ierr = DMGetLabel(dm, "pyop2_ghost", &ghost);CHKERRQ(ierr); 835 ierr = DMLabelCreateIndex(ghost, pStart, pEnd);CHKERRQ(ierr); 836 } else { 837 PetscSF sf; 838 839 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 840 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 841 nleaves = PetscMax(nleaves, 0); 842 } 843 844 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts);CHKERRQ(ierr); 845 ierr = PetscObjectSetName((PetscObject) patch->cellCounts, "Patch Cell Layout");CHKERRQ(ierr); 846 cellCounts = patch->cellCounts; 847 ierr = PetscSectionSetChart(cellCounts, vStart, vEnd);CHKERRQ(ierr); 848 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts);CHKERRQ(ierr); 849 ierr = PetscObjectSetName((PetscObject) patch->pointCounts, "Patch Point Layout");CHKERRQ(ierr); 850 pointCounts = patch->pointCounts; 851 ierr = PetscSectionSetChart(pointCounts, vStart, vEnd);CHKERRQ(ierr); 852 /* Count cells and points in the patch surrounding each entity */ 853 for (v = vStart; v < vEnd; ++v) { 854 PetscHashIter hi; 855 PetscInt chtSize, loc = -1; 856 PetscBool flg; 857 858 if (!patch->user_patches && patch->ctype != PC_PATCH_PARDECOMP) { 859 if (ghost) {ierr = DMLabelHasPoint(ghost, v, &flg);CHKERRQ(ierr);} 860 else {ierr = PetscFindInt(v, nleaves, leaves, &loc);CHKERRQ(ierr); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;} 861 /* Not an owned entity, don't make a cell patch. */ 862 if (flg) continue; 863 } 864 865 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 866 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 867 ierr = PetscHSetIGetSize(cht, &chtSize);CHKERRQ(ierr); 868 /* empty patch, continue */ 869 if (chtSize == 0) continue; 870 871 /* safe because size(cht) > 0 from above */ 872 PetscHashIterBegin(cht, hi); 873 while (!PetscHashIterAtEnd(cht, hi)) { 874 PetscInt point, pdof; 875 876 PetscHashIterGetKey(cht, hi, point); 877 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 878 if (pdof) {ierr = PetscSectionAddDof(pointCounts, v, 1);CHKERRQ(ierr);} 879 if (point >= cStart && point < cEnd) {ierr = PetscSectionAddDof(cellCounts, v, 1);CHKERRQ(ierr);} 880 PetscHashIterNext(cht, hi); 881 } 882 } 883 if (isFiredrake) {ierr = DMLabelDestroyIndex(ghost);CHKERRQ(ierr);} 884 885 ierr = PetscSectionSetUp(cellCounts);CHKERRQ(ierr); 886 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 887 ierr = PetscMalloc1(numCells, &cellsArray);CHKERRQ(ierr); 888 ierr = PetscSectionSetUp(pointCounts);CHKERRQ(ierr); 889 ierr = PetscSectionGetStorageSize(pointCounts, &numPoints);CHKERRQ(ierr); 890 ierr = PetscMalloc1(numPoints, &pointsArray);CHKERRQ(ierr); 891 892 /* Now that we know how much space we need, run through again and actually remember the cells. */ 893 for (v = vStart; v < vEnd; v++ ) { 894 PetscHashIter hi; 895 PetscInt dof, off, cdof, coff, pdof, n = 0, cn = 0; 896 897 ierr = PetscSectionGetDof(pointCounts, v, &dof);CHKERRQ(ierr); 898 ierr = PetscSectionGetOffset(pointCounts, v, &off);CHKERRQ(ierr); 899 ierr = PetscSectionGetDof(cellCounts, v, &cdof);CHKERRQ(ierr); 900 ierr = PetscSectionGetOffset(cellCounts, v, &coff);CHKERRQ(ierr); 901 if (dof <= 0) continue; 902 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 903 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 904 PetscHashIterBegin(cht, hi); 905 while (!PetscHashIterAtEnd(cht, hi)) { 906 PetscInt point; 907 908 PetscHashIterGetKey(cht, hi, point); 909 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 910 if (pdof) {pointsArray[off + n++] = point;} 911 if (point >= cStart && point < cEnd) {cellsArray[coff + cn++] = point;} 912 PetscHashIterNext(cht, hi); 913 } 914 if (cn != cdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of cells in patch %D is %D, but should be %D", v, cn, cdof); 915 if (n != dof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of points in patch %D is %D, but should be %D", v, n, dof); 916 } 917 ierr = PetscHSetIDestroy(&ht);CHKERRQ(ierr); 918 ierr = PetscHSetIDestroy(&cht);CHKERRQ(ierr); 919 ierr = DMDestroy(&plex);CHKERRQ(ierr); 920 921 ierr = ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells);CHKERRQ(ierr); 922 ierr = PetscObjectSetName((PetscObject) patch->cells, "Patch Cells");CHKERRQ(ierr); 923 if (patch->viewCells) { 924 ierr = ObjectView((PetscObject) patch->cellCounts, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 925 ierr = ObjectView((PetscObject) patch->cells, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 926 } 927 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points);CHKERRQ(ierr); 928 ierr = PetscObjectSetName((PetscObject) patch->points, "Patch Points");CHKERRQ(ierr); 929 if (patch->viewPoints) { 930 ierr = ObjectView((PetscObject) patch->pointCounts, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 931 ierr = ObjectView((PetscObject) patch->points, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 932 } 933 PetscFunctionReturn(0); 934 } 935 936 /* 937 * PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches 938 * 939 * Input Parameters: 940 * + dm - The DMPlex object defining the mesh 941 * . cellCounts - Section with counts of cells around each vertex 942 * . cells - IS of the cell point indices of cells in each patch 943 * . cellNumbering - Section mapping plex cell points to Firedrake cell indices. 944 * . nodesPerCell - number of nodes per cell. 945 * - cellNodeMap - map from cells to node indices (nodesPerCell * numCells) 946 * 947 * Output Parameters: 948 * + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering 949 * . gtolCounts - Section with counts of dofs per cell patch 950 * - gtol - IS mapping from global dofs to local dofs for each patch. 951 */ 952 static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc) 953 { 954 PC_PATCH *patch = (PC_PATCH *) pc->data; 955 PetscSection cellCounts = patch->cellCounts; 956 PetscSection pointCounts = patch->pointCounts; 957 PetscSection gtolCounts, gtolCountsWithArtificial = NULL, gtolCountsWithAll = NULL; 958 IS cells = patch->cells; 959 IS points = patch->points; 960 PetscSection cellNumbering = patch->cellNumbering; 961 PetscInt Nf = patch->nsubspaces; 962 PetscInt numCells, numPoints; 963 PetscInt numDofs; 964 PetscInt numGlobalDofs, numGlobalDofsWithArtificial, numGlobalDofsWithAll; 965 PetscInt totalDofsPerCell = patch->totalDofsPerCell; 966 PetscInt vStart, vEnd, v; 967 const PetscInt *cellsArray, *pointsArray; 968 PetscInt *newCellsArray = NULL; 969 PetscInt *dofsArray = NULL; 970 PetscInt *dofsArrayWithArtificial = NULL; 971 PetscInt *dofsArrayWithAll = NULL; 972 PetscInt *offsArray = NULL; 973 PetscInt *offsArrayWithArtificial = NULL; 974 PetscInt *offsArrayWithAll = NULL; 975 PetscInt *asmArray = NULL; 976 PetscInt *asmArrayWithArtificial = NULL; 977 PetscInt *asmArrayWithAll = NULL; 978 PetscInt *globalDofsArray = NULL; 979 PetscInt *globalDofsArrayWithArtificial = NULL; 980 PetscInt *globalDofsArrayWithAll = NULL; 981 PetscInt globalIndex = 0; 982 PetscInt key = 0; 983 PetscInt asmKey = 0; 984 DM dm = NULL; 985 const PetscInt *bcNodes = NULL; 986 PetscHMapI ht; 987 PetscHMapI htWithArtificial; 988 PetscHMapI htWithAll; 989 PetscHSetI globalBcs; 990 PetscInt numBcs; 991 PetscHSetI ownedpts, seenpts, owneddofs, seendofs, artificialbcs; 992 PetscInt pStart, pEnd, p, i; 993 char option[PETSC_MAX_PATH_LEN]; 994 PetscBool isNonlinear; 995 PetscErrorCode ierr; 996 997 PetscFunctionBegin; 998 999 ierr = PCGetDM(pc, &dm); CHKERRQ(ierr); 1000 /* dofcounts section is cellcounts section * dofPerCell */ 1001 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 1002 ierr = PetscSectionGetStorageSize(patch->pointCounts, &numPoints);CHKERRQ(ierr); 1003 numDofs = numCells * totalDofsPerCell; 1004 ierr = PetscMalloc1(numDofs, &dofsArray);CHKERRQ(ierr); 1005 ierr = PetscMalloc1(numPoints*Nf, &offsArray);CHKERRQ(ierr); 1006 ierr = PetscMalloc1(numDofs, &asmArray);CHKERRQ(ierr); 1007 ierr = PetscMalloc1(numCells, &newCellsArray);CHKERRQ(ierr); 1008 ierr = PetscSectionGetChart(cellCounts, &vStart, &vEnd);CHKERRQ(ierr); 1009 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);CHKERRQ(ierr); 1010 gtolCounts = patch->gtolCounts; 1011 ierr = PetscSectionSetChart(gtolCounts, vStart, vEnd);CHKERRQ(ierr); 1012 ierr = PetscObjectSetName((PetscObject) patch->gtolCounts, "Patch Global Index Section");CHKERRQ(ierr); 1013 1014 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) 1015 { 1016 ierr = PetscMalloc1(numPoints*Nf, &offsArrayWithArtificial);CHKERRQ(ierr); 1017 ierr = PetscMalloc1(numDofs, &asmArrayWithArtificial);CHKERRQ(ierr); 1018 ierr = PetscMalloc1(numDofs, &dofsArrayWithArtificial);CHKERRQ(ierr); 1019 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial);CHKERRQ(ierr); 1020 gtolCountsWithArtificial = patch->gtolCountsWithArtificial; 1021 ierr = PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd);CHKERRQ(ierr); 1022 ierr = PetscObjectSetName((PetscObject) patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs");CHKERRQ(ierr); 1023 } 1024 1025 isNonlinear = patch->isNonlinear; 1026 if(isNonlinear) 1027 { 1028 ierr = PetscMalloc1(numPoints*Nf, &offsArrayWithAll);CHKERRQ(ierr); 1029 ierr = PetscMalloc1(numDofs, &asmArrayWithAll);CHKERRQ(ierr); 1030 ierr = PetscMalloc1(numDofs, &dofsArrayWithAll);CHKERRQ(ierr); 1031 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll);CHKERRQ(ierr); 1032 gtolCountsWithAll = patch->gtolCountsWithAll; 1033 ierr = PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd);CHKERRQ(ierr); 1034 ierr = PetscObjectSetName((PetscObject) patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs");CHKERRQ(ierr); 1035 } 1036 1037 /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet 1038 conditions */ 1039 ierr = PetscHSetICreate(&globalBcs);CHKERRQ(ierr); 1040 ierr = ISGetIndices(patch->ghostBcNodes, &bcNodes); CHKERRQ(ierr); 1041 ierr = ISGetSize(patch->ghostBcNodes, &numBcs); CHKERRQ(ierr); 1042 for (i = 0; i < numBcs; ++i) { 1043 ierr = PetscHSetIAdd(globalBcs, bcNodes[i]);CHKERRQ(ierr); /* these are already in concatenated numbering */ 1044 } 1045 ierr = ISRestoreIndices(patch->ghostBcNodes, &bcNodes); CHKERRQ(ierr); 1046 ierr = ISDestroy(&patch->ghostBcNodes); CHKERRQ(ierr); /* memory optimisation */ 1047 1048 /* Hash tables for artificial BC construction */ 1049 ierr = PetscHSetICreate(&ownedpts);CHKERRQ(ierr); 1050 ierr = PetscHSetICreate(&seenpts);CHKERRQ(ierr); 1051 ierr = PetscHSetICreate(&owneddofs);CHKERRQ(ierr); 1052 ierr = PetscHSetICreate(&seendofs);CHKERRQ(ierr); 1053 ierr = PetscHSetICreate(&artificialbcs);CHKERRQ(ierr); 1054 1055 ierr = ISGetIndices(cells, &cellsArray);CHKERRQ(ierr); 1056 ierr = ISGetIndices(points, &pointsArray);CHKERRQ(ierr); 1057 ierr = PetscHMapICreate(&ht);CHKERRQ(ierr); 1058 ierr = PetscHMapICreate(&htWithArtificial);CHKERRQ(ierr); 1059 ierr = PetscHMapICreate(&htWithAll);CHKERRQ(ierr); 1060 for (v = vStart; v < vEnd; ++v) { 1061 PetscInt localIndex = 0; 1062 PetscInt localIndexWithArtificial = 0; 1063 PetscInt localIndexWithAll = 0; 1064 PetscInt dof, off, i, j, k, l; 1065 1066 ierr = PetscHMapIClear(ht);CHKERRQ(ierr); 1067 ierr = PetscHMapIClear(htWithArtificial);CHKERRQ(ierr); 1068 ierr = PetscHMapIClear(htWithAll);CHKERRQ(ierr); 1069 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 1070 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 1071 if (dof <= 0) continue; 1072 1073 /* Calculate the global numbers of the artificial BC dofs here first */ 1074 ierr = patch->patchconstructop((void*)patch, dm, v, ownedpts); CHKERRQ(ierr); 1075 ierr = PCPatchCompleteCellPatch(pc, ownedpts, seenpts); CHKERRQ(ierr); 1076 ierr = PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude); CHKERRQ(ierr); 1077 ierr = PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL); CHKERRQ(ierr); 1078 ierr = PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs); CHKERRQ(ierr); 1079 if (patch->viewPatches) { 1080 PetscHSetI globalbcdofs; 1081 PetscHashIter hi; 1082 MPI_Comm comm = PetscObjectComm((PetscObject)pc); 1083 1084 ierr = PetscHSetICreate(&globalbcdofs);CHKERRQ(ierr); 1085 ierr = PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v); CHKERRQ(ierr); 1086 PetscHashIterBegin(owneddofs, hi); 1087 while (!PetscHashIterAtEnd(owneddofs, hi)) { 1088 PetscInt globalDof; 1089 1090 PetscHashIterGetKey(owneddofs, hi, globalDof); 1091 PetscHashIterNext(owneddofs, hi); 1092 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 1093 } 1094 ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr); 1095 ierr = PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v); CHKERRQ(ierr); 1096 PetscHashIterBegin(seendofs, hi); 1097 while (!PetscHashIterAtEnd(seendofs, hi)) { 1098 PetscInt globalDof; 1099 PetscBool flg; 1100 1101 PetscHashIterGetKey(seendofs, hi, globalDof); 1102 PetscHashIterNext(seendofs, hi); 1103 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 1104 1105 ierr = PetscHSetIHas(globalBcs, globalDof, &flg);CHKERRQ(ierr); 1106 if (flg) {ierr = PetscHSetIAdd(globalbcdofs, globalDof);CHKERRQ(ierr);} 1107 } 1108 ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr); 1109 ierr = PetscSynchronizedPrintf(comm, "Patch %d: global BCs:\n", v); CHKERRQ(ierr); 1110 ierr = PetscHSetIGetSize(globalbcdofs, &numBcs);CHKERRQ(ierr); 1111 if (numBcs > 0) { 1112 PetscHashIterBegin(globalbcdofs, hi); 1113 while (!PetscHashIterAtEnd(globalbcdofs, hi)) { 1114 PetscInt globalDof; 1115 PetscHashIterGetKey(globalbcdofs, hi, globalDof); 1116 PetscHashIterNext(globalbcdofs, hi); 1117 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof);CHKERRQ(ierr); 1118 } 1119 } 1120 ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr); 1121 ierr = PetscSynchronizedPrintf(comm, "Patch %d: artificial BCs:\n", v);CHKERRQ(ierr); 1122 ierr = PetscHSetIGetSize(artificialbcs, &numBcs);CHKERRQ(ierr); 1123 if (numBcs > 0) { 1124 PetscHashIterBegin(artificialbcs, hi); 1125 while (!PetscHashIterAtEnd(artificialbcs, hi)) { 1126 PetscInt globalDof; 1127 PetscHashIterGetKey(artificialbcs, hi, globalDof); 1128 PetscHashIterNext(artificialbcs, hi); 1129 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 1130 } 1131 } 1132 ierr = PetscSynchronizedPrintf(comm, "\n\n"); CHKERRQ(ierr); 1133 ierr = PetscHSetIDestroy(&globalbcdofs);CHKERRQ(ierr); 1134 } 1135 for (k = 0; k < patch->nsubspaces; ++k) { 1136 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1137 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1138 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1139 PetscInt bs = patch->bs[k]; 1140 1141 for (i = off; i < off + dof; ++i) { 1142 /* Walk over the cells in this patch. */ 1143 const PetscInt c = cellsArray[i]; 1144 PetscInt cell = c; 1145 1146 /* TODO Change this to an IS */ 1147 if (cellNumbering) { 1148 ierr = PetscSectionGetDof(cellNumbering, c, &cell);CHKERRQ(ierr); 1149 if (cell <= 0) SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %D doesn't appear in cell numbering map", c); 1150 ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr); 1151 } 1152 newCellsArray[i] = cell; 1153 for (j = 0; j < nodesPerCell; ++j) { 1154 /* For each global dof, map it into contiguous local storage. */ 1155 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + subspaceOffset; 1156 /* finally, loop over block size */ 1157 for (l = 0; l < bs; ++l) { 1158 PetscInt localDof; 1159 PetscBool isGlobalBcDof, isArtificialBcDof; 1160 1161 /* first, check if this is either a globally enforced or locally enforced BC dof */ 1162 ierr = PetscHSetIHas(globalBcs, globalDof + l, &isGlobalBcDof);CHKERRQ(ierr); 1163 ierr = PetscHSetIHas(artificialbcs, globalDof + l, &isArtificialBcDof);CHKERRQ(ierr); 1164 1165 /* if it's either, don't ever give it a local dof number */ 1166 if (isGlobalBcDof || isArtificialBcDof) { 1167 dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */ 1168 } else { 1169 ierr = PetscHMapIGet(ht, globalDof + l, &localDof);CHKERRQ(ierr); 1170 if (localDof == -1) { 1171 localDof = localIndex++; 1172 ierr = PetscHMapISet(ht, globalDof + l, localDof);CHKERRQ(ierr); 1173 } 1174 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 1175 /* And store. */ 1176 dofsArray[globalIndex] = localDof; 1177 } 1178 1179 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1180 if (isGlobalBcDof) { 1181 dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */ 1182 } else { 1183 ierr = PetscHMapIGet(htWithArtificial, globalDof + l, &localDof);CHKERRQ(ierr); 1184 if (localDof == -1) { 1185 localDof = localIndexWithArtificial++; 1186 ierr = PetscHMapISet(htWithArtificial, globalDof + l, localDof);CHKERRQ(ierr); 1187 } 1188 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 1189 /* And store.*/ 1190 dofsArrayWithArtificial[globalIndex] = localDof; 1191 } 1192 } 1193 1194 if(isNonlinear) { 1195 /* Build the dofmap for the function space with _all_ dofs, 1196 including those in any kind of boundary condition */ 1197 ierr = PetscHMapIGet(htWithAll, globalDof + l, &localDof);CHKERRQ(ierr); 1198 if (localDof == -1) { 1199 localDof = localIndexWithAll++; 1200 ierr = PetscHMapISet(htWithAll, globalDof + l, localDof);CHKERRQ(ierr); 1201 } 1202 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 1203 /* And store.*/ 1204 dofsArrayWithAll[globalIndex] = localDof; 1205 } 1206 globalIndex++; 1207 } 1208 } 1209 } 1210 } 1211 /*How many local dofs in this patch? */ 1212 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1213 ierr = PetscHMapIGetSize(htWithArtificial, &dof);CHKERRQ(ierr); 1214 ierr = PetscSectionSetDof(gtolCountsWithArtificial, v, dof);CHKERRQ(ierr); 1215 } 1216 if (isNonlinear) { 1217 ierr = PetscHMapIGetSize(htWithAll, &dof);CHKERRQ(ierr); 1218 ierr = PetscSectionSetDof(gtolCountsWithAll, v, dof);CHKERRQ(ierr); 1219 } 1220 ierr = PetscHMapIGetSize(ht, &dof);CHKERRQ(ierr); 1221 ierr = PetscSectionSetDof(gtolCounts, v, dof);CHKERRQ(ierr); 1222 } 1223 if (globalIndex != numDofs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%d) doesn't match found number (%d)", numDofs, globalIndex); 1224 ierr = PetscSectionSetUp(gtolCounts);CHKERRQ(ierr); 1225 ierr = PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);CHKERRQ(ierr); 1226 ierr = PetscMalloc1(numGlobalDofs, &globalDofsArray);CHKERRQ(ierr); 1227 1228 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1229 ierr = PetscSectionSetUp(gtolCountsWithArtificial);CHKERRQ(ierr); 1230 ierr = PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial);CHKERRQ(ierr); 1231 ierr = PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial);CHKERRQ(ierr); 1232 } 1233 if (isNonlinear) { 1234 ierr = PetscSectionSetUp(gtolCountsWithAll);CHKERRQ(ierr); 1235 ierr = PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll);CHKERRQ(ierr); 1236 ierr = PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll);CHKERRQ(ierr); 1237 } 1238 /* Now populate the global to local map. This could be merged into the above loop if we were willing to deal with reallocs. */ 1239 for (v = vStart; v < vEnd; ++v) { 1240 PetscHashIter hi; 1241 PetscInt dof, off, Np, ooff, i, j, k, l; 1242 1243 ierr = PetscHMapIClear(ht);CHKERRQ(ierr); 1244 ierr = PetscHMapIClear(htWithArtificial);CHKERRQ(ierr); 1245 ierr = PetscHMapIClear(htWithAll);CHKERRQ(ierr); 1246 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 1247 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 1248 ierr = PetscSectionGetDof(pointCounts, v, &Np);CHKERRQ(ierr); 1249 ierr = PetscSectionGetOffset(pointCounts, v, &ooff);CHKERRQ(ierr); 1250 if (dof <= 0) continue; 1251 1252 for (k = 0; k < patch->nsubspaces; ++k) { 1253 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1254 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1255 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1256 PetscInt bs = patch->bs[k]; 1257 PetscInt goff; 1258 1259 for (i = off; i < off + dof; ++i) { 1260 /* Reconstruct mapping of global-to-local on this patch. */ 1261 const PetscInt c = cellsArray[i]; 1262 PetscInt cell = c; 1263 1264 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 1265 for (j = 0; j < nodesPerCell; ++j) { 1266 for (l = 0; l < bs; ++l) { 1267 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1268 const PetscInt localDof = dofsArray[key]; 1269 if (localDof >= 0) {ierr = PetscHMapISet(ht, globalDof, localDof);CHKERRQ(ierr);} 1270 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1271 const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key]; 1272 if (localDofWithArtificial >= 0) { 1273 ierr = PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial);CHKERRQ(ierr); 1274 } 1275 } 1276 if (isNonlinear) { 1277 const PetscInt localDofWithAll = dofsArrayWithAll[key]; 1278 if (localDofWithAll >= 0) { 1279 ierr = PetscHMapISet(htWithAll, globalDof, localDofWithAll);CHKERRQ(ierr); 1280 } 1281 } 1282 key++; 1283 } 1284 } 1285 } 1286 1287 /* Shove it in the output data structure. */ 1288 ierr = PetscSectionGetOffset(gtolCounts, v, &goff);CHKERRQ(ierr); 1289 PetscHashIterBegin(ht, hi); 1290 while (!PetscHashIterAtEnd(ht, hi)) { 1291 PetscInt globalDof, localDof; 1292 1293 PetscHashIterGetKey(ht, hi, globalDof); 1294 PetscHashIterGetVal(ht, hi, localDof); 1295 if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof; 1296 PetscHashIterNext(ht, hi); 1297 } 1298 1299 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1300 ierr = PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff);CHKERRQ(ierr); 1301 PetscHashIterBegin(htWithArtificial, hi); 1302 while (!PetscHashIterAtEnd(htWithArtificial, hi)) { 1303 PetscInt globalDof, localDof; 1304 PetscHashIterGetKey(htWithArtificial, hi, globalDof); 1305 PetscHashIterGetVal(htWithArtificial, hi, localDof); 1306 if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof; 1307 PetscHashIterNext(htWithArtificial, hi); 1308 } 1309 } 1310 if (isNonlinear) { 1311 ierr = PetscSectionGetOffset(gtolCountsWithAll, v, &goff);CHKERRQ(ierr); 1312 PetscHashIterBegin(htWithAll, hi); 1313 while (!PetscHashIterAtEnd(htWithAll, hi)) { 1314 PetscInt globalDof, localDof; 1315 PetscHashIterGetKey(htWithAll, hi, globalDof); 1316 PetscHashIterGetVal(htWithAll, hi, localDof); 1317 if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof; 1318 PetscHashIterNext(htWithAll, hi); 1319 } 1320 } 1321 1322 for (p = 0; p < Np; ++p) { 1323 const PetscInt point = pointsArray[ooff + p]; 1324 PetscInt globalDof, localDof; 1325 1326 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);CHKERRQ(ierr); 1327 ierr = PetscHMapIGet(ht, globalDof, &localDof);CHKERRQ(ierr); 1328 offsArray[(ooff + p)*Nf + k] = localDof; 1329 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1330 ierr = PetscHMapIGet(htWithArtificial, globalDof, &localDof);CHKERRQ(ierr); 1331 offsArrayWithArtificial[(ooff + p)*Nf + k] = localDof; 1332 } 1333 if (isNonlinear) { 1334 ierr = PetscHMapIGet(htWithAll, globalDof, &localDof);CHKERRQ(ierr); 1335 offsArrayWithAll[(ooff + p)*Nf + k] = localDof; 1336 } 1337 } 1338 } 1339 1340 ierr = PetscHSetIDestroy(&globalBcs);CHKERRQ(ierr); 1341 ierr = PetscHSetIDestroy(&ownedpts);CHKERRQ(ierr); 1342 ierr = PetscHSetIDestroy(&seenpts);CHKERRQ(ierr); 1343 ierr = PetscHSetIDestroy(&owneddofs);CHKERRQ(ierr); 1344 ierr = PetscHSetIDestroy(&seendofs);CHKERRQ(ierr); 1345 ierr = PetscHSetIDestroy(&artificialbcs);CHKERRQ(ierr); 1346 1347 /* At this point, we have a hash table ht built that maps globalDof -> localDof. 1348 We need to create the dof table laid out cellwise first, then by subspace, 1349 as the assembler assembles cell-wise and we need to stuff the different 1350 contributions of the different function spaces to the right places. So we loop 1351 over cells, then over subspaces. */ 1352 if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */ 1353 for (i = off; i < off + dof; ++i) { 1354 const PetscInt c = cellsArray[i]; 1355 PetscInt cell = c; 1356 1357 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 1358 for (k = 0; k < patch->nsubspaces; ++k) { 1359 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1360 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1361 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1362 PetscInt bs = patch->bs[k]; 1363 1364 for (j = 0; j < nodesPerCell; ++j) { 1365 for (l = 0; l < bs; ++l) { 1366 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1367 PetscInt localDof; 1368 1369 ierr = PetscHMapIGet(ht, globalDof, &localDof);CHKERRQ(ierr); 1370 /* If it's not in the hash table, i.e. is a BC dof, 1371 then the PetscHSetIMap above gives -1, which matches 1372 exactly the convention for PETSc's matrix assembly to 1373 ignore the dof. So we don't need to do anything here */ 1374 asmArray[asmKey] = localDof; 1375 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1376 ierr = PetscHMapIGet(htWithArtificial, globalDof, &localDof);CHKERRQ(ierr); 1377 asmArrayWithArtificial[asmKey] = localDof; 1378 } 1379 if (isNonlinear) { 1380 ierr = PetscHMapIGet(htWithAll, globalDof, &localDof);CHKERRQ(ierr); 1381 asmArrayWithAll[asmKey] = localDof; 1382 } 1383 asmKey++; 1384 } 1385 } 1386 } 1387 } 1388 } 1389 } 1390 if (1 == patch->nsubspaces) { 1391 ierr = PetscMemcpy(asmArray, dofsArray, numDofs * sizeof(PetscInt));CHKERRQ(ierr); 1392 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1393 ierr = PetscMemcpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs * sizeof(PetscInt));CHKERRQ(ierr); 1394 } 1395 if (isNonlinear) { 1396 ierr = PetscMemcpy(asmArrayWithAll, dofsArrayWithAll, numDofs * sizeof(PetscInt));CHKERRQ(ierr); 1397 } 1398 } 1399 1400 ierr = PetscHMapIDestroy(&ht);CHKERRQ(ierr); 1401 ierr = PetscHMapIDestroy(&htWithArtificial);CHKERRQ(ierr); 1402 ierr = PetscHMapIDestroy(&htWithAll);CHKERRQ(ierr); 1403 ierr = ISRestoreIndices(cells, &cellsArray);CHKERRQ(ierr); 1404 ierr = ISRestoreIndices(points, &pointsArray);CHKERRQ(ierr); 1405 ierr = PetscFree(dofsArray);CHKERRQ(ierr); 1406 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1407 ierr = PetscFree(dofsArrayWithArtificial);CHKERRQ(ierr); 1408 } 1409 if (isNonlinear) { 1410 ierr = PetscFree(dofsArrayWithAll);CHKERRQ(ierr); 1411 } 1412 /* Create placeholder section for map from points to patch dofs */ 1413 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);CHKERRQ(ierr); 1414 ierr = PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);CHKERRQ(ierr); 1415 if (patch->combined) { 1416 PetscInt numFields; 1417 ierr = PetscSectionGetNumFields(patch->dofSection[0], &numFields);CHKERRQ(ierr); 1418 if (numFields != patch->nsubspaces) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Mismatch between number of section fields %D and number of subspaces %D", numFields, patch->nsubspaces); 1419 ierr = PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);CHKERRQ(ierr); 1420 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1421 for (p = pStart; p < pEnd; ++p) { 1422 PetscInt dof, fdof, f; 1423 1424 ierr = PetscSectionGetDof(patch->dofSection[0], p, &dof);CHKERRQ(ierr); 1425 ierr = PetscSectionSetDof(patch->patchSection, p, dof);CHKERRQ(ierr); 1426 for (f = 0; f < patch->nsubspaces; ++f) { 1427 ierr = PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof);CHKERRQ(ierr); 1428 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1429 } 1430 } 1431 } else { 1432 PetscInt pStartf, pEndf, f; 1433 pStart = PETSC_MAX_INT; 1434 pEnd = PETSC_MIN_INT; 1435 for (f = 0; f < patch->nsubspaces; ++f) { 1436 ierr = PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);CHKERRQ(ierr); 1437 pStart = PetscMin(pStart, pStartf); 1438 pEnd = PetscMax(pEnd, pEndf); 1439 } 1440 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1441 for (f = 0; f < patch->nsubspaces; ++f) { 1442 ierr = PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);CHKERRQ(ierr); 1443 for (p = pStartf; p < pEndf; ++p) { 1444 PetscInt fdof; 1445 ierr = PetscSectionGetDof(patch->dofSection[f], p, &fdof);CHKERRQ(ierr); 1446 ierr = PetscSectionAddDof(patch->patchSection, p, fdof);CHKERRQ(ierr); 1447 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1448 } 1449 } 1450 } 1451 ierr = PetscSectionSetUp(patch->patchSection);CHKERRQ(ierr); 1452 ierr = PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE);CHKERRQ(ierr); 1453 /* Replace cell indices with firedrake-numbered ones. */ 1454 ierr = ISGeneralSetIndices(cells, numCells, (const PetscInt *) newCellsArray, PETSC_OWN_POINTER);CHKERRQ(ierr); 1455 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol);CHKERRQ(ierr); 1456 ierr = PetscObjectSetName((PetscObject) patch->gtol, "Global Indices");CHKERRQ(ierr); 1457 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_g2l_view", patch->classname);CHKERRQ(ierr); 1458 ierr = PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject) pc, option);CHKERRQ(ierr); 1459 ierr = ISViewFromOptions(patch->gtol, (PetscObject) pc, option);CHKERRQ(ierr); 1460 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs);CHKERRQ(ierr); 1461 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArray, PETSC_OWN_POINTER, &patch->offs);CHKERRQ(ierr); 1462 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1463 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithArtificial, globalDofsArrayWithArtificial, PETSC_OWN_POINTER, &patch->gtolWithArtificial);CHKERRQ(ierr); 1464 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithArtificial, PETSC_OWN_POINTER, &patch->dofsWithArtificial);CHKERRQ(ierr); 1465 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithArtificial, PETSC_OWN_POINTER, &patch->offsWithArtificial);CHKERRQ(ierr); 1466 } 1467 if (isNonlinear) { 1468 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithAll, globalDofsArrayWithAll, PETSC_OWN_POINTER, &patch->gtolWithAll);CHKERRQ(ierr); 1469 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithAll, PETSC_OWN_POINTER, &patch->dofsWithAll);CHKERRQ(ierr); 1470 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithAll, PETSC_OWN_POINTER, &patch->offsWithAll);CHKERRQ(ierr); 1471 } 1472 PetscFunctionReturn(0); 1473 } 1474 1475 static PetscErrorCode PCPatchZeroFillMatrix_Private(Mat mat, const PetscInt ncell, const PetscInt ndof, const PetscInt *dof) 1476 { 1477 PetscScalar *values = NULL; 1478 PetscInt rows, c, i; 1479 PetscErrorCode ierr; 1480 1481 PetscFunctionBegin; 1482 ierr = PetscCalloc1(ndof*ndof, &values);CHKERRQ(ierr); 1483 for (c = 0; c < ncell; ++c) { 1484 const PetscInt *idx = &dof[ndof*c]; 1485 ierr = MatSetValues(mat, ndof, idx, ndof, idx, values, INSERT_VALUES);CHKERRQ(ierr); 1486 } 1487 ierr = MatGetLocalSize(mat, &rows, NULL);CHKERRQ(ierr); 1488 for (i = 0; i < rows; ++i) { 1489 ierr = MatSetValues(mat, 1, &i, 1, &i, values, INSERT_VALUES);CHKERRQ(ierr); 1490 } 1491 ierr = MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1492 ierr = MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1493 ierr = PetscFree(values);CHKERRQ(ierr); 1494 PetscFunctionReturn(0); 1495 } 1496 1497 static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial) 1498 { 1499 PC_PATCH *patch = (PC_PATCH *) pc->data; 1500 Vec x, y; 1501 PetscBool flg; 1502 PetscInt csize, rsize; 1503 const char *prefix = NULL; 1504 PetscErrorCode ierr; 1505 1506 PetscFunctionBegin; 1507 if(withArtificial) { 1508 /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */ 1509 x = patch->patchRHSWithArtificial[point]; 1510 y = patch->patchRHSWithArtificial[point]; 1511 } else { 1512 x = patch->patchRHS[point]; 1513 y = patch->patchUpdate[point]; 1514 } 1515 1516 ierr = VecGetSize(x, &csize);CHKERRQ(ierr); 1517 ierr = VecGetSize(y, &rsize);CHKERRQ(ierr); 1518 ierr = MatCreate(PETSC_COMM_SELF, mat);CHKERRQ(ierr); 1519 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1520 ierr = MatSetOptionsPrefix(*mat, prefix);CHKERRQ(ierr); 1521 ierr = MatAppendOptionsPrefix(*mat, "pc_patch_sub_");CHKERRQ(ierr); 1522 if (patch->sub_mat_type) {ierr = MatSetType(*mat, patch->sub_mat_type);CHKERRQ(ierr);} 1523 else if (!patch->sub_mat_type) {ierr = MatSetType(*mat, MATDENSE);CHKERRQ(ierr);} 1524 ierr = MatSetSizes(*mat, rsize, csize, rsize, csize);CHKERRQ(ierr); 1525 ierr = PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);CHKERRQ(ierr); 1526 if (!flg) {ierr = PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);CHKERRQ(ierr);} 1527 /* Sparse patch matrices */ 1528 if (!flg) { 1529 PetscBT bt; 1530 PetscInt *dnnz = NULL; 1531 const PetscInt *dofsArray = NULL; 1532 PetscInt pStart, pEnd, ncell, offset, c, i, j; 1533 1534 if(withArtificial) { 1535 ierr = ISGetIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 1536 } else { 1537 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1538 } 1539 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1540 point += pStart; 1541 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1542 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1543 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1544 ierr = PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1545 /* A PetscBT uses N^2 bits to store the sparsity pattern on a 1546 * patch. This is probably OK if the patches are not too big, 1547 * but uses too much memory. We therefore switch based on rsize. */ 1548 if (rsize < 3000) { /* FIXME: I picked this switch value out of my hat */ 1549 ierr = PetscCalloc1(rsize, &dnnz);CHKERRQ(ierr); 1550 ierr = PetscBTCreate(rsize*rsize, &bt);CHKERRQ(ierr); 1551 for (c = 0; c < ncell; ++c) { 1552 const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell; 1553 for (i = 0; i < patch->totalDofsPerCell; ++i) { 1554 const PetscInt row = idx[i]; 1555 if (row < 0) continue; 1556 for (j = 0; j < patch->totalDofsPerCell; ++j) { 1557 const PetscInt col = idx[j]; 1558 const PetscInt key = row*rsize + col; 1559 if (col < 0) continue; 1560 if (!PetscBTLookupSet(bt, key)) ++dnnz[row]; 1561 } 1562 } 1563 } 1564 ierr = PetscBTDestroy(&bt);CHKERRQ(ierr); 1565 ierr = MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);CHKERRQ(ierr); 1566 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1567 ierr = PCPatchZeroFillMatrix_Private(*mat, ncell, patch->totalDofsPerCell, &dofsArray[offset*patch->totalDofsPerCell]);CHKERRQ(ierr); 1568 } else { /* rsize too big, use MATPREALLOCATOR */ 1569 Mat preallocator; 1570 PetscScalar* vals; 1571 1572 ierr = PetscCalloc1(patch->totalDofsPerCell*patch->totalDofsPerCell, &vals);CHKERRQ(ierr); 1573 ierr = MatCreate(PETSC_COMM_SELF, &preallocator);CHKERRQ(ierr); 1574 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1575 ierr = MatSetSizes(preallocator, rsize, rsize, rsize, rsize);CHKERRQ(ierr); 1576 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1577 for (c = 0; c < ncell; ++c) { 1578 const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell; 1579 ierr = MatSetValues(preallocator, patch->totalDofsPerCell, idx, patch->totalDofsPerCell, idx, vals, INSERT_VALUES);CHKERRQ(ierr); 1580 } 1581 ierr = PetscFree(vals);CHKERRQ(ierr); 1582 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1583 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1584 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, *mat);CHKERRQ(ierr); 1585 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1586 } 1587 ierr = PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1588 if(withArtificial) { 1589 ierr = ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 1590 } else { 1591 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1592 } 1593 } 1594 ierr = MatSetUp(*mat);CHKERRQ(ierr); 1595 PetscFunctionReturn(0); 1596 } 1597 1598 static PetscErrorCode PCPatchComputeFunction_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Vec F, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx) 1599 { 1600 PC_PATCH *patch = (PC_PATCH *) pc->data; 1601 DM dm; 1602 PetscSection s; 1603 const PetscInt *parray, *oarray; 1604 PetscInt Nf = patch->nsubspaces, Np, poff, p, f; 1605 PetscErrorCode ierr; 1606 1607 PetscFunctionBegin; 1608 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1609 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1610 /* Set offset into patch */ 1611 ierr = PetscSectionGetDof(patch->pointCounts, patchNum, &Np);CHKERRQ(ierr); 1612 ierr = PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);CHKERRQ(ierr); 1613 ierr = ISGetIndices(patch->points, &parray);CHKERRQ(ierr); 1614 ierr = ISGetIndices(patch->offs, &oarray);CHKERRQ(ierr); 1615 for (f = 0; f < Nf; ++f) { 1616 for (p = 0; p < Np; ++p) { 1617 const PetscInt point = parray[poff+p]; 1618 PetscInt dof; 1619 1620 ierr = PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);CHKERRQ(ierr); 1621 ierr = PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr); 1622 if (patch->nsubspaces == 1) {ierr = PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr);} 1623 else {ierr = PetscSectionSetOffset(patch->patchSection, point, -1);CHKERRQ(ierr);} 1624 } 1625 } 1626 ierr = ISRestoreIndices(patch->points, &parray);CHKERRQ(ierr); 1627 ierr = ISRestoreIndices(patch->offs, &oarray);CHKERRQ(ierr); 1628 if (patch->viewSection) {ierr = ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);CHKERRQ(ierr);} 1629 ierr = DMPlexComputeResidual_Patch_Internal(pc->dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx);CHKERRQ(ierr); 1630 PetscFunctionReturn(0); 1631 } 1632 1633 PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point) 1634 { 1635 PC_PATCH *patch = (PC_PATCH *) pc->data; 1636 const PetscInt *dofsArray; 1637 const PetscInt *dofsArrayWithAll; 1638 const PetscInt *cellsArray; 1639 PetscInt ncell, offset, pStart, pEnd; 1640 PetscErrorCode ierr; 1641 1642 PetscFunctionBegin; 1643 ierr = PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1644 if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n"); 1645 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1646 ierr = ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 1647 ierr = ISGetIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1648 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1649 1650 point += pStart; 1651 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1652 1653 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1654 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1655 if (ncell <= 0) { 1656 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1657 PetscFunctionReturn(0); 1658 } 1659 PetscStackPush("PCPatch user callback"); 1660 /* Cannot reuse the same IS because the geometry info is being cached in it */ 1661 ierr = ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);CHKERRQ(ierr); 1662 ierr = patch->usercomputef(pc, point, x, F, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, 1663 dofsArrayWithAll + offset*patch->totalDofsPerCell, 1664 patch->usercomputefctx);CHKERRQ(ierr); 1665 PetscStackPop; 1666 ierr = ISDestroy(&patch->cellIS);CHKERRQ(ierr); 1667 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1668 ierr = ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 1669 ierr = ISRestoreIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1670 if (patch->viewMatrix) { 1671 char name[PETSC_MAX_PATH_LEN]; 1672 1673 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch vector for Point %D", point);CHKERRQ(ierr); 1674 ierr = PetscObjectSetName((PetscObject) F, name);CHKERRQ(ierr); 1675 ierr = ObjectView((PetscObject) F, patch->viewerMatrix, patch->formatMatrix);CHKERRQ(ierr); 1676 } 1677 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1678 PetscFunctionReturn(0); 1679 } 1680 1681 static PetscErrorCode PCPatchComputeOperator_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Mat J, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx) 1682 { 1683 PC_PATCH *patch = (PC_PATCH *) pc->data; 1684 DM dm; 1685 PetscSection s; 1686 const PetscInt *parray, *oarray; 1687 PetscInt Nf = patch->nsubspaces, Np, poff, p, f; 1688 PetscErrorCode ierr; 1689 1690 PetscFunctionBegin; 1691 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1692 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1693 /* Set offset into patch */ 1694 ierr = PetscSectionGetDof(patch->pointCounts, patchNum, &Np);CHKERRQ(ierr); 1695 ierr = PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);CHKERRQ(ierr); 1696 ierr = ISGetIndices(patch->points, &parray);CHKERRQ(ierr); 1697 ierr = ISGetIndices(patch->offs, &oarray);CHKERRQ(ierr); 1698 for (f = 0; f < Nf; ++f) { 1699 for (p = 0; p < Np; ++p) { 1700 const PetscInt point = parray[poff+p]; 1701 PetscInt dof; 1702 1703 ierr = PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);CHKERRQ(ierr); 1704 ierr = PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr); 1705 if (patch->nsubspaces == 1) {ierr = PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr);} 1706 else {ierr = PetscSectionSetOffset(patch->patchSection, point, -1);CHKERRQ(ierr);} 1707 } 1708 } 1709 ierr = ISRestoreIndices(patch->points, &parray);CHKERRQ(ierr); 1710 ierr = ISRestoreIndices(patch->offs, &oarray);CHKERRQ(ierr); 1711 if (patch->viewSection) {ierr = ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);CHKERRQ(ierr);} 1712 /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */ 1713 ierr = DMPlexComputeJacobian_Patch_Internal(pc->dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx);CHKERRQ(ierr); 1714 PetscFunctionReturn(0); 1715 } 1716 1717 PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial) 1718 { 1719 PC_PATCH *patch = (PC_PATCH *) pc->data; 1720 const PetscInt *dofsArray; 1721 const PetscInt *dofsArrayWithArtificial = NULL; 1722 const PetscInt *dofsArrayWithAll = NULL; 1723 const PetscInt *cellsArray; 1724 PetscInt ncell, offset, pStart, pEnd; 1725 PetscBool isNonlinear; 1726 PetscErrorCode ierr; 1727 1728 PetscFunctionBegin; 1729 ierr = PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1730 isNonlinear = patch->isNonlinear; 1731 if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n"); 1732 if(withArtificial) { 1733 ierr = ISGetIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 1734 } else { 1735 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1736 } 1737 if (isNonlinear) { 1738 ierr = ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 1739 } 1740 ierr = ISGetIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1741 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1742 1743 point += pStart; 1744 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1745 1746 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1747 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1748 if (ncell <= 0) { 1749 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1750 PetscFunctionReturn(0); 1751 } 1752 PetscStackPush("PCPatch user callback"); 1753 /* Cannot reuse the same IS because the geometry info is being cached in it */ 1754 ierr = ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);CHKERRQ(ierr); 1755 ierr = patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, dofsArrayWithAll ? dofsArrayWithAll + offset*patch->totalDofsPerCell : NULL, patch->usercomputeopctx);CHKERRQ(ierr); 1756 PetscStackPop; 1757 ierr = ISDestroy(&patch->cellIS);CHKERRQ(ierr); 1758 if(withArtificial) { 1759 ierr = ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 1760 } else { 1761 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1762 } 1763 if (isNonlinear) { 1764 ierr = ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 1765 } 1766 ierr = ISRestoreIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1767 if (patch->viewMatrix) { 1768 char name[PETSC_MAX_PATH_LEN]; 1769 1770 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);CHKERRQ(ierr); 1771 ierr = PetscObjectSetName((PetscObject) mat, name);CHKERRQ(ierr); 1772 ierr = ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);CHKERRQ(ierr); 1773 } 1774 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1775 PetscFunctionReturn(0); 1776 } 1777 1778 PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype) 1779 { 1780 PC_PATCH *patch = (PC_PATCH *) pc->data; 1781 const PetscScalar *xArray = NULL; 1782 PetscScalar *yArray = NULL; 1783 const PetscInt *gtolArray = NULL; 1784 PetscInt dof, offset, lidx; 1785 PetscErrorCode ierr; 1786 1787 PetscFunctionBeginHot; 1788 ierr = PetscLogEventBegin(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 1789 ierr = VecGetArrayRead(x, &xArray);CHKERRQ(ierr); 1790 ierr = VecGetArray(y, &yArray);CHKERRQ(ierr); 1791 if (scattertype == SCATTER_WITHARTIFICIAL) { 1792 ierr = PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);CHKERRQ(ierr); 1793 ierr = PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset);CHKERRQ(ierr); 1794 ierr = ISGetIndices(patch->gtolWithArtificial, >olArray);CHKERRQ(ierr); 1795 } else if (scattertype == SCATTER_WITHALL) { 1796 ierr = PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof);CHKERRQ(ierr); 1797 ierr = PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset);CHKERRQ(ierr); 1798 ierr = ISGetIndices(patch->gtolWithAll, >olArray);CHKERRQ(ierr); 1799 } else { 1800 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 1801 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 1802 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1803 } 1804 if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n"); 1805 if (mode == ADD_VALUES && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n"); 1806 for (lidx = 0; lidx < dof; ++lidx) { 1807 const PetscInt gidx = gtolArray[offset+lidx]; 1808 1809 if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */ 1810 else yArray[gidx] += xArray[lidx]; /* Reverse */ 1811 } 1812 if (scattertype == SCATTER_WITHARTIFICIAL) { 1813 ierr = ISRestoreIndices(patch->gtolWithArtificial, >olArray);CHKERRQ(ierr); 1814 } else if (scattertype == SCATTER_WITHALL) { 1815 ierr = ISRestoreIndices(patch->gtolWithAll, >olArray);CHKERRQ(ierr); 1816 } else { 1817 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1818 } 1819 ierr = VecRestoreArrayRead(x, &xArray);CHKERRQ(ierr); 1820 ierr = VecRestoreArray(y, &yArray);CHKERRQ(ierr); 1821 ierr = PetscLogEventEnd(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 1822 PetscFunctionReturn(0); 1823 } 1824 1825 static PetscErrorCode PCSetUp_PATCH_Linear(PC pc) 1826 { 1827 PC_PATCH *patch = (PC_PATCH *) pc->data; 1828 const char *prefix; 1829 PetscInt i; 1830 PetscErrorCode ierr; 1831 1832 PetscFunctionBegin; 1833 if (!pc->setupcalled) { 1834 ierr = PetscMalloc1(patch->npatch, &patch->solver);CHKERRQ(ierr); 1835 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1836 for (i = 0; i < patch->npatch; ++i) { 1837 KSP ksp; 1838 PC subpc; 1839 1840 ierr = KSPCreate(PETSC_COMM_SELF, &ksp);CHKERRQ(ierr); 1841 ierr = KSPSetOptionsPrefix(ksp, prefix);CHKERRQ(ierr); 1842 ierr = KSPAppendOptionsPrefix(ksp, "sub_");CHKERRQ(ierr); 1843 ierr = PetscObjectIncrementTabLevel((PetscObject) ksp, (PetscObject) pc, 1);CHKERRQ(ierr); 1844 ierr = KSPGetPC(ksp, &subpc);CHKERRQ(ierr); 1845 ierr = PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);CHKERRQ(ierr); 1846 ierr = PetscLogObjectParent((PetscObject) pc, (PetscObject) ksp);CHKERRQ(ierr); 1847 patch->solver[i] = (PetscObject) ksp; 1848 } 1849 } 1850 if (patch->save_operators) { 1851 for (i = 0; i < patch->npatch; ++i) { 1852 ierr = MatZeroEntries(patch->mat[i]);CHKERRQ(ierr); 1853 ierr = PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE);CHKERRQ(ierr); 1854 ierr = KSPSetOperators((KSP) patch->solver[i], patch->mat[i], patch->mat[i]);CHKERRQ(ierr); 1855 } 1856 } 1857 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1858 for (i = 0; i < patch->npatch; ++i) { 1859 /* Instead of padding patch->patchUpdate with zeros to get */ 1860 /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */ 1861 /* just get rid of the columns that correspond to the dofs with */ 1862 /* artificial bcs. That's of course fairly inefficient, hopefully we */ 1863 /* can just assemble the rectangular matrix in the first place. */ 1864 Mat matSquare; 1865 IS rowis; 1866 PetscInt dof; 1867 1868 ierr = MatGetSize(patch->mat[i], &dof, NULL);CHKERRQ(ierr); 1869 if (dof == 0) { 1870 patch->matWithArtificial[i] = NULL; 1871 continue; 1872 } 1873 1874 ierr = PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);CHKERRQ(ierr); 1875 ierr = MatZeroEntries(matSquare);CHKERRQ(ierr); 1876 ierr = PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);CHKERRQ(ierr); 1877 1878 ierr = MatGetSize(matSquare, &dof, NULL);CHKERRQ(ierr); 1879 ierr = ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis); CHKERRQ(ierr); 1880 if(pc->setupcalled) { 1881 ierr = MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]); CHKERRQ(ierr); 1882 } else { 1883 ierr = MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]); CHKERRQ(ierr); 1884 } 1885 ierr = ISDestroy(&rowis); CHKERRQ(ierr); 1886 ierr = MatDestroy(&matSquare);CHKERRQ(ierr); 1887 } 1888 } 1889 PetscFunctionReturn(0); 1890 } 1891 1892 static PetscErrorCode PCSetUp_PATCH(PC pc) 1893 { 1894 PC_PATCH *patch = (PC_PATCH *) pc->data; 1895 PetscInt i; 1896 PetscBool isNonlinear; 1897 PetscErrorCode ierr; 1898 1899 PetscFunctionBegin; 1900 if (!pc->setupcalled) { 1901 PetscInt pStart, pEnd, p; 1902 PetscInt localSize; 1903 1904 ierr = PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 1905 1906 isNonlinear = patch->isNonlinear; 1907 if (!patch->nsubspaces) { 1908 DM dm; 1909 PetscDS prob; 1910 PetscSection s; 1911 PetscInt cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, totNb = 0, **cellDofs; 1912 1913 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1914 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()"); 1915 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1916 ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr); 1917 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr); 1918 for (p = pStart; p < pEnd; ++p) { 1919 PetscInt cdof; 1920 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 1921 numGlobalBcs += cdof; 1922 } 1923 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1924 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1925 ierr = PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);CHKERRQ(ierr); 1926 for (f = 0; f < Nf; ++f) { 1927 PetscFE fe; 1928 PetscDualSpace sp; 1929 PetscInt cdoff = 0; 1930 1931 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1932 /* ierr = PetscFEGetNumComponents(fe, &Nc[f]);CHKERRQ(ierr); */ 1933 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 1934 ierr = PetscDualSpaceGetDimension(sp, &Nb[f]);CHKERRQ(ierr); 1935 totNb += Nb[f]; 1936 1937 ierr = PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);CHKERRQ(ierr); 1938 for (c = cStart; c < cEnd; ++c) { 1939 PetscInt *closure = NULL; 1940 PetscInt clSize = 0, cl; 1941 1942 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1943 for (cl = 0; cl < clSize*2; cl += 2) { 1944 const PetscInt p = closure[cl]; 1945 PetscInt fdof, d, foff; 1946 1947 ierr = PetscSectionGetFieldDof(s, p, f, &fdof);CHKERRQ(ierr); 1948 ierr = PetscSectionGetFieldOffset(s, p, f, &foff);CHKERRQ(ierr); 1949 for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d; 1950 } 1951 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1952 } 1953 if (cdoff != (cEnd-cStart)*Nb[f]) SETERRQ4(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_SIZ, "Total number of cellDofs %D for field %D should be Nc (%D) * cellDof (%D)", cdoff, f, cEnd-cStart, Nb[f]); 1954 } 1955 numGlobalBcs = 0; 1956 for (p = pStart; p < pEnd; ++p) { 1957 const PetscInt *ind; 1958 PetscInt off, cdof, d; 1959 1960 ierr = PetscSectionGetOffset(s, p, &off);CHKERRQ(ierr); 1961 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 1962 ierr = PetscSectionGetConstraintIndices(s, p, &ind);CHKERRQ(ierr); 1963 for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d]; 1964 } 1965 1966 ierr = PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);CHKERRQ(ierr); 1967 for (f = 0; f < Nf; ++f) { 1968 ierr = PetscFree(cellDofs[f]);CHKERRQ(ierr); 1969 } 1970 ierr = PetscFree3(Nb, cellDofs, globalBcs);CHKERRQ(ierr); 1971 ierr = PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL);CHKERRQ(ierr); 1972 ierr = PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);CHKERRQ(ierr); 1973 } 1974 1975 localSize = patch->subspaceOffsets[patch->nsubspaces]; 1976 ierr = VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS);CHKERRQ(ierr); 1977 ierr = VecSetUp(patch->localRHS);CHKERRQ(ierr); 1978 ierr = VecDuplicate(patch->localRHS, &patch->localUpdate);CHKERRQ(ierr); 1979 ierr = PCPatchCreateCellPatches(pc);CHKERRQ(ierr); 1980 ierr = PCPatchCreateCellPatchDiscretisationInfo(pc);CHKERRQ(ierr); 1981 1982 /* OK, now build the work vectors */ 1983 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd);CHKERRQ(ierr); 1984 ierr = PetscMalloc1(patch->npatch, &patch->patchRHS);CHKERRQ(ierr); 1985 ierr = PetscMalloc1(patch->npatch, &patch->patchUpdate);CHKERRQ(ierr); 1986 1987 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1988 ierr = PetscMalloc1(patch->npatch, &patch->patchRHSWithArtificial);CHKERRQ(ierr); 1989 ierr = PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial);CHKERRQ(ierr); 1990 } 1991 if (isNonlinear) { 1992 ierr = PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll);CHKERRQ(ierr); 1993 } 1994 for (p = pStart; p < pEnd; ++p) { 1995 PetscInt dof; 1996 1997 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 1998 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchRHS[p-pStart]);CHKERRQ(ierr); 1999 ierr = VecSetUp(patch->patchRHS[p-pStart]);CHKERRQ(ierr); 2000 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchUpdate[p-pStart]);CHKERRQ(ierr); 2001 ierr = VecSetUp(patch->patchUpdate[p-pStart]);CHKERRQ(ierr); 2002 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 2003 const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL; 2004 PetscInt numPatchDofs, offset; 2005 PetscInt numPatchDofsWithArtificial, offsetWithArtificial; 2006 PetscInt dofWithoutArtificialCounter = 0; 2007 PetscInt *patchWithoutArtificialToWithArtificialArray; 2008 2009 ierr = PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);CHKERRQ(ierr); 2010 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchRHSWithArtificial[p-pStart]);CHKERRQ(ierr); 2011 ierr = VecSetUp(patch->patchRHSWithArtificial[p-pStart]);CHKERRQ(ierr); 2012 2013 /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */ 2014 /* the index in the patch with all dofs */ 2015 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 2016 2017 ierr = PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);CHKERRQ(ierr); 2018 if (numPatchDofs == 0) { 2019 patch->dofMappingWithoutToWithArtificial[p-pStart] = NULL; 2020 continue; 2021 } 2022 2023 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 2024 ierr = ISGetIndices(patch->gtolWithArtificial, >olArrayWithArtificial);CHKERRQ(ierr); 2025 ierr = PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial);CHKERRQ(ierr); 2026 ierr = PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial);CHKERRQ(ierr); 2027 2028 ierr = PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray);CHKERRQ(ierr); 2029 for (i=0; i<numPatchDofsWithArtificial; i++) { 2030 if (gtolArrayWithArtificial[i+offsetWithArtificial] == gtolArray[offset+dofWithoutArtificialCounter]) { 2031 patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i; 2032 dofWithoutArtificialCounter++; 2033 if (dofWithoutArtificialCounter == numPatchDofs) 2034 break; 2035 } 2036 } 2037 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p-pStart]);CHKERRQ(ierr); 2038 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 2039 ierr = ISRestoreIndices(patch->gtolWithArtificial, >olArrayWithArtificial);CHKERRQ(ierr); 2040 } 2041 if (isNonlinear) { 2042 const PetscInt *gtolArray, *gtolArrayWithAll = NULL; 2043 PetscInt numPatchDofs, offset; 2044 PetscInt numPatchDofsWithAll, offsetWithAll; 2045 PetscInt dofWithoutAllCounter = 0; 2046 PetscInt *patchWithoutAllToWithAllArray; 2047 2048 /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */ 2049 /* the index in the patch with all dofs */ 2050 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 2051 2052 ierr = PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);CHKERRQ(ierr); 2053 if (numPatchDofs == 0) { 2054 patch->dofMappingWithoutToWithAll[p-pStart] = NULL; 2055 continue; 2056 } 2057 2058 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 2059 ierr = ISGetIndices(patch->gtolWithAll, >olArrayWithAll);CHKERRQ(ierr); 2060 ierr = PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll);CHKERRQ(ierr); 2061 ierr = PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll);CHKERRQ(ierr); 2062 2063 ierr = PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);CHKERRQ(ierr); 2064 2065 for (i=0; i<numPatchDofsWithAll; i++) { 2066 if (gtolArrayWithAll[i+offsetWithAll] == gtolArray[offset+dofWithoutAllCounter]) { 2067 patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i; 2068 dofWithoutAllCounter++; 2069 if (dofWithoutAllCounter == numPatchDofs) 2070 break; 2071 } 2072 } 2073 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p-pStart]);CHKERRQ(ierr); 2074 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 2075 ierr = ISRestoreIndices(patch->gtolWithAll, >olArrayWithAll);CHKERRQ(ierr); 2076 } 2077 } 2078 if (patch->save_operators) { 2079 ierr = PetscMalloc1(patch->npatch, &patch->mat);CHKERRQ(ierr); 2080 for (i = 0; i < patch->npatch; ++i) { 2081 ierr = PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE);CHKERRQ(ierr); 2082 } 2083 } 2084 ierr = PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 2085 2086 /* If desired, calculate weights for dof multiplicity */ 2087 if (patch->partition_of_unity) { 2088 PetscScalar *input = NULL; 2089 PetscScalar *output = NULL; 2090 Vec global; 2091 2092 ierr = VecDuplicate(patch->localRHS, &patch->dof_weights);CHKERRQ(ierr); 2093 if(patch->local_composition_type == PC_COMPOSITE_ADDITIVE) { 2094 for (i = 0; i < patch->npatch; ++i) { 2095 PetscInt dof; 2096 2097 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);CHKERRQ(ierr); 2098 if (dof <= 0) continue; 2099 ierr = VecSet(patch->patchRHS[i], 1.0);CHKERRQ(ierr); 2100 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchRHS[i], patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);CHKERRQ(ierr); 2101 } 2102 } else { 2103 /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */ 2104 ierr = VecSet(patch->dof_weights, 1.0);CHKERRQ(ierr); 2105 } 2106 2107 VecDuplicate(patch->dof_weights, &global); 2108 VecSet(global, 0.); 2109 2110 ierr = VecGetArray(patch->dof_weights, &input);CHKERRQ(ierr); 2111 ierr = VecGetArray(global, &output);CHKERRQ(ierr); 2112 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);CHKERRQ(ierr); 2113 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);CHKERRQ(ierr); 2114 ierr = VecRestoreArray(patch->dof_weights, &input);CHKERRQ(ierr); 2115 ierr = VecRestoreArray(global, &output);CHKERRQ(ierr); 2116 2117 ierr = VecReciprocal(global);CHKERRQ(ierr); 2118 2119 ierr = VecGetArray(patch->dof_weights, &output);CHKERRQ(ierr); 2120 ierr = VecGetArray(global, &input);CHKERRQ(ierr); 2121 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, input, output);CHKERRQ(ierr); 2122 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, input, output);CHKERRQ(ierr); 2123 ierr = VecRestoreArray(patch->dof_weights, &output);CHKERRQ(ierr); 2124 ierr = VecRestoreArray(global, &input);CHKERRQ(ierr); 2125 ierr = VecDestroy(&global);CHKERRQ(ierr); 2126 } 2127 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators) { 2128 ierr = PetscMalloc1(patch->npatch, &patch->matWithArtificial);CHKERRQ(ierr); 2129 } 2130 } 2131 ierr = (*patch->setupsolver)(pc);CHKERRQ(ierr); 2132 PetscFunctionReturn(0); 2133 } 2134 2135 static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y) 2136 { 2137 PC_PATCH *patch = (PC_PATCH *) pc->data; 2138 KSP ksp = (KSP) patch->solver[i]; 2139 PetscErrorCode ierr; 2140 2141 PetscFunctionBegin; 2142 if (!patch->save_operators) { 2143 Mat mat; 2144 2145 ierr = PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE);CHKERRQ(ierr); 2146 /* Populate operator here. */ 2147 ierr = PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE);CHKERRQ(ierr); 2148 ierr = KSPSetOperators(ksp, mat, mat);CHKERRQ(ierr); 2149 /* Drop reference so the KSPSetOperators below will blow it away. */ 2150 ierr = MatDestroy(&mat);CHKERRQ(ierr); 2151 } 2152 ierr = PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 2153 if (!ksp->setfromoptionscalled) { 2154 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 2155 } 2156 ierr = KSPSolve(ksp, x, y);CHKERRQ(ierr); 2157 ierr = PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 2158 if (!patch->save_operators) { 2159 PC pc; 2160 ierr = KSPSetOperators(ksp, NULL, NULL);CHKERRQ(ierr); 2161 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 2162 /* Destroy PC context too, otherwise the factored matrix hangs around. */ 2163 ierr = PCReset(pc);CHKERRQ(ierr); 2164 } 2165 PetscFunctionReturn(0); 2166 } 2167 2168 static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart) 2169 { 2170 PC_PATCH *patch = (PC_PATCH *) pc->data; 2171 Mat multMat; 2172 PetscErrorCode ierr; 2173 2174 PetscFunctionBegin; 2175 2176 if (patch->save_operators) { 2177 multMat = patch->matWithArtificial[i]; 2178 } else { 2179 /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/ 2180 Mat matSquare; 2181 PetscInt dof; 2182 IS rowis; 2183 ierr = PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);CHKERRQ(ierr); 2184 ierr = MatZeroEntries(matSquare);CHKERRQ(ierr); 2185 ierr = PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);CHKERRQ(ierr); 2186 ierr = MatGetSize(matSquare, &dof, NULL);CHKERRQ(ierr); 2187 ierr = ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis); CHKERRQ(ierr); 2188 ierr = MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat); CHKERRQ(ierr); 2189 ierr = MatDestroy(&matSquare);CHKERRQ(ierr); 2190 ierr = ISDestroy(&rowis); CHKERRQ(ierr); 2191 } 2192 ierr = MatMult(multMat, patch->patchUpdate[i], patch->patchRHSWithArtificial[i]); CHKERRQ(ierr); 2193 ierr = VecScale(patch->patchRHSWithArtificial[i], -1.0); CHKERRQ(ierr); 2194 ierr = PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial[i], patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL); CHKERRQ(ierr); 2195 if (!patch->save_operators) { 2196 ierr = MatDestroy(&multMat); CHKERRQ(ierr); 2197 } 2198 PetscFunctionReturn(0); 2199 } 2200 2201 static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y) 2202 { 2203 PC_PATCH *patch = (PC_PATCH *) pc->data; 2204 const PetscScalar *globalRHS = NULL; 2205 PetscScalar *localRHS = NULL; 2206 PetscScalar *globalUpdate = NULL; 2207 const PetscInt *bcNodes = NULL; 2208 PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1; 2209 PetscInt start[2] = {0, 0}; 2210 PetscInt end[2] = {-1, -1}; 2211 const PetscInt inc[2] = {1, -1}; 2212 const PetscScalar *localUpdate; 2213 const PetscInt *iterationSet; 2214 PetscInt pStart, numBcs, n, sweep, bc, j; 2215 PetscErrorCode ierr; 2216 2217 PetscFunctionBegin; 2218 ierr = PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 2219 ierr = PetscOptionsPushGetViewerOff(PETSC_TRUE);CHKERRQ(ierr); 2220 /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */ 2221 end[0] = patch->npatch; 2222 start[1] = patch->npatch-1; 2223 if (patch->user_patches) { 2224 ierr = ISGetLocalSize(patch->iterationSet, &end[0]);CHKERRQ(ierr); 2225 start[1] = end[0] - 1; 2226 ierr = ISGetIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr); 2227 } 2228 /* Scatter from global space into overlapped local spaces */ 2229 ierr = VecGetArrayRead(x, &globalRHS);CHKERRQ(ierr); 2230 ierr = VecGetArray(patch->localRHS, &localRHS);CHKERRQ(ierr); 2231 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, globalRHS, localRHS);CHKERRQ(ierr); 2232 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, globalRHS, localRHS);CHKERRQ(ierr); 2233 ierr = VecRestoreArrayRead(x, &globalRHS);CHKERRQ(ierr); 2234 ierr = VecRestoreArray(patch->localRHS, &localRHS);CHKERRQ(ierr); 2235 2236 ierr = VecSet(patch->localUpdate, 0.0);CHKERRQ(ierr); 2237 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);CHKERRQ(ierr); 2238 for (sweep = 0; sweep < nsweep; sweep++) { 2239 for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) { 2240 PetscInt i = patch->user_patches ? iterationSet[j] : j; 2241 PetscInt start, len; 2242 2243 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);CHKERRQ(ierr); 2244 ierr = PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);CHKERRQ(ierr); 2245 /* TODO: Squash out these guys in the setup as well. */ 2246 if (len <= 0) continue; 2247 /* TODO: Do we need different scatters for X and Y? */ 2248 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localRHS, patch->patchRHS[i], INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR);CHKERRQ(ierr); 2249 ierr = (*patch->applysolver)(pc, i, patch->patchRHS[i], patch->patchUpdate[i]);CHKERRQ(ierr); 2250 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchUpdate[i], patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);CHKERRQ(ierr); 2251 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 2252 ierr = (*patch->updatemultiplicative)(pc, i, pStart);CHKERRQ(ierr); 2253 } 2254 } 2255 } 2256 if (patch->user_patches) {ierr = ISRestoreIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr);} 2257 /* XXX: should we do this on the global vector? */ 2258 if (patch->partition_of_unity) { 2259 ierr = VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights);CHKERRQ(ierr); 2260 } 2261 /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */ 2262 ierr = VecSet(y, 0.0);CHKERRQ(ierr); 2263 ierr = VecGetArray(y, &globalUpdate);CHKERRQ(ierr); 2264 ierr = VecGetArrayRead(patch->localUpdate, &localUpdate);CHKERRQ(ierr); 2265 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);CHKERRQ(ierr); 2266 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);CHKERRQ(ierr); 2267 ierr = VecRestoreArrayRead(patch->localUpdate, &localUpdate);CHKERRQ(ierr); 2268 2269 /* Now we need to send the global BC values through */ 2270 ierr = VecGetArrayRead(x, &globalRHS);CHKERRQ(ierr); 2271 ierr = ISGetSize(patch->globalBcNodes, &numBcs);CHKERRQ(ierr); 2272 ierr = ISGetIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 2273 ierr = VecGetLocalSize(x, &n);CHKERRQ(ierr); 2274 for (bc = 0; bc < numBcs; ++bc) { 2275 const PetscInt idx = bcNodes[bc]; 2276 if (idx < n) globalUpdate[idx] = globalRHS[idx]; 2277 } 2278 2279 ierr = ISRestoreIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 2280 ierr = VecRestoreArrayRead(x, &globalRHS);CHKERRQ(ierr); 2281 ierr = VecRestoreArray(y, &globalUpdate);CHKERRQ(ierr); 2282 2283 ierr = PetscOptionsPopGetViewerOff();CHKERRQ(ierr); 2284 ierr = PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 2285 PetscFunctionReturn(0); 2286 } 2287 2288 static PetscErrorCode PCReset_PATCH_Linear(PC pc) 2289 { 2290 PC_PATCH *patch = (PC_PATCH *) pc->data; 2291 PetscInt i; 2292 PetscErrorCode ierr; 2293 2294 PetscFunctionBegin; 2295 if (patch->solver) { 2296 for (i = 0; i < patch->npatch; ++i) {ierr = KSPReset((KSP) patch->solver[i]);CHKERRQ(ierr);} 2297 } 2298 PetscFunctionReturn(0); 2299 } 2300 2301 static PetscErrorCode PCReset_PATCH(PC pc) 2302 { 2303 PC_PATCH *patch = (PC_PATCH *) pc->data; 2304 PetscInt i; 2305 PetscErrorCode ierr; 2306 2307 PetscFunctionBegin; 2308 /* TODO: Get rid of all these ifs */ 2309 ierr = PetscSFDestroy(&patch->defaultSF);CHKERRQ(ierr); 2310 ierr = PetscSectionDestroy(&patch->cellCounts);CHKERRQ(ierr); 2311 ierr = PetscSectionDestroy(&patch->pointCounts);CHKERRQ(ierr); 2312 ierr = PetscSectionDestroy(&patch->cellNumbering);CHKERRQ(ierr); 2313 ierr = PetscSectionDestroy(&patch->gtolCounts);CHKERRQ(ierr); 2314 ierr = ISDestroy(&patch->gtol);CHKERRQ(ierr); 2315 ierr = ISDestroy(&patch->cells);CHKERRQ(ierr); 2316 ierr = ISDestroy(&patch->points);CHKERRQ(ierr); 2317 ierr = ISDestroy(&patch->dofs);CHKERRQ(ierr); 2318 ierr = ISDestroy(&patch->offs);CHKERRQ(ierr); 2319 ierr = PetscSectionDestroy(&patch->patchSection);CHKERRQ(ierr); 2320 ierr = ISDestroy(&patch->ghostBcNodes);CHKERRQ(ierr); 2321 ierr = ISDestroy(&patch->globalBcNodes);CHKERRQ(ierr); 2322 ierr = PetscSectionDestroy(&patch->gtolCountsWithArtificial);CHKERRQ(ierr); 2323 ierr = ISDestroy(&patch->gtolWithArtificial);CHKERRQ(ierr); 2324 ierr = ISDestroy(&patch->dofsWithArtificial);CHKERRQ(ierr); 2325 ierr = ISDestroy(&patch->offsWithArtificial);CHKERRQ(ierr); 2326 ierr = PetscSectionDestroy(&patch->gtolCountsWithAll);CHKERRQ(ierr); 2327 ierr = ISDestroy(&patch->gtolWithAll);CHKERRQ(ierr); 2328 ierr = ISDestroy(&patch->dofsWithAll);CHKERRQ(ierr); 2329 ierr = ISDestroy(&patch->offsWithAll);CHKERRQ(ierr); 2330 2331 2332 if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscSectionDestroy(&patch->dofSection[i]);CHKERRQ(ierr);} 2333 ierr = PetscFree(patch->dofSection);CHKERRQ(ierr); 2334 ierr = PetscFree(patch->bs);CHKERRQ(ierr); 2335 ierr = PetscFree(patch->nodesPerCell);CHKERRQ(ierr); 2336 if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscFree(patch->cellNodeMap[i]);CHKERRQ(ierr);} 2337 ierr = PetscFree(patch->cellNodeMap);CHKERRQ(ierr); 2338 ierr = PetscFree(patch->subspaceOffsets);CHKERRQ(ierr); 2339 2340 ierr = (*patch->resetsolver)(pc);CHKERRQ(ierr); 2341 2342 if (patch->subspaces_to_exclude) { 2343 PetscHSetIDestroy(&patch->subspaces_to_exclude); 2344 } 2345 2346 ierr = VecDestroy(&patch->localRHS);CHKERRQ(ierr); 2347 ierr = VecDestroy(&patch->localUpdate);CHKERRQ(ierr); 2348 if (patch->patchRHS) { 2349 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchRHS[i]);CHKERRQ(ierr);} 2350 ierr = PetscFree(patch->patchRHS);CHKERRQ(ierr); 2351 } 2352 if (patch->patchUpdate) { 2353 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchUpdate[i]);CHKERRQ(ierr);} 2354 ierr = PetscFree(patch->patchUpdate);CHKERRQ(ierr); 2355 } 2356 ierr = VecDestroy(&patch->dof_weights);CHKERRQ(ierr); 2357 if (patch->patch_dof_weights) { 2358 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patch_dof_weights[i]);CHKERRQ(ierr);} 2359 ierr = PetscFree(patch->patch_dof_weights);CHKERRQ(ierr); 2360 } 2361 if (patch->mat) { 2362 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->mat[i]);CHKERRQ(ierr);} 2363 ierr = PetscFree(patch->mat);CHKERRQ(ierr); 2364 } 2365 if (patch->matWithArtificial) { 2366 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->matWithArtificial[i]);CHKERRQ(ierr);} 2367 ierr = PetscFree(patch->matWithArtificial);CHKERRQ(ierr); 2368 } 2369 if (patch->patchRHSWithArtificial) { 2370 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchRHSWithArtificial[i]);CHKERRQ(ierr);} 2371 ierr = PetscFree(patch->patchRHSWithArtificial);CHKERRQ(ierr); 2372 } 2373 if(patch->dofMappingWithoutToWithArtificial) { 2374 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]);CHKERRQ(ierr);} 2375 ierr = PetscFree(patch->dofMappingWithoutToWithArtificial);CHKERRQ(ierr); 2376 2377 } 2378 if(patch->dofMappingWithoutToWithAll) { 2379 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->dofMappingWithoutToWithAll[i]);CHKERRQ(ierr);} 2380 ierr = PetscFree(patch->dofMappingWithoutToWithAll);CHKERRQ(ierr); 2381 2382 } 2383 ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr); 2384 if (patch->userIS) { 2385 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->userIS[i]);CHKERRQ(ierr);} 2386 ierr = PetscFree(patch->userIS);CHKERRQ(ierr); 2387 } 2388 patch->bs = 0; 2389 patch->cellNodeMap = NULL; 2390 patch->nsubspaces = 0; 2391 ierr = ISDestroy(&patch->iterationSet);CHKERRQ(ierr); 2392 2393 ierr = PetscViewerDestroy(&patch->viewerSection);CHKERRQ(ierr); 2394 PetscFunctionReturn(0); 2395 } 2396 2397 static PetscErrorCode PCDestroy_PATCH_Linear(PC pc) 2398 { 2399 PC_PATCH *patch = (PC_PATCH *) pc->data; 2400 PetscInt i; 2401 PetscErrorCode ierr; 2402 2403 PetscFunctionBegin; 2404 if (patch->solver) { 2405 for (i = 0; i < patch->npatch; ++i) {ierr = KSPDestroy((KSP *) &patch->solver[i]);CHKERRQ(ierr);} 2406 ierr = PetscFree(patch->solver);CHKERRQ(ierr); 2407 } 2408 PetscFunctionReturn(0); 2409 } 2410 2411 static PetscErrorCode PCDestroy_PATCH(PC pc) 2412 { 2413 PC_PATCH *patch = (PC_PATCH *) pc->data; 2414 PetscErrorCode ierr; 2415 2416 PetscFunctionBegin; 2417 ierr = PCReset_PATCH(pc);CHKERRQ(ierr); 2418 ierr = (*patch->destroysolver)(pc);CHKERRQ(ierr); 2419 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2420 PetscFunctionReturn(0); 2421 } 2422 2423 static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc) 2424 { 2425 PC_PATCH *patch = (PC_PATCH *) pc->data; 2426 PCPatchConstructType patchConstructionType = PC_PATCH_STAR; 2427 char sub_mat_type[PETSC_MAX_PATH_LEN]; 2428 char option[PETSC_MAX_PATH_LEN]; 2429 const char *prefix; 2430 PetscBool flg, dimflg, codimflg; 2431 MPI_Comm comm; 2432 PetscInt *ifields, nfields, k; 2433 PetscErrorCode ierr; 2434 PCCompositeType loctype = PC_COMPOSITE_ADDITIVE; 2435 2436 PetscFunctionBegin; 2437 ierr = PetscObjectGetComm((PetscObject) pc, &comm);CHKERRQ(ierr); 2438 ierr = PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);CHKERRQ(ierr); 2439 ierr = PetscOptionsHead(PetscOptionsObject, "Patch solver options");CHKERRQ(ierr); 2440 2441 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_save_operators", patch->classname); 2442 ierr = PetscOptionsBool(option, "Store all patch operators for lifetime of object?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg);CHKERRQ(ierr); 2443 2444 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_partition_of_unity", patch->classname); 2445 ierr = PetscOptionsBool(option, "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg);CHKERRQ(ierr); 2446 2447 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname); 2448 ierr = PetscOptionsEnum(option,"Type of local solver composition (additive or multiplicative)","PCPatchSetLocalComposition",PCCompositeTypes,(PetscEnum)loctype,(PetscEnum*)&loctype,&flg);CHKERRQ(ierr); 2449 if(flg) { ierr = PCPatchSetLocalComposition(pc, loctype);CHKERRQ(ierr);} 2450 2451 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname); 2452 ierr = PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg);CHKERRQ(ierr); 2453 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname); 2454 ierr = PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg);CHKERRQ(ierr); 2455 if (dimflg && codimflg) {SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");CHKERRQ(ierr);} 2456 2457 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_type", patch->classname); 2458 ierr = PetscOptionsEnum(option, "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum) patchConstructionType, (PetscEnum *) &patchConstructionType, &flg);CHKERRQ(ierr); 2459 if (flg) {ierr = PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL);CHKERRQ(ierr);} 2460 2461 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_vanka_dim", patch->classname); 2462 ierr = PetscOptionsInt(option, "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg);CHKERRQ(ierr); 2463 2464 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_ignore_dim", patch->classname); 2465 ierr = PetscOptionsInt(option, "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg);CHKERRQ(ierr); 2466 2467 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_pardecomp_overlap", patch->classname); 2468 ierr = PetscOptionsInt(option, "What overlap should we use in construct type pardecomp?", "PCPATCH", patch->pardecomp_overlap, &patch->pardecomp_overlap, &flg);CHKERRQ(ierr); 2469 2470 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_sub_mat_type", patch->classname); 2471 ierr = PetscOptionsFList(option, "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg);CHKERRQ(ierr); 2472 if (flg) {ierr = PCPatchSetSubMatType(pc, sub_mat_type);CHKERRQ(ierr);} 2473 2474 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_symmetrise_sweep", patch->classname); 2475 ierr = PetscOptionsBool(option, "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg);CHKERRQ(ierr); 2476 2477 /* If the user has set the number of subspaces, use that for the buffer size, 2478 otherwise use a large number */ 2479 if (patch->nsubspaces <= 0) { 2480 nfields = 128; 2481 } else { 2482 nfields = patch->nsubspaces; 2483 } 2484 ierr = PetscMalloc1(nfields, &ifields);CHKERRQ(ierr); 2485 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname); 2486 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,option,ifields,&nfields,&flg);CHKERRQ(ierr); 2487 if (flg && (patchConstructionType == PC_PATCH_USER)) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "We cannot support excluding a subspace with user patches because we do not index patches with a mesh point"); 2488 if (flg) { 2489 PetscHSetIClear(patch->subspaces_to_exclude); 2490 for (k = 0; k < nfields; k++) { 2491 PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]); 2492 } 2493 } 2494 ierr = PetscFree(ifields);CHKERRQ(ierr); 2495 2496 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname); 2497 ierr = PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);CHKERRQ(ierr); 2498 2499 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname); 2500 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells);CHKERRQ(ierr); 2501 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname); 2502 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);CHKERRQ(ierr); 2503 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname); 2504 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection);CHKERRQ(ierr); 2505 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname); 2506 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);CHKERRQ(ierr); 2507 ierr = PetscOptionsTail();CHKERRQ(ierr); 2508 patch->optionsSet = PETSC_TRUE; 2509 PetscFunctionReturn(0); 2510 } 2511 2512 static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc) 2513 { 2514 PC_PATCH *patch = (PC_PATCH*) pc->data; 2515 KSPConvergedReason reason; 2516 PetscInt i; 2517 PetscErrorCode ierr; 2518 2519 PetscFunctionBegin; 2520 if (!patch->save_operators) { 2521 /* Can't do this here because the sub KSPs don't have an operator attached yet. */ 2522 PetscFunctionReturn(0); 2523 } 2524 for (i = 0; i < patch->npatch; ++i) { 2525 if (!((KSP) patch->solver[i])->setfromoptionscalled) { 2526 ierr = KSPSetFromOptions((KSP) patch->solver[i]);CHKERRQ(ierr); 2527 } 2528 ierr = KSPSetUp((KSP) patch->solver[i]);CHKERRQ(ierr); 2529 ierr = KSPGetConvergedReason((KSP) patch->solver[i], &reason);CHKERRQ(ierr); 2530 if (reason == KSP_DIVERGED_PCSETUP_FAILED) pc->failedreason = PC_SUBPC_ERROR; 2531 } 2532 PetscFunctionReturn(0); 2533 } 2534 2535 static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer) 2536 { 2537 PC_PATCH *patch = (PC_PATCH *) pc->data; 2538 PetscViewer sviewer; 2539 PetscBool isascii; 2540 PetscMPIInt rank; 2541 PetscErrorCode ierr; 2542 2543 PetscFunctionBegin; 2544 /* TODO Redo tabbing with set tbas in new style */ 2545 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 2546 if (!isascii) PetscFunctionReturn(0); 2547 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);CHKERRQ(ierr); 2548 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 2549 ierr = PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);CHKERRQ(ierr); 2550 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 2551 ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n");CHKERRQ(ierr); 2552 } else { 2553 ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");CHKERRQ(ierr); 2554 } 2555 if (patch->partition_of_unity) {ierr = PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");CHKERRQ(ierr);} 2556 else {ierr = PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");CHKERRQ(ierr);} 2557 if (patch->symmetrise_sweep) {ierr = PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");CHKERRQ(ierr);} 2558 else {ierr = PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");CHKERRQ(ierr);} 2559 if (!patch->save_operators) {ierr = PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");CHKERRQ(ierr);} 2560 else {ierr = PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");CHKERRQ(ierr);} 2561 if (patch->patchconstructop == PCPatchConstruct_Star) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");CHKERRQ(ierr);} 2562 else if (patch->patchconstructop == PCPatchConstruct_Vanka) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");CHKERRQ(ierr);} 2563 else if (patch->patchconstructop == PCPatchConstruct_User) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");CHKERRQ(ierr);} 2564 else {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");CHKERRQ(ierr);} 2565 ierr = PetscViewerASCIIPrintf(viewer, "Solver on patches (all same):\n");CHKERRQ(ierr); 2566 if (patch->solver) { 2567 ierr = PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 2568 if (!rank) { 2569 ierr = PetscViewerASCIIPushTab(sviewer);CHKERRQ(ierr); 2570 ierr = PetscObjectView(patch->solver[0], sviewer);CHKERRQ(ierr); 2571 ierr = PetscViewerASCIIPopTab(sviewer);CHKERRQ(ierr); 2572 } 2573 ierr = PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 2574 } else { 2575 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 2576 ierr = PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n");CHKERRQ(ierr); 2577 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 2578 } 2579 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 2580 PetscFunctionReturn(0); 2581 } 2582 2583 /*MC 2584 PCPATCH - A PC object that encapsulates flexible definition of blocks for overlapping and non-overlapping 2585 small block additive preconditioners. Block definition is based on topology from 2586 a DM and equation numbering from a PetscSection. 2587 2588 Options Database Keys: 2589 + -pc_patch_cells_view - Views the process local cell numbers for each patch 2590 . -pc_patch_points_view - Views the process local mesh point numbers for each patch 2591 . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch 2592 . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary 2593 - -pc_patch_sub_mat_view - Views the matrix associated with each patch 2594 2595 Level: intermediate 2596 2597 .seealso: PCType, PCCreate(), PCSetType() 2598 M*/ 2599 PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc) 2600 { 2601 PC_PATCH *patch; 2602 PetscErrorCode ierr; 2603 2604 PetscFunctionBegin; 2605 ierr = PetscNewLog(pc, &patch);CHKERRQ(ierr); 2606 2607 if (patch->subspaces_to_exclude) { 2608 PetscHSetIDestroy(&patch->subspaces_to_exclude); 2609 } 2610 PetscHSetICreate(&patch->subspaces_to_exclude); 2611 2612 patch->classname = "pc"; 2613 patch->isNonlinear = PETSC_FALSE; 2614 2615 /* Set some defaults */ 2616 patch->combined = PETSC_FALSE; 2617 patch->save_operators = PETSC_TRUE; 2618 patch->local_composition_type = PC_COMPOSITE_ADDITIVE; 2619 patch->partition_of_unity = PETSC_FALSE; 2620 patch->codim = -1; 2621 patch->dim = -1; 2622 patch->vankadim = -1; 2623 patch->ignoredim = -1; 2624 patch->pardecomp_overlap = 0; 2625 patch->patchconstructop = PCPatchConstruct_Star; 2626 patch->symmetrise_sweep = PETSC_FALSE; 2627 patch->npatch = 0; 2628 patch->userIS = NULL; 2629 patch->optionsSet = PETSC_FALSE; 2630 patch->iterationSet = NULL; 2631 patch->user_patches = PETSC_FALSE; 2632 ierr = PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 2633 patch->viewPatches = PETSC_FALSE; 2634 patch->viewCells = PETSC_FALSE; 2635 patch->viewPoints = PETSC_FALSE; 2636 patch->viewSection = PETSC_FALSE; 2637 patch->viewMatrix = PETSC_FALSE; 2638 patch->setupsolver = PCSetUp_PATCH_Linear; 2639 patch->applysolver = PCApply_PATCH_Linear; 2640 patch->resetsolver = PCReset_PATCH_Linear; 2641 patch->destroysolver = PCDestroy_PATCH_Linear; 2642 patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear; 2643 patch->dofMappingWithoutToWithArtificial = NULL; 2644 patch->dofMappingWithoutToWithAll = NULL; 2645 2646 pc->data = (void *) patch; 2647 pc->ops->apply = PCApply_PATCH; 2648 pc->ops->applytranspose = 0; /* PCApplyTranspose_PATCH; */ 2649 pc->ops->setup = PCSetUp_PATCH; 2650 pc->ops->reset = PCReset_PATCH; 2651 pc->ops->destroy = PCDestroy_PATCH; 2652 pc->ops->setfromoptions = PCSetFromOptions_PATCH; 2653 pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH; 2654 pc->ops->view = PCView_PATCH; 2655 pc->ops->applyrichardson = 0; 2656 2657 PetscFunctionReturn(0); 2658 } 2659