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