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