1 #include <petsc/private/pcpatchimpl.h> /*I "petscpc.h" I*/ 2 #include <petsc/private/dmpleximpl.h> /* For DMPlexComputeJacobian_Patch_Internal() */ 3 #include <petscsf.h> 4 #include <petscbt.h> 5 #include <petscds.h> 6 7 PetscLogEvent PC_Patch_CreatePatches, PC_Patch_ComputeOp, PC_Patch_Solve, PC_Patch_Scatter, PC_Patch_Apply, PC_Patch_Prealloc; 8 9 PETSC_STATIC_INLINE PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format) 10 { 11 PetscErrorCode ierr; 12 13 ierr = PetscViewerPushFormat(viewer, format);CHKERRQ(ierr); 14 ierr = PetscObjectView(obj, viewer);CHKERRQ(ierr); 15 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 16 return(0); 17 } 18 19 static PetscErrorCode PCPatchConstruct_Star(void *vpatch, DM dm, PetscInt point, PetscHSetI ht) 20 { 21 PetscInt starSize; 22 PetscInt *star = NULL, si; 23 PetscErrorCode ierr; 24 25 PetscFunctionBegin; 26 PetscHSetIClear(ht); 27 /* To start with, add the point we care about */ 28 ierr = PetscHSetIAdd(ht, point);CHKERRQ(ierr); 29 /* Loop over all the points that this point connects to */ 30 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 31 for (si = 0; si < starSize*2; si += 2) {ierr = PetscHSetIAdd(ht, star[si]);CHKERRQ(ierr);} 32 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 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 PetscErrorCode ierr; 44 45 PetscFunctionBegin; 46 ierr = PetscHSetIClear(ht);CHKERRQ(ierr); 47 /* To start with, add the point we care about */ 48 ierr = PetscHSetIAdd(ht, point);CHKERRQ(ierr); 49 /* Should we ignore any points of a certain dimension? */ 50 if (patch->vankadim >= 0) { 51 shouldIgnore = PETSC_TRUE; 52 ierr = DMPlexGetDepthStratum(dm, patch->vankadim, &iStart, &iEnd);CHKERRQ(ierr); 53 } 54 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 55 /* Loop over all the cells that this point connects to */ 56 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 57 for (si = 0; si < starSize*2; si += 2) { 58 const PetscInt cell = star[si]; 59 PetscInt closureSize; 60 PetscInt *closure = NULL, ci; 61 62 if (cell < cStart || cell >= cEnd) continue; 63 /* now loop over all entities in the closure of that cell */ 64 ierr = DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 65 for (ci = 0; ci < closureSize*2; ci += 2) { 66 const PetscInt newpoint = closure[ci]; 67 68 /* We've been told to ignore entities of this type.*/ 69 if (shouldIgnore && newpoint >= iStart && newpoint < iEnd) continue; 70 ierr = PetscHSetIAdd(ht, newpoint);CHKERRQ(ierr); 71 } 72 ierr = DMPlexRestoreTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 73 } 74 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 /* The user's already set the patches in patch->userIS. Build the hash tables */ 79 static PetscErrorCode PCPatchConstruct_User(void *vpatch, DM dm, PetscInt point, PetscHSetI ht) 80 { 81 PC_PATCH *patch = (PC_PATCH *) vpatch; 82 IS patchis = patch->userIS[point]; 83 PetscInt n; 84 const PetscInt *patchdata; 85 PetscInt pStart, pEnd, i; 86 PetscErrorCode ierr; 87 88 PetscFunctionBegin; 89 ierr = PetscHSetIClear(ht);CHKERRQ(ierr); 90 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 91 ierr = ISGetLocalSize(patchis, &n);CHKERRQ(ierr); 92 ierr = ISGetIndices(patchis, &patchdata);CHKERRQ(ierr); 93 for (i = 0; i < n; ++i) { 94 const PetscInt ownedpoint = patchdata[i]; 95 96 if (ownedpoint < pStart || ownedpoint >= pEnd) { 97 SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D was not in [%D, %D)", ownedpoint, pStart, pEnd); 98 } 99 ierr = PetscHSetIAdd(ht, ownedpoint);CHKERRQ(ierr); 100 } 101 ierr = ISRestoreIndices(patchis, &patchdata);CHKERRQ(ierr); 102 PetscFunctionReturn(0); 103 } 104 105 static PetscErrorCode PCPatchCreateDefaultSF_Private(PC pc, PetscInt n, const PetscSF *sf, const PetscInt *bs) 106 { 107 PC_PATCH *patch = (PC_PATCH *) pc->data; 108 PetscInt i; 109 PetscErrorCode ierr; 110 111 PetscFunctionBegin; 112 if (n == 1 && bs[0] == 1) { 113 patch->defaultSF = sf[0]; 114 ierr = PetscObjectReference((PetscObject) patch->defaultSF);CHKERRQ(ierr); 115 } else { 116 PetscInt allRoots = 0, allLeaves = 0; 117 PetscInt leafOffset = 0; 118 PetscInt *ilocal = NULL; 119 PetscSFNode *iremote = NULL; 120 PetscInt *remoteOffsets = NULL; 121 PetscInt index = 0; 122 PetscHMapI rankToIndex; 123 PetscInt numRanks = 0; 124 PetscSFNode *remote = NULL; 125 PetscSF rankSF; 126 PetscInt *ranks = NULL; 127 PetscInt *offsets = NULL; 128 MPI_Datatype contig; 129 PetscHSetI ranksUniq; 130 131 /* First figure out how many dofs there are in the concatenated numbering. 132 * allRoots: number of owned global dofs; 133 * allLeaves: number of visible dofs (global + ghosted). 134 */ 135 for (i = 0; i < n; ++i) { 136 PetscInt nroots, nleaves; 137 138 ierr = PetscSFGetGraph(sf[i], &nroots, &nleaves, NULL, NULL);CHKERRQ(ierr); 139 allRoots += nroots * bs[i]; 140 allLeaves += nleaves * bs[i]; 141 } 142 ierr = PetscMalloc1(allLeaves, &ilocal);CHKERRQ(ierr); 143 ierr = PetscMalloc1(allLeaves, &iremote);CHKERRQ(ierr); 144 /* Now build an SF that just contains process connectivity. */ 145 ierr = PetscHSetICreate(&ranksUniq);CHKERRQ(ierr); 146 for (i = 0; i < n; ++i) { 147 const PetscMPIInt *ranks = NULL; 148 PetscInt nranks, j; 149 150 ierr = PetscSFSetUp(sf[i]);CHKERRQ(ierr); 151 ierr = PetscSFGetRanks(sf[i], &nranks, &ranks, NULL, NULL, NULL);CHKERRQ(ierr); 152 /* These are all the ranks who communicate with me. */ 153 for (j = 0; j < nranks; ++j) { 154 ierr = PetscHSetIAdd(ranksUniq, (PetscInt) ranks[j]);CHKERRQ(ierr); 155 } 156 } 157 ierr = PetscHSetIGetSize(ranksUniq, &numRanks);CHKERRQ(ierr); 158 ierr = PetscMalloc1(numRanks, &remote);CHKERRQ(ierr); 159 ierr = PetscMalloc1(numRanks, &ranks);CHKERRQ(ierr); 160 ierr = PetscHSetIGetElems(ranksUniq, &index, ranks);CHKERRQ(ierr); 161 162 ierr = PetscHMapICreate(&rankToIndex);CHKERRQ(ierr); 163 for (i = 0; i < numRanks; ++i) { 164 remote[i].rank = ranks[i]; 165 remote[i].index = 0; 166 ierr = PetscHMapISet(rankToIndex, ranks[i], i);CHKERRQ(ierr); 167 } 168 ierr = PetscFree(ranks);CHKERRQ(ierr); 169 ierr = PetscHSetIDestroy(&ranksUniq);CHKERRQ(ierr); 170 ierr = PetscSFCreate(PetscObjectComm((PetscObject) pc), &rankSF);CHKERRQ(ierr); 171 ierr = PetscSFSetGraph(rankSF, 1, numRanks, NULL, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);CHKERRQ(ierr); 172 ierr = PetscSFSetUp(rankSF);CHKERRQ(ierr); 173 /* OK, use it to communicate the root offset on the remote 174 * processes for each subspace. */ 175 ierr = PetscMalloc1(n, &offsets);CHKERRQ(ierr); 176 ierr = PetscMalloc1(n*numRanks, &remoteOffsets);CHKERRQ(ierr); 177 178 offsets[0] = 0; 179 for (i = 1; i < n; ++i) { 180 PetscInt nroots; 181 182 ierr = PetscSFGetGraph(sf[i-1], &nroots, NULL, NULL, NULL);CHKERRQ(ierr); 183 offsets[i] = offsets[i-1] + nroots*bs[i-1]; 184 } 185 /* Offsets are the offsets on the current process of the 186 * global dof numbering for the subspaces. */ 187 ierr = MPI_Type_contiguous(n, MPIU_INT, &contig);CHKERRQ(ierr); 188 ierr = MPI_Type_commit(&contig);CHKERRQ(ierr); 189 190 ierr = PetscSFBcastBegin(rankSF, contig, offsets, remoteOffsets);CHKERRQ(ierr); 191 ierr = PetscSFBcastEnd(rankSF, contig, offsets, remoteOffsets);CHKERRQ(ierr); 192 ierr = MPI_Type_free(&contig);CHKERRQ(ierr); 193 ierr = PetscFree(offsets);CHKERRQ(ierr); 194 ierr = PetscSFDestroy(&rankSF);CHKERRQ(ierr); 195 /* Now remoteOffsets contains the offsets on the remote 196 * processes who communicate with me. So now we can 197 * concatenate the list of SFs into a single one. */ 198 index = 0; 199 for (i = 0; i < n; ++i) { 200 const PetscSFNode *remote = NULL; 201 const PetscInt *local = NULL; 202 PetscInt nroots, nleaves, j; 203 204 ierr = PetscSFGetGraph(sf[i], &nroots, &nleaves, &local, &remote);CHKERRQ(ierr); 205 for (j = 0; j < nleaves; ++j) { 206 PetscInt rank = remote[j].rank; 207 PetscInt idx, rootOffset, k; 208 209 ierr = PetscHMapIGet(rankToIndex, rank, &idx);CHKERRQ(ierr); 210 if (idx == -1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Didn't find rank, huh?"); 211 /* Offset on given rank for ith subspace */ 212 rootOffset = remoteOffsets[n*idx + i]; 213 for (k = 0; k < bs[i]; ++k) { 214 ilocal[index] = (local ? local[j] : j)*bs[i] + k + leafOffset; 215 iremote[index].rank = remote[j].rank; 216 iremote[index].index = remote[j].index*bs[i] + k + rootOffset; 217 ++index; 218 } 219 } 220 leafOffset += nleaves * bs[i]; 221 } 222 ierr = PetscHMapIDestroy(&rankToIndex);CHKERRQ(ierr); 223 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 224 ierr = PetscSFCreate(PetscObjectComm((PetscObject)pc), &patch->defaultSF);CHKERRQ(ierr); 225 ierr = PetscSFSetGraph(patch->defaultSF, allRoots, allLeaves, ilocal, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER);CHKERRQ(ierr); 226 } 227 PetscFunctionReturn(0); 228 } 229 230 /* TODO: Docs */ 231 PetscErrorCode PCPatchSetIgnoreDim(PC pc, PetscInt dim) 232 { 233 PC_PATCH *patch = (PC_PATCH *) pc->data; 234 PetscFunctionBegin; 235 patch->ignoredim = dim; 236 PetscFunctionReturn(0); 237 } 238 239 /* TODO: Docs */ 240 PetscErrorCode PCPatchGetIgnoreDim(PC pc, PetscInt *dim) 241 { 242 PC_PATCH *patch = (PC_PATCH *) pc->data; 243 PetscFunctionBegin; 244 *dim = patch->ignoredim; 245 PetscFunctionReturn(0); 246 } 247 248 /* TODO: Docs */ 249 PetscErrorCode PCPatchSetSaveOperators(PC pc, PetscBool flg) 250 { 251 PC_PATCH *patch = (PC_PATCH *) pc->data; 252 PetscFunctionBegin; 253 patch->save_operators = flg; 254 PetscFunctionReturn(0); 255 } 256 257 /* TODO: Docs */ 258 PetscErrorCode PCPatchGetSaveOperators(PC pc, PetscBool *flg) 259 { 260 PC_PATCH *patch = (PC_PATCH *) pc->data; 261 PetscFunctionBegin; 262 *flg = patch->save_operators; 263 PetscFunctionReturn(0); 264 } 265 266 /* TODO: Docs */ 267 PetscErrorCode PCPatchSetPartitionOfUnity(PC pc, PetscBool flg) 268 { 269 PC_PATCH *patch = (PC_PATCH *) pc->data; 270 PetscFunctionBegin; 271 patch->partition_of_unity = flg; 272 PetscFunctionReturn(0); 273 } 274 275 /* TODO: Docs */ 276 PetscErrorCode PCPatchGetPartitionOfUnity(PC pc, PetscBool *flg) 277 { 278 PC_PATCH *patch = (PC_PATCH *) pc->data; 279 PetscFunctionBegin; 280 *flg = patch->partition_of_unity; 281 PetscFunctionReturn(0); 282 } 283 284 /* TODO: Docs */ 285 PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type) 286 { 287 PC_PATCH *patch = (PC_PATCH *) pc->data; 288 PetscErrorCode ierr; 289 290 PetscFunctionBegin; 291 if (patch->sub_mat_type) {ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr);} 292 ierr = PetscStrallocpy(sub_mat_type, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295 296 /* TODO: Docs */ 297 PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type) 298 { 299 PC_PATCH *patch = (PC_PATCH *) pc->data; 300 PetscFunctionBegin; 301 *sub_mat_type = patch->sub_mat_type; 302 PetscFunctionReturn(0); 303 } 304 305 /* TODO: Docs */ 306 PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering) 307 { 308 PC_PATCH *patch = (PC_PATCH *) pc->data; 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 patch->cellNumbering = cellNumbering; 313 ierr = PetscObjectReference((PetscObject) cellNumbering);CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316 317 /* TODO: Docs */ 318 PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering) 319 { 320 PC_PATCH *patch = (PC_PATCH *) pc->data; 321 PetscFunctionBegin; 322 *cellNumbering = patch->cellNumbering; 323 PetscFunctionReturn(0); 324 } 325 326 /* TODO: Docs */ 327 PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx) 328 { 329 PC_PATCH *patch = (PC_PATCH *) pc->data; 330 331 PetscFunctionBegin; 332 patch->ctype = ctype; 333 switch (ctype) { 334 case PC_PATCH_STAR: 335 patch->user_patches = PETSC_FALSE; 336 patch->patchconstructop = PCPatchConstruct_Star; 337 break; 338 case PC_PATCH_VANKA: 339 patch->user_patches = PETSC_FALSE; 340 patch->patchconstructop = PCPatchConstruct_Vanka; 341 break; 342 case PC_PATCH_USER: 343 case PC_PATCH_PYTHON: 344 patch->user_patches = PETSC_TRUE; 345 patch->patchconstructop = PCPatchConstruct_User; 346 if (func) { 347 patch->userpatchconstructionop = func; 348 patch->userpatchconstructctx = ctx; 349 } 350 break; 351 default: 352 SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype); 353 } 354 PetscFunctionReturn(0); 355 } 356 357 /* TODO: Docs */ 358 PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx) 359 { 360 PC_PATCH *patch = (PC_PATCH *) pc->data; 361 362 PetscFunctionBegin; 363 *ctype = patch->ctype; 364 switch (patch->ctype) { 365 case PC_PATCH_STAR: 366 case PC_PATCH_VANKA: 367 break; 368 case PC_PATCH_USER: 369 case PC_PATCH_PYTHON: 370 *func = patch->userpatchconstructionop; 371 *ctx = patch->userpatchconstructctx; 372 break; 373 default: 374 SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype); 375 } 376 PetscFunctionReturn(0); 377 } 378 379 /* TODO: Docs */ 380 PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, 381 const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 382 { 383 PC_PATCH *patch = (PC_PATCH *) pc->data; 384 DM dm; 385 PetscSF *sfs; 386 PetscInt cStart, cEnd, i, j; 387 PetscErrorCode ierr; 388 389 PetscFunctionBegin; 390 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 391 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 392 ierr = PetscMalloc1(nsubspaces, &sfs);CHKERRQ(ierr); 393 ierr = PetscMalloc1(nsubspaces, &patch->dofSection);CHKERRQ(ierr); 394 ierr = PetscMalloc1(nsubspaces, &patch->bs);CHKERRQ(ierr); 395 ierr = PetscMalloc1(nsubspaces, &patch->nodesPerCell);CHKERRQ(ierr); 396 ierr = PetscMalloc1(nsubspaces, &patch->cellNodeMap);CHKERRQ(ierr); 397 ierr = PetscMalloc1(nsubspaces+1, &patch->subspaceOffsets);CHKERRQ(ierr); 398 399 patch->nsubspaces = nsubspaces; 400 patch->totalDofsPerCell = 0; 401 for (i = 0; i < nsubspaces; ++i) { 402 ierr = DMGetDefaultSection(dms[i], &patch->dofSection[i]);CHKERRQ(ierr); 403 ierr = PetscObjectReference((PetscObject) patch->dofSection[i]);CHKERRQ(ierr); 404 ierr = DMGetDefaultSF(dms[i], &sfs[i]);CHKERRQ(ierr); 405 patch->bs[i] = bs[i]; 406 patch->nodesPerCell[i] = nodesPerCell[i]; 407 patch->totalDofsPerCell += nodesPerCell[i]*bs[i]; 408 ierr = PetscMalloc1((cEnd-cStart)*nodesPerCell[i]*bs[i], &patch->cellNodeMap[i]);CHKERRQ(ierr); 409 for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]*bs[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j]; 410 patch->subspaceOffsets[i] = subspaceOffsets[i]; 411 } 412 ierr = PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs);CHKERRQ(ierr); 413 ierr = PetscFree(sfs);CHKERRQ(ierr); 414 415 patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces]; 416 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);CHKERRQ(ierr); 417 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);CHKERRQ(ierr); 418 PetscFunctionReturn(0); 419 } 420 421 /* TODO: Docs */ 422 PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 423 { 424 PC_PATCH *patch = (PC_PATCH *) pc->data; 425 PetscInt cStart, cEnd, i, j; 426 PetscErrorCode ierr; 427 428 PetscFunctionBegin; 429 patch->combined = PETSC_TRUE; 430 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 431 ierr = DMGetNumFields(dm, &patch->nsubspaces);CHKERRQ(ierr); 432 ierr = PetscCalloc1(patch->nsubspaces, &patch->dofSection);CHKERRQ(ierr); 433 ierr = PetscMalloc1(patch->nsubspaces, &patch->bs);CHKERRQ(ierr); 434 ierr = PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell);CHKERRQ(ierr); 435 ierr = PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap);CHKERRQ(ierr); 436 ierr = PetscCalloc1(patch->nsubspaces+1, &patch->subspaceOffsets);CHKERRQ(ierr); 437 ierr = DMGetDefaultSection(dm, &patch->dofSection[0]);CHKERRQ(ierr); 438 ierr = PetscObjectReference((PetscObject) patch->dofSection[0]);CHKERRQ(ierr); 439 ierr = PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]);CHKERRQ(ierr); 440 patch->totalDofsPerCell = 0; 441 for (i = 0; i < patch->nsubspaces; ++i) { 442 patch->bs[i] = 1; 443 patch->nodesPerCell[i] = nodesPerCell[i]; 444 patch->totalDofsPerCell += nodesPerCell[i]; 445 ierr = PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);CHKERRQ(ierr); 446 for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j]; 447 } 448 ierr = DMGetDefaultSF(dm, &patch->defaultSF);CHKERRQ(ierr); 449 ierr = PetscObjectReference((PetscObject) patch->defaultSF);CHKERRQ(ierr); 450 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);CHKERRQ(ierr); 451 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);CHKERRQ(ierr); 452 PetscFunctionReturn(0); 453 } 454 455 /*@C 456 457 PCPatchSetComputeOperator - Set the callback used to compute patch matrices 458 459 Input Parameters: 460 + pc - The PC 461 . func - The callback 462 - ctx - The user context 463 464 Level: advanced 465 466 Note: 467 The callback has signature: 468 + usercomputeop(pc, point, mat, cellIS, n, u, ctx) 469 + pc - The PC 470 + point - The point 471 + mat - The patch matrix 472 + cellIS - An array of the cell numbers 473 + n - The size of g2l 474 + g2l - The global to local dof translation table 475 + ctx - The user context 476 and can assume that the matrix entries have been set to zero before the call. 477 478 .seealso: PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo() 479 @*/ 480 PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC, PetscInt, Mat, IS, PetscInt, const PetscInt *, void *), void *ctx) 481 { 482 PC_PATCH *patch = (PC_PATCH *) pc->data; 483 484 PetscFunctionBegin; 485 patch->usercomputeop = func; 486 patch->usercomputectx = ctx; 487 PetscFunctionReturn(0); 488 } 489 490 /* On entry, ht contains the topological entities whose dofs we are responsible for solving for; 491 on exit, cht contains all the topological entities we need to compute their residuals. 492 In full generality this should incorporate knowledge of the sparsity pattern of the matrix; 493 here we assume a standard FE sparsity pattern.*/ 494 /* TODO: Use DMPlexGetAdjacency() */ 495 static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht) 496 { 497 DM dm; 498 PetscHashIter hi; 499 PetscInt point; 500 PetscInt *star = NULL, *closure = NULL; 501 PetscInt ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci; 502 PetscErrorCode ierr; 503 504 PetscFunctionBegin; 505 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 506 ierr = PCPatchGetIgnoreDim(pc, &ignoredim);CHKERRQ(ierr); 507 if (ignoredim >= 0) {ierr = DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd);CHKERRQ(ierr);} 508 ierr = PetscHSetIClear(cht);CHKERRQ(ierr); 509 PetscHashIterBegin(ht, hi); 510 while (!PetscHashIterAtEnd(ht, hi)) { 511 512 PetscHashIterGetKey(ht, hi, point); 513 PetscHashIterNext(ht, hi); 514 515 /* Loop over all the cells that this point connects to */ 516 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 517 for (si = 0; si < starSize*2; si += 2) { 518 const PetscInt ownedpoint = star[si]; 519 /* TODO Check for point in cht before running through closure again */ 520 /* now loop over all entities in the closure of that cell */ 521 ierr = DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 522 for (ci = 0; ci < closureSize*2; ci += 2) { 523 const PetscInt seenpoint = closure[ci]; 524 if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue; 525 ierr = PetscHSetIAdd(cht, seenpoint);CHKERRQ(ierr); 526 } 527 } 528 } 529 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);CHKERRQ(ierr); 530 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_FALSE, NULL, &star);CHKERRQ(ierr); 531 PetscFunctionReturn(0); 532 } 533 534 static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off) 535 { 536 PetscErrorCode ierr; 537 538 PetscFunctionBegin; 539 if (combined) { 540 if (f < 0) { 541 if (dof) {ierr = PetscSectionGetDof(dofSection[0], p, dof);CHKERRQ(ierr);} 542 if (off) {ierr = PetscSectionGetOffset(dofSection[0], p, off);CHKERRQ(ierr);} 543 } else { 544 if (dof) {ierr = PetscSectionGetFieldDof(dofSection[0], p, f, dof);CHKERRQ(ierr);} 545 if (off) {ierr = PetscSectionGetFieldOffset(dofSection[0], p, f, off);CHKERRQ(ierr);} 546 } 547 } else { 548 if (f < 0) { 549 PC_PATCH *patch = (PC_PATCH *) pc->data; 550 PetscInt fdof, g; 551 552 if (dof) { 553 *dof = 0; 554 for (g = 0; g < patch->nsubspaces; ++g) { 555 ierr = PetscSectionGetDof(dofSection[g], p, &fdof);CHKERRQ(ierr); 556 *dof += fdof; 557 } 558 } 559 if (off) {ierr = PetscSectionGetOffset(dofSection[0], p, off);CHKERRQ(ierr);} 560 } else { 561 if (dof) {ierr = PetscSectionGetDof(dofSection[f], p, dof);CHKERRQ(ierr);} 562 if (off) {ierr = PetscSectionGetOffset(dofSection[f], p, off);CHKERRQ(ierr);} 563 } 564 } 565 PetscFunctionReturn(0); 566 } 567 568 /* Given a hash table with a set of topological entities (pts), compute the degrees of 569 freedom in global concatenated numbering on those entities. 570 For Vanka smoothing, this needs to do something special: ignore dofs of the 571 constraint subspace on entities that aren't the base entity we're building the patch 572 around. */ 573 static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI* subspaces_to_exclude) 574 { 575 PC_PATCH *patch = (PC_PATCH *) pc->data; 576 PetscHashIter hi; 577 PetscInt ldof, loff; 578 PetscInt k, p; 579 PetscErrorCode ierr; 580 581 PetscFunctionBegin; 582 ierr = PetscHSetIClear(dofs);CHKERRQ(ierr); 583 for (k = 0; k < patch->nsubspaces; ++k) { 584 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 585 PetscInt bs = patch->bs[k]; 586 PetscInt j, l; 587 588 if (subspaces_to_exclude != NULL) { 589 PetscBool should_exclude_k = PETSC_FALSE; 590 PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k); 591 if (should_exclude_k) { 592 /* only get this subspace dofs at the base entity, not any others */ 593 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff);CHKERRQ(ierr); 594 if (0 == ldof) continue; 595 for (j = loff; j < ldof + loff; ++j) { 596 for (l = 0; l < bs; ++l) { 597 PetscInt dof = bs*j + l + subspaceOffset; 598 ierr = PetscHSetIAdd(dofs, dof);CHKERRQ(ierr); 599 } 600 } 601 continue; /* skip the other dofs of this subspace */ 602 } 603 } 604 605 PetscHashIterBegin(pts, hi); 606 while (!PetscHashIterAtEnd(pts, hi)) { 607 PetscHashIterGetKey(pts, hi, p); 608 PetscHashIterNext(pts, hi); 609 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);CHKERRQ(ierr); 610 if (0 == ldof) continue; 611 for (j = loff; j < ldof + loff; ++j) { 612 for (l = 0; l < bs; ++l) { 613 PetscInt dof = bs*j + l + subspaceOffset; 614 ierr = PetscHSetIAdd(dofs, dof);CHKERRQ(ierr); 615 } 616 } 617 } 618 } 619 PetscFunctionReturn(0); 620 } 621 622 /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */ 623 static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C) 624 { 625 PetscHashIter hi; 626 PetscInt key; 627 PetscBool flg; 628 PetscErrorCode ierr; 629 630 PetscFunctionBegin; 631 ierr = PetscHSetIClear(C);CHKERRQ(ierr); 632 PetscHashIterBegin(B, hi); 633 while (!PetscHashIterAtEnd(B, hi)) { 634 PetscHashIterGetKey(B, hi, key); 635 PetscHashIterNext(B, hi); 636 ierr = PetscHSetIHas(A, key, &flg);CHKERRQ(ierr); 637 if (!flg) {ierr = PetscHSetIAdd(C, key);CHKERRQ(ierr);} 638 } 639 PetscFunctionReturn(0); 640 } 641 642 /* 643 * PCPatchCreateCellPatches - create patches. 644 * 645 * Input Parameters: 646 * + dm - The DMPlex object defining the mesh 647 * 648 * Output Parameters: 649 * + cellCounts - Section with counts of cells around each vertex 650 * . cells - IS of the cell point indices of cells in each patch 651 * . pointCounts - Section with counts of cells around each vertex 652 * - point - IS of the cell point indices of cells in each patch 653 */ 654 static PetscErrorCode PCPatchCreateCellPatches(PC pc) 655 { 656 PC_PATCH *patch = (PC_PATCH *) pc->data; 657 DMLabel ghost = NULL; 658 DM dm, plex; 659 PetscHSetI ht, cht; 660 PetscSection cellCounts, pointCounts; 661 PetscInt *cellsArray, *pointsArray; 662 PetscInt numCells, numPoints; 663 const PetscInt *leaves; 664 PetscInt nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, v; 665 PetscBool isFiredrake; 666 PetscErrorCode ierr; 667 668 PetscFunctionBegin; 669 /* Used to keep track of the cells in the patch. */ 670 ierr = PetscHSetICreate(&ht);CHKERRQ(ierr); 671 ierr = PetscHSetICreate(&cht);CHKERRQ(ierr); 672 673 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 674 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC\n"); 675 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 676 ierr = DMPlexGetChart(plex, &pStart, &pEnd);CHKERRQ(ierr); 677 ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr); 678 679 if (patch->user_patches) { 680 ierr = patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);CHKERRQ(ierr); 681 vStart = 0; vEnd = patch->npatch; 682 } else if (patch->codim < 0) { 683 if (patch->dim < 0) {ierr = DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);CHKERRQ(ierr);} 684 else {ierr = DMPlexGetDepthStratum(plex, patch->dim, &vStart, &vEnd);CHKERRQ(ierr);} 685 } else {ierr = DMPlexGetHeightStratum(plex, patch->codim, &vStart, &vEnd);CHKERRQ(ierr);} 686 patch->npatch = vEnd - vStart; 687 688 /* These labels mark the owned points. We only create patches around points that this process owns. */ 689 ierr = DMHasLabel(dm, "pyop2_ghost", &isFiredrake);CHKERRQ(ierr); 690 if (isFiredrake) { 691 ierr = DMGetLabel(dm, "pyop2_ghost", &ghost);CHKERRQ(ierr); 692 ierr = DMLabelCreateIndex(ghost, pStart, pEnd);CHKERRQ(ierr); 693 } else { 694 PetscSF sf; 695 696 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 697 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 698 nleaves = PetscMax(nleaves, 0); 699 } 700 701 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts);CHKERRQ(ierr); 702 ierr = PetscObjectSetName((PetscObject) patch->cellCounts, "Patch Cell Layout");CHKERRQ(ierr); 703 cellCounts = patch->cellCounts; 704 ierr = PetscSectionSetChart(cellCounts, vStart, vEnd);CHKERRQ(ierr); 705 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts);CHKERRQ(ierr); 706 ierr = PetscObjectSetName((PetscObject) patch->pointCounts, "Patch Point Layout");CHKERRQ(ierr); 707 pointCounts = patch->pointCounts; 708 ierr = PetscSectionSetChart(pointCounts, vStart, vEnd);CHKERRQ(ierr); 709 /* Count cells and points in the patch surrounding each entity */ 710 for (v = vStart; v < vEnd; ++v) { 711 PetscHashIter hi; 712 PetscInt chtSize, loc = -1; 713 PetscBool flg; 714 715 if (!patch->user_patches) { 716 if (ghost) {ierr = DMLabelHasPoint(ghost, v, &flg);CHKERRQ(ierr);} 717 else {ierr = PetscFindInt(v, nleaves, leaves, &loc);CHKERRQ(ierr); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;} 718 /* Not an owned entity, don't make a cell patch. */ 719 if (flg) continue; 720 } 721 722 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 723 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 724 ierr = PetscHSetIGetSize(cht, &chtSize);CHKERRQ(ierr); 725 /* empty patch, continue */ 726 if (chtSize == 0) continue; 727 728 /* safe because size(cht) > 0 from above */ 729 PetscHashIterBegin(cht, hi); 730 while (!PetscHashIterAtEnd(cht, hi)) { 731 PetscInt point, pdof; 732 733 PetscHashIterGetKey(cht, hi, point); 734 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 735 if (pdof) {ierr = PetscSectionAddDof(pointCounts, v, 1);CHKERRQ(ierr);} 736 if (point >= cStart && point < cEnd) {ierr = PetscSectionAddDof(cellCounts, v, 1);CHKERRQ(ierr);} 737 PetscHashIterNext(cht, hi); 738 } 739 } 740 if (isFiredrake) {ierr = DMLabelDestroyIndex(ghost);CHKERRQ(ierr);} 741 742 ierr = PetscSectionSetUp(cellCounts);CHKERRQ(ierr); 743 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 744 ierr = PetscMalloc1(numCells, &cellsArray);CHKERRQ(ierr); 745 ierr = PetscSectionSetUp(pointCounts);CHKERRQ(ierr); 746 ierr = PetscSectionGetStorageSize(pointCounts, &numPoints);CHKERRQ(ierr); 747 ierr = PetscMalloc1(numPoints, &pointsArray);CHKERRQ(ierr); 748 749 /* Now that we know how much space we need, run through again and actually remember the cells. */ 750 for (v = vStart; v < vEnd; v++ ) { 751 PetscHashIter hi; 752 PetscInt dof, off, cdof, coff, pdof, n = 0, cn = 0; 753 754 ierr = PetscSectionGetDof(pointCounts, v, &dof);CHKERRQ(ierr); 755 ierr = PetscSectionGetOffset(pointCounts, v, &off);CHKERRQ(ierr); 756 ierr = PetscSectionGetDof(cellCounts, v, &cdof);CHKERRQ(ierr); 757 ierr = PetscSectionGetOffset(cellCounts, v, &coff);CHKERRQ(ierr); 758 if (dof <= 0) continue; 759 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 760 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 761 PetscHashIterBegin(cht, hi); 762 while (!PetscHashIterAtEnd(cht, hi)) { 763 PetscInt point; 764 765 PetscHashIterGetKey(cht, hi, point); 766 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 767 if (pdof) {pointsArray[off + n++] = point;} 768 if (point >= cStart && point < cEnd) {cellsArray[coff + cn++] = point;} 769 PetscHashIterNext(cht, hi); 770 } 771 if (cn != cdof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of cells in patch %D is %D, but should be %D", v, cn, cdof); 772 if (n != dof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of points in patch %D is %D, but should be %D", v, n, dof); 773 } 774 ierr = PetscHSetIDestroy(&ht);CHKERRQ(ierr); 775 ierr = PetscHSetIDestroy(&cht);CHKERRQ(ierr); 776 ierr = DMDestroy(&plex);CHKERRQ(ierr); 777 778 ierr = ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells);CHKERRQ(ierr); 779 ierr = PetscObjectSetName((PetscObject) patch->cells, "Patch Cells");CHKERRQ(ierr); 780 if (patch->viewCells) { 781 ierr = ObjectView((PetscObject) patch->cellCounts, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 782 ierr = ObjectView((PetscObject) patch->cells, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 783 } 784 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points);CHKERRQ(ierr); 785 ierr = PetscObjectSetName((PetscObject) patch->points, "Patch Points");CHKERRQ(ierr); 786 if (patch->viewPoints) { 787 ierr = ObjectView((PetscObject) patch->pointCounts, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 788 ierr = ObjectView((PetscObject) patch->points, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 789 } 790 PetscFunctionReturn(0); 791 } 792 793 /* 794 * PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches 795 * 796 * Input Parameters: 797 * + dm - The DMPlex object defining the mesh 798 * . cellCounts - Section with counts of cells around each vertex 799 * . cells - IS of the cell point indices of cells in each patch 800 * . cellNumbering - Section mapping plex cell points to Firedrake cell indices. 801 * . nodesPerCell - number of nodes per cell. 802 * - cellNodeMap - map from cells to node indices (nodesPerCell * numCells) 803 * 804 * Output Parameters: 805 * + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering 806 * . gtolCounts - Section with counts of dofs per cell patch 807 * - gtol - IS mapping from global dofs to local dofs for each patch. 808 */ 809 static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc) 810 { 811 PC_PATCH *patch = (PC_PATCH *) pc->data; 812 PetscSection cellCounts = patch->cellCounts; 813 PetscSection pointCounts = patch->pointCounts; 814 PetscSection gtolCounts; 815 IS cells = patch->cells; 816 IS points = patch->points; 817 PetscSection cellNumbering = patch->cellNumbering; 818 PetscInt Nf = patch->nsubspaces; 819 PetscInt numCells, numPoints; 820 PetscInt numDofs; 821 PetscInt numGlobalDofs; 822 PetscInt totalDofsPerCell = patch->totalDofsPerCell; 823 PetscInt vStart, vEnd, v; 824 const PetscInt *cellsArray, *pointsArray; 825 PetscInt *newCellsArray = NULL; 826 PetscInt *dofsArray = NULL; 827 PetscInt *offsArray = NULL; 828 PetscInt *asmArray = NULL; 829 PetscInt *globalDofsArray = NULL; 830 PetscInt globalIndex = 0; 831 PetscInt key = 0; 832 PetscInt asmKey = 0; 833 DM dm = NULL; 834 const PetscInt *bcNodes = NULL; 835 PetscHMapI ht; 836 PetscHSetI globalBcs; 837 PetscInt numBcs; 838 PetscHSetI ownedpts, seenpts, owneddofs, seendofs, artificialbcs; 839 PetscInt pStart, pEnd, p, i; 840 PetscErrorCode ierr; 841 842 PetscFunctionBegin; 843 844 ierr = PCGetDM(pc, &dm); CHKERRQ(ierr); 845 /* dofcounts section is cellcounts section * dofPerCell */ 846 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 847 ierr = PetscSectionGetStorageSize(patch->pointCounts, &numPoints);CHKERRQ(ierr); 848 numDofs = numCells * totalDofsPerCell; 849 ierr = PetscMalloc1(numDofs, &dofsArray);CHKERRQ(ierr); 850 ierr = PetscMalloc1(numPoints*Nf, &offsArray);CHKERRQ(ierr); 851 ierr = PetscMalloc1(numDofs, &asmArray);CHKERRQ(ierr); 852 ierr = PetscMalloc1(numCells, &newCellsArray);CHKERRQ(ierr); 853 ierr = PetscSectionGetChart(cellCounts, &vStart, &vEnd);CHKERRQ(ierr); 854 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);CHKERRQ(ierr); 855 gtolCounts = patch->gtolCounts; 856 ierr = PetscSectionSetChart(gtolCounts, vStart, vEnd);CHKERRQ(ierr); 857 ierr = PetscObjectSetName((PetscObject) patch->gtolCounts, "Patch Global Index Section");CHKERRQ(ierr); 858 859 /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet 860 conditions */ 861 ierr = PetscHSetICreate(&globalBcs);CHKERRQ(ierr); 862 ierr = ISGetIndices(patch->ghostBcNodes, &bcNodes); CHKERRQ(ierr); 863 ierr = ISGetSize(patch->ghostBcNodes, &numBcs); CHKERRQ(ierr); 864 for (i = 0; i < numBcs; ++i) { 865 ierr = PetscHSetIAdd(globalBcs, bcNodes[i]);CHKERRQ(ierr); /* these are already in concatenated numbering */ 866 } 867 ierr = ISRestoreIndices(patch->ghostBcNodes, &bcNodes); CHKERRQ(ierr); 868 ierr = ISDestroy(&patch->ghostBcNodes); CHKERRQ(ierr); /* memory optimisation */ 869 870 /* Hash tables for artificial BC construction */ 871 ierr = PetscHSetICreate(&ownedpts);CHKERRQ(ierr); 872 ierr = PetscHSetICreate(&seenpts);CHKERRQ(ierr); 873 ierr = PetscHSetICreate(&owneddofs);CHKERRQ(ierr); 874 ierr = PetscHSetICreate(&seendofs);CHKERRQ(ierr); 875 ierr = PetscHSetICreate(&artificialbcs);CHKERRQ(ierr); 876 877 ierr = ISGetIndices(cells, &cellsArray);CHKERRQ(ierr); 878 ierr = ISGetIndices(points, &pointsArray);CHKERRQ(ierr); 879 ierr = PetscHMapICreate(&ht);CHKERRQ(ierr); 880 for (v = vStart; v < vEnd; ++v) { 881 PetscInt localIndex = 0; 882 PetscInt dof, off, i, j, k, l; 883 884 ierr = PetscHMapIClear(ht);CHKERRQ(ierr); 885 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 886 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 887 if (dof <= 0) continue; 888 889 /* Calculate the global numbers of the artificial BC dofs here first */ 890 ierr = patch->patchconstructop((void*)patch, dm, v, ownedpts); CHKERRQ(ierr); 891 ierr = PCPatchCompleteCellPatch(pc, ownedpts, seenpts); CHKERRQ(ierr); 892 ierr = PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude); CHKERRQ(ierr); 893 ierr = PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL); CHKERRQ(ierr); 894 ierr = PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs); CHKERRQ(ierr); 895 if (patch->viewPatches) { 896 PetscHSetI globalbcdofs; 897 PetscHashIter hi; 898 MPI_Comm comm = PetscObjectComm((PetscObject)pc); 899 900 ierr = PetscHSetICreate(&globalbcdofs);CHKERRQ(ierr); 901 ierr = PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v); CHKERRQ(ierr); 902 PetscHashIterBegin(owneddofs, hi); 903 while (!PetscHashIterAtEnd(owneddofs, hi)) { 904 PetscInt globalDof; 905 906 PetscHashIterGetKey(owneddofs, hi, globalDof); 907 PetscHashIterNext(owneddofs, hi); 908 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 909 } 910 ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr); 911 ierr = PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v); CHKERRQ(ierr); 912 PetscHashIterBegin(seendofs, hi); 913 while (!PetscHashIterAtEnd(seendofs, hi)) { 914 PetscInt globalDof; 915 PetscBool flg; 916 917 PetscHashIterGetKey(seendofs, hi, globalDof); 918 PetscHashIterNext(seendofs, hi); 919 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 920 921 ierr = PetscHSetIHas(globalBcs, globalDof, &flg);CHKERRQ(ierr); 922 if (flg) {ierr = PetscHSetIAdd(globalbcdofs, globalDof);CHKERRQ(ierr);} 923 } 924 ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr); 925 ierr = PetscSynchronizedPrintf(comm, "Patch %d: global BCs:\n", v); CHKERRQ(ierr); 926 ierr = PetscHSetIGetSize(globalbcdofs, &numBcs);CHKERRQ(ierr); 927 if (numBcs > 0) { 928 PetscHashIterBegin(globalbcdofs, hi); 929 while (!PetscHashIterAtEnd(globalbcdofs, hi)) { 930 PetscInt globalDof; 931 PetscHashIterGetKey(globalbcdofs, hi, globalDof); 932 PetscHashIterNext(globalbcdofs, hi); 933 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof);CHKERRQ(ierr); 934 } 935 } 936 ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr); 937 ierr = PetscSynchronizedPrintf(comm, "Patch %d: artificial BCs:\n", v);CHKERRQ(ierr); 938 ierr = PetscHSetIGetSize(artificialbcs, &numBcs);CHKERRQ(ierr); 939 if (numBcs > 0) { 940 PetscHashIterBegin(artificialbcs, hi); 941 while (!PetscHashIterAtEnd(artificialbcs, hi)) { 942 PetscInt globalDof; 943 PetscHashIterGetKey(artificialbcs, hi, globalDof); 944 PetscHashIterNext(artificialbcs, hi); 945 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 946 } 947 } 948 ierr = PetscSynchronizedPrintf(comm, "\n\n"); CHKERRQ(ierr); 949 ierr = PetscHSetIDestroy(&globalbcdofs);CHKERRQ(ierr); 950 } 951 for (k = 0; k < patch->nsubspaces; ++k) { 952 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 953 PetscInt nodesPerCell = patch->nodesPerCell[k]; 954 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 955 PetscInt bs = patch->bs[k]; 956 957 for (i = off; i < off + dof; ++i) { 958 /* Walk over the cells in this patch. */ 959 const PetscInt c = cellsArray[i]; 960 PetscInt cell = c; 961 962 /* TODO Change this to an IS */ 963 if (cellNumbering) { 964 ierr = PetscSectionGetDof(cellNumbering, c, &cell);CHKERRQ(ierr); 965 if (cell <= 0) SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %D doesn't appear in cell numbering map", c); 966 ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr); 967 } 968 newCellsArray[i] = cell; 969 for (j = 0; j < nodesPerCell; ++j) { 970 /* For each global dof, map it into contiguous local storage. */ 971 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + subspaceOffset; 972 /* finally, loop over block size */ 973 for (l = 0; l < bs; ++l) { 974 PetscInt localDof; 975 PetscBool isGlobalBcDof, isArtificialBcDof; 976 977 /* first, check if this is either a globally enforced or locally enforced BC dof */ 978 ierr = PetscHSetIHas(globalBcs, globalDof + l, &isGlobalBcDof);CHKERRQ(ierr); 979 ierr = PetscHSetIHas(artificialbcs, globalDof + l, &isArtificialBcDof);CHKERRQ(ierr); 980 981 /* if it's either, don't ever give it a local dof number */ 982 if (isGlobalBcDof || isArtificialBcDof) { 983 dofsArray[globalIndex++] = -1; /* don't use this in assembly in this patch */ 984 } else { 985 ierr = PetscHMapIGet(ht, globalDof + l, &localDof);CHKERRQ(ierr); 986 if (localDof == -1) { 987 localDof = localIndex++; 988 ierr = PetscHMapISet(ht, globalDof + l, localDof);CHKERRQ(ierr); 989 } 990 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 991 /* And store. */ 992 dofsArray[globalIndex++] = localDof; 993 } 994 } 995 } 996 } 997 } 998 /* How many local dofs in this patch? */ 999 ierr = PetscHMapIGetSize(ht, &dof);CHKERRQ(ierr); 1000 ierr = PetscSectionSetDof(gtolCounts, v, dof);CHKERRQ(ierr); 1001 } 1002 if (globalIndex != numDofs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%d) doesn't match found number (%d)", numDofs, globalIndex); 1003 ierr = PetscSectionSetUp(gtolCounts);CHKERRQ(ierr); 1004 ierr = PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);CHKERRQ(ierr); 1005 ierr = PetscMalloc1(numGlobalDofs, &globalDofsArray);CHKERRQ(ierr); 1006 1007 /* Now populate the global to local map. This could be merged into the above loop if we were willing to deal with reallocs. */ 1008 for (v = vStart; v < vEnd; ++v) { 1009 PetscHashIter hi; 1010 PetscInt dof, off, Np, ooff, i, j, k, l; 1011 1012 ierr = PetscHMapIClear(ht);CHKERRQ(ierr); 1013 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 1014 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 1015 ierr = PetscSectionGetDof(pointCounts, v, &Np);CHKERRQ(ierr); 1016 ierr = PetscSectionGetOffset(pointCounts, v, &ooff);CHKERRQ(ierr); 1017 if (dof <= 0) continue; 1018 1019 for (k = 0; k < patch->nsubspaces; ++k) { 1020 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1021 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1022 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1023 PetscInt bs = patch->bs[k]; 1024 PetscInt goff; 1025 1026 for (i = off; i < off + dof; ++i) { 1027 /* Reconstruct mapping of global-to-local on this patch. */ 1028 const PetscInt c = cellsArray[i]; 1029 PetscInt cell = c; 1030 1031 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 1032 for (j = 0; j < nodesPerCell; ++j) { 1033 for (l = 0; l < bs; ++l) { 1034 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1035 const PetscInt localDof = dofsArray[key++]; 1036 1037 if (localDof >= 0) {ierr = PetscHMapISet(ht, globalDof, localDof);CHKERRQ(ierr);} 1038 } 1039 } 1040 } 1041 1042 /* Shove it in the output data structure. */ 1043 ierr = PetscSectionGetOffset(gtolCounts, v, &goff);CHKERRQ(ierr); 1044 PetscHashIterBegin(ht, hi); 1045 while (!PetscHashIterAtEnd(ht, hi)) { 1046 PetscInt globalDof, localDof; 1047 1048 PetscHashIterGetKey(ht, hi, globalDof); 1049 PetscHashIterGetVal(ht, hi, localDof); 1050 if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof; 1051 PetscHashIterNext(ht, hi); 1052 } 1053 1054 for (p = 0; p < Np; ++p) { 1055 const PetscInt point = pointsArray[ooff + p]; 1056 PetscInt globalDof, localDof; 1057 1058 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);CHKERRQ(ierr); 1059 ierr = PetscHMapIGet(ht, globalDof, &localDof);CHKERRQ(ierr); 1060 offsArray[(ooff + p)*Nf + k] = localDof; 1061 } 1062 } 1063 1064 ierr = PetscHSetIDestroy(&globalBcs);CHKERRQ(ierr); 1065 ierr = PetscHSetIDestroy(&ownedpts);CHKERRQ(ierr); 1066 ierr = PetscHSetIDestroy(&seenpts);CHKERRQ(ierr); 1067 ierr = PetscHSetIDestroy(&owneddofs);CHKERRQ(ierr); 1068 ierr = PetscHSetIDestroy(&seendofs);CHKERRQ(ierr); 1069 ierr = PetscHSetIDestroy(&artificialbcs);CHKERRQ(ierr); 1070 1071 /* At this point, we have a hash table ht built that maps globalDof -> localDof. 1072 We need to create the dof table laid out cellwise first, then by subspace, 1073 as the assembler assembles cell-wise and we need to stuff the different 1074 contributions of the different function spaces to the right places. So we loop 1075 over cells, then over subspaces. */ 1076 if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */ 1077 for (i = off; i < off + dof; ++i) { 1078 const PetscInt c = cellsArray[i]; 1079 PetscInt cell = c; 1080 1081 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 1082 for (k = 0; k < patch->nsubspaces; ++k) { 1083 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1084 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1085 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1086 PetscInt bs = patch->bs[k]; 1087 1088 for (j = 0; j < nodesPerCell; ++j) { 1089 for (l = 0; l < bs; ++l) { 1090 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1091 PetscInt localDof; 1092 1093 ierr = PetscHMapIGet(ht, globalDof, &localDof);CHKERRQ(ierr); 1094 /* If it's not in the hash table, i.e. is a BC dof, 1095 then the PetscHSetIMap above gives -1, which matches 1096 exactly the convention for PETSc's matrix assembly to 1097 ignore the dof. So we don't need to do anything here */ 1098 asmArray[asmKey++] = localDof; 1099 } 1100 } 1101 } 1102 } 1103 } 1104 } 1105 if (1 == patch->nsubspaces) {ierr = PetscMemcpy(asmArray, dofsArray, numDofs * sizeof(PetscInt));CHKERRQ(ierr);} 1106 1107 ierr = PetscHMapIDestroy(&ht);CHKERRQ(ierr); 1108 ierr = ISRestoreIndices(cells, &cellsArray);CHKERRQ(ierr); 1109 ierr = ISRestoreIndices(points, &pointsArray);CHKERRQ(ierr); 1110 ierr = PetscFree(dofsArray);CHKERRQ(ierr); 1111 /* Create placeholder section for map from points to patch dofs */ 1112 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);CHKERRQ(ierr); 1113 ierr = PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);CHKERRQ(ierr); 1114 if (patch->combined) { 1115 PetscInt numFields; 1116 ierr = PetscSectionGetNumFields(patch->dofSection[0], &numFields);CHKERRQ(ierr); 1117 if (numFields != patch->nsubspaces) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Mismatch between number of section fields %D and number of subspaces %D", numFields, patch->nsubspaces); 1118 ierr = PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);CHKERRQ(ierr); 1119 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1120 for (p = pStart; p < pEnd; ++p) { 1121 PetscInt dof, fdof, f; 1122 1123 ierr = PetscSectionGetDof(patch->dofSection[0], p, &dof);CHKERRQ(ierr); 1124 ierr = PetscSectionSetDof(patch->patchSection, p, dof);CHKERRQ(ierr); 1125 for (f = 0; f < patch->nsubspaces; ++f) { 1126 ierr = PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof);CHKERRQ(ierr); 1127 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1128 } 1129 } 1130 } else { 1131 PetscInt pStartf, pEndf, f; 1132 pStart = PETSC_MAX_INT; 1133 pEnd = PETSC_MIN_INT; 1134 for (f = 0; f < patch->nsubspaces; ++f) { 1135 ierr = PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);CHKERRQ(ierr); 1136 pStart = PetscMin(pStart, pStartf); 1137 pEnd = PetscMax(pEnd, pEndf); 1138 } 1139 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1140 for (f = 0; f < patch->nsubspaces; ++f) { 1141 ierr = PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);CHKERRQ(ierr); 1142 for (p = pStartf; p < pEndf; ++p) { 1143 PetscInt fdof; 1144 ierr = PetscSectionGetDof(patch->dofSection[f], p, &fdof);CHKERRQ(ierr); 1145 ierr = PetscSectionAddDof(patch->patchSection, p, fdof);CHKERRQ(ierr); 1146 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1147 } 1148 } 1149 } 1150 ierr = PetscSectionSetUp(patch->patchSection);CHKERRQ(ierr); 1151 ierr = PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE);CHKERRQ(ierr); 1152 /* Replace cell indices with firedrake-numbered ones. */ 1153 ierr = ISGeneralSetIndices(cells, numCells, (const PetscInt *) newCellsArray, PETSC_OWN_POINTER);CHKERRQ(ierr); 1154 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol);CHKERRQ(ierr); 1155 ierr = PetscObjectSetName((PetscObject) patch->gtol, "Global Indices");CHKERRQ(ierr); 1156 ierr = PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject) pc, "-pc_patch_g2l_view");CHKERRQ(ierr); 1157 ierr = ISViewFromOptions(patch->gtol, (PetscObject) pc, "-pc_patch_g2l_view");CHKERRQ(ierr); 1158 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs);CHKERRQ(ierr); 1159 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArray, PETSC_OWN_POINTER, &patch->offs);CHKERRQ(ierr); 1160 PetscFunctionReturn(0); 1161 } 1162 1163 static PetscErrorCode PCPatchZeroFillMatrix_Private(Mat mat, const PetscInt ncell, const PetscInt ndof, const PetscInt *dof) 1164 { 1165 PetscScalar *values = NULL; 1166 PetscInt rows, c, i; 1167 PetscErrorCode ierr; 1168 1169 PetscFunctionBegin; 1170 ierr = PetscCalloc1(ndof*ndof, &values);CHKERRQ(ierr); 1171 for (c = 0; c < ncell; ++c) { 1172 const PetscInt *idx = &dof[ndof*c]; 1173 ierr = MatSetValues(mat, ndof, idx, ndof, idx, values, INSERT_VALUES);CHKERRQ(ierr); 1174 } 1175 ierr = MatGetLocalSize(mat, &rows, NULL);CHKERRQ(ierr); 1176 for (i = 0; i < rows; ++i) { 1177 ierr = MatSetValues(mat, 1, &i, 1, &i, values, INSERT_VALUES);CHKERRQ(ierr); 1178 } 1179 ierr = MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1180 ierr = MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1181 ierr = PetscFree(values);CHKERRQ(ierr); 1182 PetscFunctionReturn(0); 1183 } 1184 1185 static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat) 1186 { 1187 PC_PATCH *patch = (PC_PATCH *) pc->data; 1188 Vec x, y; 1189 PetscBool flg; 1190 PetscInt csize, rsize; 1191 const char *prefix = NULL; 1192 PetscErrorCode ierr; 1193 1194 PetscFunctionBegin; 1195 x = patch->patchX[point]; 1196 y = patch->patchY[point]; 1197 ierr = VecGetSize(x, &csize);CHKERRQ(ierr); 1198 ierr = VecGetSize(y, &rsize);CHKERRQ(ierr); 1199 ierr = MatCreate(PETSC_COMM_SELF, mat);CHKERRQ(ierr); 1200 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1201 ierr = MatSetOptionsPrefix(*mat, prefix);CHKERRQ(ierr); 1202 ierr = MatAppendOptionsPrefix(*mat, "pc_patch_sub_");CHKERRQ(ierr); 1203 if (patch->sub_mat_type) {ierr = MatSetType(*mat, patch->sub_mat_type);CHKERRQ(ierr);} 1204 else if (!patch->sub_mat_type) {ierr = MatSetType(*mat, MATDENSE);CHKERRQ(ierr);} 1205 ierr = MatSetSizes(*mat, rsize, csize, rsize, csize);CHKERRQ(ierr); 1206 ierr = PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);CHKERRQ(ierr); 1207 if (!flg) {ierr = PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);CHKERRQ(ierr);} 1208 /* Sparse patch matrices */ 1209 if (!flg) { 1210 PetscBT bt; 1211 PetscInt *dnnz = NULL; 1212 const PetscInt *dofsArray = NULL; 1213 PetscInt pStart, pEnd, ncell, offset, c, i, j; 1214 1215 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1216 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1217 point += pStart; 1218 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1219 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1220 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1221 ierr = PetscCalloc1(rsize, &dnnz);CHKERRQ(ierr); 1222 ierr = PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1223 /* XXX: This uses N^2 bits to store the sparsity pattern on a 1224 * patch. This is probably OK if the patches are not too big, 1225 * but could use quite a bit of memory for planes in 3D. 1226 * Should we switch based on the value of rsize to a 1227 * hash-table (slower, but more memory efficient) approach? */ 1228 ierr = PetscBTCreate(rsize*rsize, &bt);CHKERRQ(ierr); 1229 for (c = 0; c < ncell; ++c) { 1230 const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell; 1231 for (i = 0; i < patch->totalDofsPerCell; ++i) { 1232 const PetscInt row = idx[i]; 1233 if (row < 0) continue; 1234 for (j = 0; j < patch->totalDofsPerCell; ++j) { 1235 const PetscInt col = idx[j]; 1236 const PetscInt key = row*rsize + col; 1237 if (col < 0) continue; 1238 if (!PetscBTLookupSet(bt, key)) ++dnnz[row]; 1239 } 1240 } 1241 } 1242 ierr = PetscBTDestroy(&bt);CHKERRQ(ierr); 1243 ierr = MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);CHKERRQ(ierr); 1244 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1245 ierr = PCPatchZeroFillMatrix_Private(*mat, ncell, patch->totalDofsPerCell, &dofsArray[offset*patch->totalDofsPerCell]);CHKERRQ(ierr); 1246 ierr = PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1247 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1248 } 1249 ierr = MatSetUp(*mat);CHKERRQ(ierr); 1250 PetscFunctionReturn(0); 1251 } 1252 1253 static PetscErrorCode PCPatchComputeOperator_DMPlex_Private(PC pc, PetscInt patchNum, Mat J, IS cellIS, PetscInt n, const PetscInt *l2p, void *ctx) 1254 { 1255 PC_PATCH *patch = (PC_PATCH *) pc->data; 1256 DM dm; 1257 PetscSection s; 1258 const PetscInt *parray, *oarray; 1259 PetscInt Nf = patch->nsubspaces, Np, poff, p, f; 1260 PetscErrorCode ierr; 1261 1262 PetscFunctionBegin; 1263 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1264 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1265 /* Set offset into patch */ 1266 ierr = PetscSectionGetDof(patch->pointCounts, patchNum, &Np);CHKERRQ(ierr); 1267 ierr = PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);CHKERRQ(ierr); 1268 ierr = ISGetIndices(patch->points, &parray);CHKERRQ(ierr); 1269 ierr = ISGetIndices(patch->offs, &oarray);CHKERRQ(ierr); 1270 for (f = 0; f < Nf; ++f) { 1271 for (p = 0; p < Np; ++p) { 1272 const PetscInt point = parray[poff+p]; 1273 PetscInt dof; 1274 1275 ierr = PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);CHKERRQ(ierr); 1276 ierr = PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr); 1277 if (patch->nsubspaces == 1) {ierr = PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr);} 1278 else {ierr = PetscSectionSetOffset(patch->patchSection, point, -1);CHKERRQ(ierr);} 1279 } 1280 } 1281 ierr = ISRestoreIndices(patch->points, &parray);CHKERRQ(ierr); 1282 ierr = ISRestoreIndices(patch->offs, &oarray);CHKERRQ(ierr); 1283 if (patch->viewSection) {ierr = ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);CHKERRQ(ierr);} 1284 /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */ 1285 ierr = DMPlexComputeJacobian_Patch_Internal(pc->dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, NULL, NULL, J, J, ctx);CHKERRQ(ierr); 1286 PetscFunctionReturn(0); 1287 } 1288 1289 static PetscErrorCode PCPatchComputeOperator_Private(PC pc, Mat mat, PetscInt point) 1290 { 1291 PC_PATCH *patch = (PC_PATCH *) pc->data; 1292 const PetscInt *dofsArray; 1293 const PetscInt *cellsArray; 1294 PetscInt ncell, offset, pStart, pEnd; 1295 PetscErrorCode ierr; 1296 1297 PetscFunctionBegin; 1298 ierr = PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1299 if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n"); 1300 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1301 ierr = ISGetIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1302 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1303 1304 point += pStart; 1305 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1306 1307 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1308 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1309 if (ncell <= 0) { 1310 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1311 PetscFunctionReturn(0); 1312 } 1313 PetscStackPush("PCPatch user callback"); 1314 /* Cannot reuse the same IS because the geometry info is being cached in it */ 1315 ierr = ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);CHKERRQ(ierr); 1316 ierr = patch->usercomputeop(pc, point, mat, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, patch->usercomputectx);CHKERRQ(ierr); 1317 PetscStackPop; 1318 ierr = ISDestroy(&patch->cellIS);CHKERRQ(ierr); 1319 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1320 ierr = ISRestoreIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1321 if (patch->viewMatrix) { 1322 char name[PETSC_MAX_PATH_LEN]; 1323 1324 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);CHKERRQ(ierr); 1325 ierr = PetscObjectSetName((PetscObject) mat, name);CHKERRQ(ierr); 1326 ierr = ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);CHKERRQ(ierr); 1327 } 1328 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1329 PetscFunctionReturn(0); 1330 } 1331 1332 static PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat) 1333 { 1334 PC_PATCH *patch = (PC_PATCH *) pc->data; 1335 const PetscScalar *xArray = NULL; 1336 PetscScalar *yArray = NULL; 1337 const PetscInt *gtolArray = NULL; 1338 PetscInt dof, offset, lidx; 1339 PetscErrorCode ierr; 1340 1341 PetscFunctionBeginHot; 1342 ierr = PetscLogEventBegin(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 1343 ierr = VecGetArrayRead(x, &xArray);CHKERRQ(ierr); 1344 ierr = VecGetArray(y, &yArray);CHKERRQ(ierr); 1345 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 1346 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 1347 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1348 if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n"); 1349 if (mode == ADD_VALUES && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n"); 1350 for (lidx = 0; lidx < dof; ++lidx) { 1351 const PetscInt gidx = gtolArray[offset+lidx]; 1352 1353 if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */ 1354 else yArray[gidx] += xArray[lidx]; /* Reverse */ 1355 } 1356 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1357 ierr = VecRestoreArrayRead(x, &xArray);CHKERRQ(ierr); 1358 ierr = VecRestoreArray(y, &yArray);CHKERRQ(ierr); 1359 ierr = PetscLogEventEnd(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 1360 PetscFunctionReturn(0); 1361 } 1362 1363 static PetscErrorCode PCSetUp_PATCH(PC pc) 1364 { 1365 PC_PATCH *patch = (PC_PATCH *) pc->data; 1366 PetscInt i; 1367 const char *prefix; 1368 PetscErrorCode ierr; 1369 1370 PetscFunctionBegin; 1371 if (!pc->setupcalled) { 1372 PetscInt pStart, pEnd, p; 1373 PetscInt localSize; 1374 1375 ierr = PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 1376 1377 if (!patch->nsubspaces) { 1378 DM dm; 1379 PetscDS prob; 1380 PetscSection s; 1381 PetscInt cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, totNb = 0, **cellDofs; 1382 1383 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1384 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()"); 1385 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1386 ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr); 1387 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr); 1388 for (p = pStart; p < pEnd; ++p) { 1389 PetscInt cdof; 1390 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 1391 numGlobalBcs += cdof; 1392 } 1393 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1394 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1395 ierr = PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);CHKERRQ(ierr); 1396 for (f = 0; f < Nf; ++f) { 1397 PetscFE fe; 1398 PetscDualSpace sp; 1399 PetscInt cdoff = 0; 1400 1401 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1402 /* ierr = PetscFEGetNumComponents(fe, &Nc[f]);CHKERRQ(ierr); */ 1403 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 1404 ierr = PetscDualSpaceGetDimension(sp, &Nb[f]);CHKERRQ(ierr); 1405 totNb += Nb[f]; 1406 1407 ierr = PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);CHKERRQ(ierr); 1408 for (c = cStart; c < cEnd; ++c) { 1409 PetscInt *closure = NULL; 1410 PetscInt clSize = 0, cl; 1411 1412 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1413 for (cl = 0; cl < clSize*2; cl += 2) { 1414 const PetscInt p = closure[cl]; 1415 PetscInt fdof, d, foff; 1416 1417 ierr = PetscSectionGetFieldDof(s, p, f, &fdof);CHKERRQ(ierr); 1418 ierr = PetscSectionGetFieldOffset(s, p, f, &foff);CHKERRQ(ierr); 1419 for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d; 1420 } 1421 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1422 } 1423 if (cdoff != (cEnd-cStart)*Nb[f]) SETERRQ4(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_SIZ, "Total number of cellDofs %D for field %D should be Nc (%D) * cellDof (%D)", cdoff, f, cEnd-cStart, Nb[f]); 1424 } 1425 numGlobalBcs = 0; 1426 for (p = pStart; p < pEnd; ++p) { 1427 const PetscInt *ind; 1428 PetscInt off, cdof, d; 1429 1430 ierr = PetscSectionGetOffset(s, p, &off);CHKERRQ(ierr); 1431 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 1432 ierr = PetscSectionGetConstraintIndices(s, p, &ind);CHKERRQ(ierr); 1433 for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d]; 1434 } 1435 1436 ierr = PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);CHKERRQ(ierr); 1437 for (f = 0; f < Nf; ++f) { 1438 ierr = PetscFree(cellDofs[f]);CHKERRQ(ierr); 1439 } 1440 ierr = PetscFree3(Nb, cellDofs, globalBcs);CHKERRQ(ierr); 1441 ierr = PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);CHKERRQ(ierr); 1442 } 1443 1444 localSize = patch->subspaceOffsets[patch->nsubspaces]; 1445 ierr = VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localX);CHKERRQ(ierr); 1446 ierr = VecSetUp(patch->localX);CHKERRQ(ierr); 1447 ierr = VecDuplicate(patch->localX, &patch->localY);CHKERRQ(ierr); 1448 ierr = PCPatchCreateCellPatches(pc);CHKERRQ(ierr); 1449 ierr = PCPatchCreateCellPatchDiscretisationInfo(pc);CHKERRQ(ierr); 1450 1451 /* OK, now build the work vectors */ 1452 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd);CHKERRQ(ierr); 1453 ierr = PetscMalloc1(patch->npatch, &patch->patchX);CHKERRQ(ierr); 1454 ierr = PetscMalloc1(patch->npatch, &patch->patchY);CHKERRQ(ierr); 1455 for (p = pStart; p < pEnd; ++p) { 1456 PetscInt dof; 1457 1458 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 1459 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchX[p-pStart]);CHKERRQ(ierr); 1460 ierr = VecSetUp(patch->patchX[p-pStart]);CHKERRQ(ierr); 1461 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchY[p-pStart]);CHKERRQ(ierr); 1462 ierr = VecSetUp(patch->patchY[p-pStart]);CHKERRQ(ierr); 1463 } 1464 ierr = PetscMalloc1(patch->npatch, &patch->ksp);CHKERRQ(ierr); 1465 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1466 for (i = 0; i < patch->npatch; ++i) { 1467 PC subpc; 1468 1469 ierr = KSPCreate(PETSC_COMM_SELF, &patch->ksp[i]);CHKERRQ(ierr); 1470 ierr = KSPSetOptionsPrefix(patch->ksp[i], prefix);CHKERRQ(ierr); 1471 ierr = KSPAppendOptionsPrefix(patch->ksp[i], "sub_");CHKERRQ(ierr); 1472 ierr = PetscObjectIncrementTabLevel((PetscObject) patch->ksp[i], (PetscObject) pc, 1);CHKERRQ(ierr); 1473 ierr = KSPGetPC(patch->ksp[i], &subpc);CHKERRQ(ierr); 1474 ierr = PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);CHKERRQ(ierr); 1475 ierr = PetscLogObjectParent((PetscObject) pc, (PetscObject) patch->ksp[i]);CHKERRQ(ierr); 1476 } 1477 if (patch->save_operators) { 1478 ierr = PetscMalloc1(patch->npatch, &patch->mat);CHKERRQ(ierr); 1479 for (i = 0; i < patch->npatch; ++i) { 1480 ierr = PCPatchCreateMatrix_Private(pc, i, &patch->mat[i]);CHKERRQ(ierr); 1481 } 1482 } 1483 ierr = PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 1484 1485 /* If desired, calculate weights for dof multiplicity */ 1486 if (patch->partition_of_unity) { 1487 ierr = VecDuplicate(patch->localX, &patch->dof_weights);CHKERRQ(ierr); 1488 for (i = 0; i < patch->npatch; ++i) { 1489 PetscInt dof; 1490 1491 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);CHKERRQ(ierr); 1492 if (dof <= 0) continue; 1493 ierr = VecSet(patch->patchX[i], 1.0);CHKERRQ(ierr); 1494 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchX[i], patch->dof_weights, ADD_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 1495 } 1496 ierr = VecReciprocal(patch->dof_weights);CHKERRQ(ierr); 1497 } 1498 } 1499 if (patch->save_operators) { 1500 for (i = 0; i < patch->npatch; ++i) { 1501 ierr = MatZeroEntries(patch->mat[i]);CHKERRQ(ierr); 1502 ierr = PCPatchComputeOperator_Private(pc, patch->mat[i], i);CHKERRQ(ierr); 1503 ierr = KSPSetOperators(patch->ksp[i], patch->mat[i], patch->mat[i]);CHKERRQ(ierr); 1504 } 1505 } 1506 if (!pc->setupcalled && patch->optionsSet) for (i = 0; i < patch->npatch; ++i) {ierr = KSPSetFromOptions(patch->ksp[i]);CHKERRQ(ierr);} 1507 PetscFunctionReturn(0); 1508 } 1509 1510 static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y) 1511 { 1512 PC_PATCH *patch = (PC_PATCH *) pc->data; 1513 const PetscScalar *globalX = NULL; 1514 PetscScalar *localX = NULL; 1515 PetscScalar *globalY = NULL; 1516 const PetscInt *bcNodes = NULL; 1517 PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1; 1518 PetscInt start[2] = {0, 0}; 1519 PetscInt end[2] = {-1, -1}; 1520 const PetscInt inc[2] = {1, -1}; 1521 const PetscScalar *localY; 1522 const PetscInt *iterationSet; 1523 PetscInt pStart, numBcs, n, sweep, bc, j; 1524 PetscErrorCode ierr; 1525 1526 PetscFunctionBegin; 1527 ierr = PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 1528 ierr = PetscOptionsPushGetViewerOff(PETSC_TRUE);CHKERRQ(ierr); 1529 end[0] = patch->npatch; 1530 start[1] = patch->npatch-1; 1531 if (patch->user_patches) { 1532 ierr = ISGetLocalSize(patch->iterationSet, &end[0]);CHKERRQ(ierr); 1533 start[1] = end[0] - 1; 1534 ierr = ISGetIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr); 1535 } 1536 /* Scatter from global space into overlapped local spaces */ 1537 ierr = VecGetArrayRead(x, &globalX);CHKERRQ(ierr); 1538 ierr = VecGetArray(patch->localX, &localX);CHKERRQ(ierr); 1539 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, globalX, localX);CHKERRQ(ierr); 1540 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, globalX, localX);CHKERRQ(ierr); 1541 ierr = VecRestoreArrayRead(x, &globalX);CHKERRQ(ierr); 1542 ierr = VecRestoreArray(patch->localX, &localX);CHKERRQ(ierr); 1543 1544 ierr = VecSet(patch->localY, 0.0);CHKERRQ(ierr); 1545 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);CHKERRQ(ierr); 1546 for (sweep = 0; sweep < nsweep; sweep++) { 1547 for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) { 1548 PetscInt i = patch->user_patches ? iterationSet[j] : j; 1549 PetscInt start, len; 1550 1551 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);CHKERRQ(ierr); 1552 ierr = PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);CHKERRQ(ierr); 1553 /* TODO: Squash out these guys in the setup as well. */ 1554 if (len <= 0) continue; 1555 /* TODO: Do we need different scatters for X and Y? */ 1556 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localX, patch->patchX[i], INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1557 if (!patch->save_operators) { 1558 Mat mat; 1559 1560 ierr = PCPatchCreateMatrix_Private(pc, i, &mat);CHKERRQ(ierr); 1561 /* Populate operator here. */ 1562 ierr = PCPatchComputeOperator_Private(pc, mat, i);CHKERRQ(ierr); 1563 ierr = KSPSetOperators(patch->ksp[i], mat, mat);CHKERRQ(ierr); 1564 /* Drop reference so the KSPSetOperators below will blow it away. */ 1565 ierr = MatDestroy(&mat);CHKERRQ(ierr); 1566 } 1567 ierr = PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 1568 ierr = KSPSolve(patch->ksp[i], patch->patchX[i], patch->patchY[i]);CHKERRQ(ierr); 1569 ierr = PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 1570 1571 if (!patch->save_operators) { 1572 PC pc; 1573 ierr = KSPSetOperators(patch->ksp[i], NULL, NULL);CHKERRQ(ierr); 1574 ierr = KSPGetPC(patch->ksp[i], &pc);CHKERRQ(ierr); 1575 /* Destroy PC context too, otherwise the factored matrix hangs around. */ 1576 ierr = PCReset(pc);CHKERRQ(ierr); 1577 } 1578 1579 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchY[i], patch->localY, ADD_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 1580 } 1581 } 1582 if (patch->user_patches) {ierr = ISRestoreIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr);} 1583 /* XXX: should we do this on the global vector? */ 1584 if (patch->partition_of_unity) { 1585 ierr = VecPointwiseMult(patch->localY, patch->localY, patch->dof_weights);CHKERRQ(ierr); 1586 } 1587 /* Now patch->localY contains the solution of the patch solves, so we need to combine them all. */ 1588 ierr = VecSet(y, 0.0);CHKERRQ(ierr); 1589 ierr = VecGetArray(y, &globalY);CHKERRQ(ierr); 1590 ierr = VecGetArrayRead(patch->localY, &localY);CHKERRQ(ierr); 1591 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, localY, globalY, MPI_SUM);CHKERRQ(ierr); 1592 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, localY, globalY, MPI_SUM);CHKERRQ(ierr); 1593 ierr = VecRestoreArrayRead(patch->localY, &localY);CHKERRQ(ierr); 1594 1595 /* Now we need to send the global BC values through */ 1596 ierr = VecGetArrayRead(x, &globalX);CHKERRQ(ierr); 1597 ierr = ISGetSize(patch->globalBcNodes, &numBcs);CHKERRQ(ierr); 1598 ierr = ISGetIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 1599 ierr = VecGetLocalSize(x, &n);CHKERRQ(ierr); 1600 for (bc = 0; bc < numBcs; ++bc) { 1601 const PetscInt idx = bcNodes[bc]; 1602 if (idx < n) globalY[idx] = globalX[idx]; 1603 } 1604 1605 ierr = ISRestoreIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 1606 ierr = VecRestoreArrayRead(x, &globalX);CHKERRQ(ierr); 1607 ierr = VecRestoreArray(y, &globalY);CHKERRQ(ierr); 1608 1609 ierr = PetscOptionsPopGetViewerOff();CHKERRQ(ierr); 1610 ierr = PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 1611 PetscFunctionReturn(0); 1612 } 1613 1614 static PetscErrorCode PCReset_PATCH(PC pc) 1615 { 1616 PC_PATCH *patch = (PC_PATCH *) pc->data; 1617 PetscInt i; 1618 PetscErrorCode ierr; 1619 1620 PetscFunctionBegin; 1621 /* TODO: Get rid of all these ifs */ 1622 ierr = PetscSFDestroy(&patch->defaultSF);CHKERRQ(ierr); 1623 ierr = PetscSectionDestroy(&patch->cellCounts);CHKERRQ(ierr); 1624 ierr = PetscSectionDestroy(&patch->pointCounts);CHKERRQ(ierr); 1625 ierr = PetscSectionDestroy(&patch->cellNumbering);CHKERRQ(ierr); 1626 ierr = PetscSectionDestroy(&patch->gtolCounts);CHKERRQ(ierr); 1627 ierr = ISDestroy(&patch->gtol);CHKERRQ(ierr); 1628 ierr = ISDestroy(&patch->cells);CHKERRQ(ierr); 1629 ierr = ISDestroy(&patch->points);CHKERRQ(ierr); 1630 ierr = ISDestroy(&patch->dofs);CHKERRQ(ierr); 1631 ierr = ISDestroy(&patch->offs);CHKERRQ(ierr); 1632 ierr = PetscSectionDestroy(&patch->patchSection);CHKERRQ(ierr); 1633 ierr = ISDestroy(&patch->ghostBcNodes);CHKERRQ(ierr); 1634 ierr = ISDestroy(&patch->globalBcNodes);CHKERRQ(ierr); 1635 1636 if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscSectionDestroy(&patch->dofSection[i]);CHKERRQ(ierr);} 1637 ierr = PetscFree(patch->dofSection);CHKERRQ(ierr); 1638 ierr = PetscFree(patch->bs);CHKERRQ(ierr); 1639 ierr = PetscFree(patch->nodesPerCell);CHKERRQ(ierr); 1640 if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscFree(patch->cellNodeMap[i]);CHKERRQ(ierr);} 1641 ierr = PetscFree(patch->cellNodeMap);CHKERRQ(ierr); 1642 ierr = PetscFree(patch->subspaceOffsets);CHKERRQ(ierr); 1643 1644 if (patch->ksp) { 1645 for (i = 0; i < patch->npatch; ++i) {ierr = KSPReset(patch->ksp[i]);CHKERRQ(ierr);} 1646 } 1647 1648 if (patch->subspaces_to_exclude) { 1649 PetscHSetIDestroy(&patch->subspaces_to_exclude); 1650 } 1651 1652 ierr = VecDestroy(&patch->localX);CHKERRQ(ierr); 1653 ierr = VecDestroy(&patch->localY);CHKERRQ(ierr); 1654 if (patch->patchX) { 1655 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchX[i]);CHKERRQ(ierr);} 1656 ierr = PetscFree(patch->patchX);CHKERRQ(ierr); 1657 } 1658 if (patch->patchY) { 1659 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchY[i]);CHKERRQ(ierr);} 1660 ierr = PetscFree(patch->patchY);CHKERRQ(ierr); 1661 } 1662 ierr = VecDestroy(&patch->dof_weights);CHKERRQ(ierr); 1663 if (patch->patch_dof_weights) { 1664 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patch_dof_weights[i]);CHKERRQ(ierr);} 1665 ierr = PetscFree(patch->patch_dof_weights);CHKERRQ(ierr); 1666 } 1667 if (patch->mat) { 1668 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->mat[i]);CHKERRQ(ierr);} 1669 ierr = PetscFree(patch->mat);CHKERRQ(ierr); 1670 } 1671 ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr); 1672 if (patch->userIS) { 1673 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->userIS[i]);CHKERRQ(ierr);} 1674 ierr = PetscFree(patch->userIS);CHKERRQ(ierr); 1675 } 1676 patch->bs = 0; 1677 patch->cellNodeMap = NULL; 1678 patch->nsubspaces = 0; 1679 ierr = ISDestroy(&patch->iterationSet);CHKERRQ(ierr); 1680 1681 ierr = PetscViewerDestroy(&patch->viewerSection);CHKERRQ(ierr); 1682 PetscFunctionReturn(0); 1683 } 1684 1685 static PetscErrorCode PCDestroy_PATCH(PC pc) 1686 { 1687 PC_PATCH *patch = (PC_PATCH *) pc->data; 1688 PetscInt i; 1689 PetscErrorCode ierr; 1690 1691 PetscFunctionBegin; 1692 ierr = PCReset_PATCH(pc);CHKERRQ(ierr); 1693 if (patch->ksp) { 1694 for (i = 0; i < patch->npatch; ++i) {ierr = KSPDestroy(&patch->ksp[i]);CHKERRQ(ierr);} 1695 ierr = PetscFree(patch->ksp);CHKERRQ(ierr); 1696 } 1697 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1698 PetscFunctionReturn(0); 1699 } 1700 1701 static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc) 1702 { 1703 PC_PATCH *patch = (PC_PATCH *) pc->data; 1704 PCPatchConstructType patchConstructionType = PC_PATCH_STAR; 1705 char sub_mat_type[PETSC_MAX_PATH_LEN]; 1706 const char *prefix; 1707 PetscBool flg, dimflg, codimflg; 1708 MPI_Comm comm; 1709 PetscInt *ifields, nfields, k; 1710 PetscErrorCode ierr; 1711 1712 PetscFunctionBegin; 1713 ierr = PetscObjectGetComm((PetscObject) pc, &comm);CHKERRQ(ierr); 1714 ierr = PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);CHKERRQ(ierr); 1715 ierr = PetscOptionsHead(PetscOptionsObject, "Vertex-patch Additive Schwarz options");CHKERRQ(ierr); 1716 ierr = PetscOptionsBool("-pc_patch_save_operators", "Store all patch operators for lifetime of PC?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg);CHKERRQ(ierr); 1717 ierr = PetscOptionsBool("-pc_patch_partition_of_unity", "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg);CHKERRQ(ierr); 1718 ierr = PetscOptionsInt("-pc_patch_construct_dim", "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg);CHKERRQ(ierr); 1719 ierr = PetscOptionsInt("-pc_patch_construct_codim", "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg);CHKERRQ(ierr); 1720 if (dimflg && codimflg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");CHKERRQ(ierr); 1721 ierr = PetscOptionsEnum("-pc_patch_construct_type", "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum) patchConstructionType, (PetscEnum *) &patchConstructionType, &flg);CHKERRQ(ierr); 1722 if (flg) {ierr = PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL);CHKERRQ(ierr);} 1723 ierr = PetscOptionsInt("-pc_patch_vanka_dim", "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg);CHKERRQ(ierr); 1724 ierr = PetscOptionsInt("-pc_patch_ignore_dim", "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg);CHKERRQ(ierr); 1725 ierr = PetscOptionsFList("-pc_patch_sub_mat_type", "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg);CHKERRQ(ierr); 1726 if (flg) {ierr = PCPatchSetSubMatType(pc, sub_mat_type);CHKERRQ(ierr);} 1727 ierr = PetscOptionsBool("-pc_patch_symmetrise_sweep", "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg);CHKERRQ(ierr); 1728 1729 /* If the user has set the number of subspaces, use that for the buffer size, 1730 otherwise use a large number */ 1731 if (patch->nsubspaces <= 0) { 1732 nfields = 128; 1733 } else { 1734 nfields = patch->nsubspaces; 1735 } 1736 ierr = PetscMalloc1(nfields, &ifields);CHKERRQ(ierr); 1737 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_patch_exclude_subspaces",ifields,&nfields,&flg);CHKERRQ(ierr); 1738 if (flg && (patchConstructionType == PC_PATCH_USER)) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "We cannot support excluding a subspace with user patches because we do not index patches with a mesh point"); 1739 if (flg) { 1740 PetscHSetIClear(patch->subspaces_to_exclude); 1741 for (k = 0; k < nfields; k++) { 1742 PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]); 1743 } 1744 } 1745 ierr = PetscFree(ifields);CHKERRQ(ierr); 1746 1747 ierr = PetscOptionsBool("-pc_patch_patches_view", "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);CHKERRQ(ierr); 1748 ierr = PetscOptionsGetViewer(comm, prefix, "-pc_patch_cells_view", &patch->viewerCells, &patch->formatCells, &patch->viewCells);CHKERRQ(ierr); 1749 ierr = PetscOptionsGetViewer(comm, prefix, "-pc_patch_points_view", &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);CHKERRQ(ierr); 1750 ierr = PetscOptionsGetViewer(comm, prefix, "-pc_patch_section_view", &patch->viewerSection, &patch->formatSection, &patch->viewSection);CHKERRQ(ierr); 1751 ierr = PetscOptionsGetViewer(comm, prefix, "-pc_patch_mat_view", &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);CHKERRQ(ierr); 1752 ierr = PetscOptionsTail();CHKERRQ(ierr); 1753 patch->optionsSet = PETSC_TRUE; 1754 PetscFunctionReturn(0); 1755 } 1756 1757 static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc) 1758 { 1759 PC_PATCH *patch = (PC_PATCH*) pc->data; 1760 KSPConvergedReason reason; 1761 PetscInt i; 1762 PetscErrorCode ierr; 1763 1764 PetscFunctionBegin; 1765 for (i = 0; i < patch->npatch; ++i) { 1766 ierr = KSPSetUp(patch->ksp[i]);CHKERRQ(ierr); 1767 ierr = KSPGetConvergedReason(patch->ksp[i], &reason);CHKERRQ(ierr); 1768 if (reason == KSP_DIVERGED_PCSETUP_FAILED) pc->failedreason = PC_SUBPC_ERROR; 1769 } 1770 PetscFunctionReturn(0); 1771 } 1772 1773 static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer) 1774 { 1775 PC_PATCH *patch = (PC_PATCH *) pc->data; 1776 PetscViewer sviewer; 1777 PetscBool isascii; 1778 PetscMPIInt rank; 1779 PetscErrorCode ierr; 1780 1781 PetscFunctionBegin; 1782 /* TODO Redo tabbing with set tbas in new style */ 1783 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 1784 if (!isascii) PetscFunctionReturn(0); 1785 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);CHKERRQ(ierr); 1786 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1787 ierr = PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);CHKERRQ(ierr); 1788 ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");CHKERRQ(ierr); 1789 if (patch->partition_of_unity) {ierr = PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");CHKERRQ(ierr);} 1790 else {ierr = PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");CHKERRQ(ierr);} 1791 if (patch->symmetrise_sweep) {ierr = PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");CHKERRQ(ierr);} 1792 else {ierr = PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");CHKERRQ(ierr);} 1793 if (!patch->save_operators) {ierr = PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");CHKERRQ(ierr);} 1794 else {ierr = PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");CHKERRQ(ierr);} 1795 if (patch->patchconstructop == PCPatchConstruct_Star) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");CHKERRQ(ierr);} 1796 else if (patch->patchconstructop == PCPatchConstruct_Vanka) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");CHKERRQ(ierr);} 1797 else if (patch->patchconstructop == PCPatchConstruct_User) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");CHKERRQ(ierr);} 1798 else {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");CHKERRQ(ierr);} 1799 ierr = PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");CHKERRQ(ierr); 1800 if (patch->ksp) { 1801 ierr = PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 1802 if (!rank) { 1803 ierr = PetscViewerASCIIPushTab(sviewer);CHKERRQ(ierr); 1804 ierr = KSPView(patch->ksp[0], sviewer);CHKERRQ(ierr); 1805 ierr = PetscViewerASCIIPopTab(sviewer);CHKERRQ(ierr); 1806 } 1807 ierr = PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 1808 } else { 1809 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1810 ierr = PetscViewerASCIIPrintf(viewer, "KSP not yet set.\n");CHKERRQ(ierr); 1811 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1812 } 1813 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1814 PetscFunctionReturn(0); 1815 } 1816 1817 /*MC 1818 PCPATCH - A PC object that encapsulates flexible definition of blocks for overlapping and non-overlapping 1819 small block additive preconditioners. Block definition is based on topology from 1820 a DM and equation numbering from a PetscSection. 1821 1822 Options Database Keys: 1823 + -pc_patch_cells_view - Views the process local cell numbers for each patch 1824 . -pc_patch_points_view - Views the process local mesh point numbers for each patch 1825 . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch 1826 . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary 1827 - -pc_patch_sub_mat_view - Views the matrix associated with each patch 1828 1829 Level: intermediate 1830 1831 .seealso: PCType, PCCreate(), PCSetType() 1832 M*/ 1833 PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc) 1834 { 1835 PC_PATCH *patch; 1836 PetscErrorCode ierr; 1837 1838 PetscFunctionBegin; 1839 ierr = PetscNewLog(pc, &patch);CHKERRQ(ierr); 1840 1841 if (patch->subspaces_to_exclude) { 1842 PetscHSetIDestroy(&patch->subspaces_to_exclude); 1843 } 1844 PetscHSetICreate(&patch->subspaces_to_exclude); 1845 1846 /* Set some defaults */ 1847 patch->combined = PETSC_FALSE; 1848 patch->save_operators = PETSC_TRUE; 1849 patch->partition_of_unity = PETSC_FALSE; 1850 patch->codim = -1; 1851 patch->dim = -1; 1852 patch->vankadim = -1; 1853 patch->ignoredim = -1; 1854 patch->patchconstructop = PCPatchConstruct_Star; 1855 patch->symmetrise_sweep = PETSC_FALSE; 1856 patch->npatch = 0; 1857 patch->userIS = NULL; 1858 patch->optionsSet = PETSC_FALSE; 1859 patch->iterationSet = NULL; 1860 patch->user_patches = PETSC_FALSE; 1861 ierr = PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 1862 patch->viewPatches = PETSC_FALSE; 1863 patch->viewCells = PETSC_FALSE; 1864 patch->viewPoints = PETSC_FALSE; 1865 patch->viewSection = PETSC_FALSE; 1866 patch->viewMatrix = PETSC_FALSE; 1867 1868 pc->data = (void *) patch; 1869 pc->ops->apply = PCApply_PATCH; 1870 pc->ops->applytranspose = 0; /* PCApplyTranspose_PATCH; */ 1871 pc->ops->setup = PCSetUp_PATCH; 1872 pc->ops->reset = PCReset_PATCH; 1873 pc->ops->destroy = PCDestroy_PATCH; 1874 pc->ops->setfromoptions = PCSetFromOptions_PATCH; 1875 pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH; 1876 pc->ops->view = PCView_PATCH; 1877 pc->ops->applyrichardson = 0; 1878 1879 PetscFunctionReturn(0); 1880 } 1881