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, mat, cellIS, n, u, ctx) 470 + pc - The PC 471 + point - The point 472 + mat - The patch matrix 473 + cellIS - An array of the cell numbers 474 + n - The size of g2l 475 + g2l - The global to local dof translation table 476 + ctx - The user context 477 and can assume that the matrix entries have been set to zero before the call. 478 479 .seealso: PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo() 480 @*/ 481 PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC, PetscInt, Mat, IS, PetscInt, const PetscInt *, void *), void *ctx) 482 { 483 PC_PATCH *patch = (PC_PATCH *) pc->data; 484 485 PetscFunctionBegin; 486 patch->usercomputeop = func; 487 patch->usercomputectx = ctx; 488 PetscFunctionReturn(0); 489 } 490 491 /* On entry, ht contains the topological entities whose dofs we are responsible for solving for; 492 on exit, cht contains all the topological entities we need to compute their residuals. 493 In full generality this should incorporate knowledge of the sparsity pattern of the matrix; 494 here we assume a standard FE sparsity pattern.*/ 495 /* TODO: Use DMPlexGetAdjacency() */ 496 static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht) 497 { 498 DM dm; 499 PetscHashIter hi; 500 PetscInt point; 501 PetscInt *star = NULL, *closure = NULL; 502 PetscInt ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci; 503 PetscErrorCode ierr; 504 505 PetscFunctionBegin; 506 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 507 ierr = PCPatchGetIgnoreDim(pc, &ignoredim);CHKERRQ(ierr); 508 if (ignoredim >= 0) {ierr = DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd);CHKERRQ(ierr);} 509 ierr = PetscHSetIClear(cht);CHKERRQ(ierr); 510 PetscHashIterBegin(ht, hi); 511 while (!PetscHashIterAtEnd(ht, hi)) { 512 513 PetscHashIterGetKey(ht, hi, point); 514 PetscHashIterNext(ht, hi); 515 516 /* Loop over all the cells that this point connects to */ 517 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 518 for (si = 0; si < starSize*2; si += 2) { 519 const PetscInt ownedpoint = star[si]; 520 /* TODO Check for point in cht before running through closure again */ 521 /* now loop over all entities in the closure of that cell */ 522 ierr = DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 523 for (ci = 0; ci < closureSize*2; ci += 2) { 524 const PetscInt seenpoint = closure[ci]; 525 if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue; 526 ierr = PetscHSetIAdd(cht, seenpoint);CHKERRQ(ierr); 527 } 528 } 529 } 530 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);CHKERRQ(ierr); 531 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_FALSE, NULL, &star);CHKERRQ(ierr); 532 PetscFunctionReturn(0); 533 } 534 535 static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off) 536 { 537 PetscErrorCode ierr; 538 539 PetscFunctionBegin; 540 if (combined) { 541 if (f < 0) { 542 if (dof) {ierr = PetscSectionGetDof(dofSection[0], p, dof);CHKERRQ(ierr);} 543 if (off) {ierr = PetscSectionGetOffset(dofSection[0], p, off);CHKERRQ(ierr);} 544 } else { 545 if (dof) {ierr = PetscSectionGetFieldDof(dofSection[0], p, f, dof);CHKERRQ(ierr);} 546 if (off) {ierr = PetscSectionGetFieldOffset(dofSection[0], p, f, off);CHKERRQ(ierr);} 547 } 548 } else { 549 if (f < 0) { 550 PC_PATCH *patch = (PC_PATCH *) pc->data; 551 PetscInt fdof, g; 552 553 if (dof) { 554 *dof = 0; 555 for (g = 0; g < patch->nsubspaces; ++g) { 556 ierr = PetscSectionGetDof(dofSection[g], p, &fdof);CHKERRQ(ierr); 557 *dof += fdof; 558 } 559 } 560 if (off) { 561 *off = 0; 562 for (g = 0; g < patch->nsubspaces; ++g) { 563 ierr = PetscSectionGetOffset(dofSection[g], p, &fdof);CHKERRQ(ierr); 564 *off += fdof; 565 } 566 } 567 } else { 568 if (dof) {ierr = PetscSectionGetDof(dofSection[f], p, dof);CHKERRQ(ierr);} 569 if (off) {ierr = PetscSectionGetOffset(dofSection[f], p, off);CHKERRQ(ierr);} 570 } 571 } 572 PetscFunctionReturn(0); 573 } 574 575 /* Given a hash table with a set of topological entities (pts), compute the degrees of 576 freedom in global concatenated numbering on those entities. 577 For Vanka smoothing, this needs to do something special: ignore dofs of the 578 constraint subspace on entities that aren't the base entity we're building the patch 579 around. */ 580 static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI* subspaces_to_exclude) 581 { 582 PC_PATCH *patch = (PC_PATCH *) pc->data; 583 PetscHashIter hi; 584 PetscInt ldof, loff; 585 PetscInt k, p; 586 PetscErrorCode ierr; 587 588 PetscFunctionBegin; 589 ierr = PetscHSetIClear(dofs);CHKERRQ(ierr); 590 for (k = 0; k < patch->nsubspaces; ++k) { 591 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 592 PetscInt bs = patch->bs[k]; 593 PetscInt j, l; 594 595 if (subspaces_to_exclude != NULL) { 596 PetscBool should_exclude_k = PETSC_FALSE; 597 PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k); 598 if (should_exclude_k) { 599 /* only get this subspace dofs at the base entity, not any others */ 600 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff);CHKERRQ(ierr); 601 if (0 == ldof) continue; 602 for (j = loff; j < ldof + loff; ++j) { 603 for (l = 0; l < bs; ++l) { 604 PetscInt dof = bs*j + l + subspaceOffset; 605 ierr = PetscHSetIAdd(dofs, dof);CHKERRQ(ierr); 606 } 607 } 608 continue; /* skip the other dofs of this subspace */ 609 } 610 } 611 612 PetscHashIterBegin(pts, hi); 613 while (!PetscHashIterAtEnd(pts, hi)) { 614 PetscHashIterGetKey(pts, hi, p); 615 PetscHashIterNext(pts, hi); 616 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);CHKERRQ(ierr); 617 if (0 == ldof) continue; 618 for (j = loff; j < ldof + loff; ++j) { 619 for (l = 0; l < bs; ++l) { 620 PetscInt dof = bs*j + l + subspaceOffset; 621 ierr = PetscHSetIAdd(dofs, dof);CHKERRQ(ierr); 622 } 623 } 624 } 625 } 626 PetscFunctionReturn(0); 627 } 628 629 /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */ 630 static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C) 631 { 632 PetscHashIter hi; 633 PetscInt key; 634 PetscBool flg; 635 PetscErrorCode ierr; 636 637 PetscFunctionBegin; 638 ierr = PetscHSetIClear(C);CHKERRQ(ierr); 639 PetscHashIterBegin(B, hi); 640 while (!PetscHashIterAtEnd(B, hi)) { 641 PetscHashIterGetKey(B, hi, key); 642 PetscHashIterNext(B, hi); 643 ierr = PetscHSetIHas(A, key, &flg);CHKERRQ(ierr); 644 if (!flg) {ierr = PetscHSetIAdd(C, key);CHKERRQ(ierr);} 645 } 646 PetscFunctionReturn(0); 647 } 648 649 /* 650 * PCPatchCreateCellPatches - create patches. 651 * 652 * Input Parameters: 653 * + dm - The DMPlex object defining the mesh 654 * 655 * Output Parameters: 656 * + cellCounts - Section with counts of cells around each vertex 657 * . cells - IS of the cell point indices of cells in each patch 658 * . pointCounts - Section with counts of cells around each vertex 659 * - point - IS of the cell point indices of cells in each patch 660 */ 661 static PetscErrorCode PCPatchCreateCellPatches(PC pc) 662 { 663 PC_PATCH *patch = (PC_PATCH *) pc->data; 664 DMLabel ghost = NULL; 665 DM dm, plex; 666 PetscHSetI ht, cht; 667 PetscSection cellCounts, pointCounts; 668 PetscInt *cellsArray, *pointsArray; 669 PetscInt numCells, numPoints; 670 const PetscInt *leaves; 671 PetscInt nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, v; 672 PetscBool isFiredrake; 673 PetscErrorCode ierr; 674 675 PetscFunctionBegin; 676 /* Used to keep track of the cells in the patch. */ 677 ierr = PetscHSetICreate(&ht);CHKERRQ(ierr); 678 ierr = PetscHSetICreate(&cht);CHKERRQ(ierr); 679 680 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 681 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC\n"); 682 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 683 ierr = DMPlexGetChart(plex, &pStart, &pEnd);CHKERRQ(ierr); 684 ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr); 685 686 if (patch->user_patches) { 687 ierr = patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);CHKERRQ(ierr); 688 vStart = 0; vEnd = patch->npatch; 689 } else if (patch->codim < 0) { 690 if (patch->dim < 0) {ierr = DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);CHKERRQ(ierr);} 691 else {ierr = DMPlexGetDepthStratum(plex, patch->dim, &vStart, &vEnd);CHKERRQ(ierr);} 692 } else {ierr = DMPlexGetHeightStratum(plex, patch->codim, &vStart, &vEnd);CHKERRQ(ierr);} 693 patch->npatch = vEnd - vStart; 694 695 /* These labels mark the owned points. We only create patches around points that this process owns. */ 696 ierr = DMHasLabel(dm, "pyop2_ghost", &isFiredrake);CHKERRQ(ierr); 697 if (isFiredrake) { 698 ierr = DMGetLabel(dm, "pyop2_ghost", &ghost);CHKERRQ(ierr); 699 ierr = DMLabelCreateIndex(ghost, pStart, pEnd);CHKERRQ(ierr); 700 } else { 701 PetscSF sf; 702 703 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 704 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 705 nleaves = PetscMax(nleaves, 0); 706 } 707 708 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts);CHKERRQ(ierr); 709 ierr = PetscObjectSetName((PetscObject) patch->cellCounts, "Patch Cell Layout");CHKERRQ(ierr); 710 cellCounts = patch->cellCounts; 711 ierr = PetscSectionSetChart(cellCounts, vStart, vEnd);CHKERRQ(ierr); 712 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts);CHKERRQ(ierr); 713 ierr = PetscObjectSetName((PetscObject) patch->pointCounts, "Patch Point Layout");CHKERRQ(ierr); 714 pointCounts = patch->pointCounts; 715 ierr = PetscSectionSetChart(pointCounts, vStart, vEnd);CHKERRQ(ierr); 716 /* Count cells and points in the patch surrounding each entity */ 717 for (v = vStart; v < vEnd; ++v) { 718 PetscHashIter hi; 719 PetscInt chtSize, loc = -1; 720 PetscBool flg; 721 722 if (!patch->user_patches) { 723 if (ghost) {ierr = DMLabelHasPoint(ghost, v, &flg);CHKERRQ(ierr);} 724 else {ierr = PetscFindInt(v, nleaves, leaves, &loc);CHKERRQ(ierr); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;} 725 /* Not an owned entity, don't make a cell patch. */ 726 if (flg) continue; 727 } 728 729 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 730 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 731 ierr = PetscHSetIGetSize(cht, &chtSize);CHKERRQ(ierr); 732 /* empty patch, continue */ 733 if (chtSize == 0) continue; 734 735 /* safe because size(cht) > 0 from above */ 736 PetscHashIterBegin(cht, hi); 737 while (!PetscHashIterAtEnd(cht, hi)) { 738 PetscInt point, pdof; 739 740 PetscHashIterGetKey(cht, hi, point); 741 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 742 if (pdof) {ierr = PetscSectionAddDof(pointCounts, v, 1);CHKERRQ(ierr);} 743 if (point >= cStart && point < cEnd) {ierr = PetscSectionAddDof(cellCounts, v, 1);CHKERRQ(ierr);} 744 PetscHashIterNext(cht, hi); 745 } 746 } 747 if (isFiredrake) {ierr = DMLabelDestroyIndex(ghost);CHKERRQ(ierr);} 748 749 ierr = PetscSectionSetUp(cellCounts);CHKERRQ(ierr); 750 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 751 ierr = PetscMalloc1(numCells, &cellsArray);CHKERRQ(ierr); 752 ierr = PetscSectionSetUp(pointCounts);CHKERRQ(ierr); 753 ierr = PetscSectionGetStorageSize(pointCounts, &numPoints);CHKERRQ(ierr); 754 ierr = PetscMalloc1(numPoints, &pointsArray);CHKERRQ(ierr); 755 756 /* Now that we know how much space we need, run through again and actually remember the cells. */ 757 for (v = vStart; v < vEnd; v++ ) { 758 PetscHashIter hi; 759 PetscInt dof, off, cdof, coff, pdof, n = 0, cn = 0; 760 761 ierr = PetscSectionGetDof(pointCounts, v, &dof);CHKERRQ(ierr); 762 ierr = PetscSectionGetOffset(pointCounts, v, &off);CHKERRQ(ierr); 763 ierr = PetscSectionGetDof(cellCounts, v, &cdof);CHKERRQ(ierr); 764 ierr = PetscSectionGetOffset(cellCounts, v, &coff);CHKERRQ(ierr); 765 if (dof <= 0) continue; 766 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 767 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 768 PetscHashIterBegin(cht, hi); 769 while (!PetscHashIterAtEnd(cht, hi)) { 770 PetscInt point; 771 772 PetscHashIterGetKey(cht, hi, point); 773 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 774 if (pdof) {pointsArray[off + n++] = point;} 775 if (point >= cStart && point < cEnd) {cellsArray[coff + cn++] = point;} 776 PetscHashIterNext(cht, hi); 777 } 778 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); 779 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); 780 } 781 ierr = PetscHSetIDestroy(&ht);CHKERRQ(ierr); 782 ierr = PetscHSetIDestroy(&cht);CHKERRQ(ierr); 783 ierr = DMDestroy(&plex);CHKERRQ(ierr); 784 785 ierr = ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells);CHKERRQ(ierr); 786 ierr = PetscObjectSetName((PetscObject) patch->cells, "Patch Cells");CHKERRQ(ierr); 787 if (patch->viewCells) { 788 ierr = ObjectView((PetscObject) patch->cellCounts, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 789 ierr = ObjectView((PetscObject) patch->cells, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 790 } 791 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points);CHKERRQ(ierr); 792 ierr = PetscObjectSetName((PetscObject) patch->points, "Patch Points");CHKERRQ(ierr); 793 if (patch->viewPoints) { 794 ierr = ObjectView((PetscObject) patch->pointCounts, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 795 ierr = ObjectView((PetscObject) patch->points, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 796 } 797 PetscFunctionReturn(0); 798 } 799 800 /* 801 * PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches 802 * 803 * Input Parameters: 804 * + dm - The DMPlex object defining the mesh 805 * . cellCounts - Section with counts of cells around each vertex 806 * . cells - IS of the cell point indices of cells in each patch 807 * . cellNumbering - Section mapping plex cell points to Firedrake cell indices. 808 * . nodesPerCell - number of nodes per cell. 809 * - cellNodeMap - map from cells to node indices (nodesPerCell * numCells) 810 * 811 * Output Parameters: 812 * + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering 813 * . gtolCounts - Section with counts of dofs per cell patch 814 * - gtol - IS mapping from global dofs to local dofs for each patch. 815 */ 816 static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc) 817 { 818 PC_PATCH *patch = (PC_PATCH *) pc->data; 819 PetscSection cellCounts = patch->cellCounts; 820 PetscSection pointCounts = patch->pointCounts; 821 PetscSection gtolCounts; 822 IS cells = patch->cells; 823 IS points = patch->points; 824 PetscSection cellNumbering = patch->cellNumbering; 825 PetscInt Nf = patch->nsubspaces; 826 PetscInt numCells, numPoints; 827 PetscInt numDofs; 828 PetscInt numGlobalDofs; 829 PetscInt totalDofsPerCell = patch->totalDofsPerCell; 830 PetscInt vStart, vEnd, v; 831 const PetscInt *cellsArray, *pointsArray; 832 PetscInt *newCellsArray = NULL; 833 PetscInt *dofsArray = NULL; 834 PetscInt *offsArray = NULL; 835 PetscInt *asmArray = NULL; 836 PetscInt *globalDofsArray = NULL; 837 PetscInt globalIndex = 0; 838 PetscInt key = 0; 839 PetscInt asmKey = 0; 840 DM dm = NULL; 841 const PetscInt *bcNodes = NULL; 842 PetscHMapI ht; 843 PetscHSetI globalBcs; 844 PetscInt numBcs; 845 PetscHSetI ownedpts, seenpts, owneddofs, seendofs, artificialbcs; 846 PetscInt pStart, pEnd, p, i; 847 PetscErrorCode ierr; 848 849 PetscFunctionBegin; 850 851 ierr = PCGetDM(pc, &dm); CHKERRQ(ierr); 852 /* dofcounts section is cellcounts section * dofPerCell */ 853 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 854 ierr = PetscSectionGetStorageSize(patch->pointCounts, &numPoints);CHKERRQ(ierr); 855 numDofs = numCells * totalDofsPerCell; 856 ierr = PetscMalloc1(numDofs, &dofsArray);CHKERRQ(ierr); 857 ierr = PetscMalloc1(numPoints*Nf, &offsArray);CHKERRQ(ierr); 858 ierr = PetscMalloc1(numDofs, &asmArray);CHKERRQ(ierr); 859 ierr = PetscMalloc1(numCells, &newCellsArray);CHKERRQ(ierr); 860 ierr = PetscSectionGetChart(cellCounts, &vStart, &vEnd);CHKERRQ(ierr); 861 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);CHKERRQ(ierr); 862 gtolCounts = patch->gtolCounts; 863 ierr = PetscSectionSetChart(gtolCounts, vStart, vEnd);CHKERRQ(ierr); 864 ierr = PetscObjectSetName((PetscObject) patch->gtolCounts, "Patch Global Index Section");CHKERRQ(ierr); 865 866 /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet 867 conditions */ 868 ierr = PetscHSetICreate(&globalBcs);CHKERRQ(ierr); 869 ierr = ISGetIndices(patch->ghostBcNodes, &bcNodes); CHKERRQ(ierr); 870 ierr = ISGetSize(patch->ghostBcNodes, &numBcs); CHKERRQ(ierr); 871 for (i = 0; i < numBcs; ++i) { 872 ierr = PetscHSetIAdd(globalBcs, bcNodes[i]);CHKERRQ(ierr); /* these are already in concatenated numbering */ 873 } 874 ierr = ISRestoreIndices(patch->ghostBcNodes, &bcNodes); CHKERRQ(ierr); 875 ierr = ISDestroy(&patch->ghostBcNodes); CHKERRQ(ierr); /* memory optimisation */ 876 877 /* Hash tables for artificial BC construction */ 878 ierr = PetscHSetICreate(&ownedpts);CHKERRQ(ierr); 879 ierr = PetscHSetICreate(&seenpts);CHKERRQ(ierr); 880 ierr = PetscHSetICreate(&owneddofs);CHKERRQ(ierr); 881 ierr = PetscHSetICreate(&seendofs);CHKERRQ(ierr); 882 ierr = PetscHSetICreate(&artificialbcs);CHKERRQ(ierr); 883 884 ierr = ISGetIndices(cells, &cellsArray);CHKERRQ(ierr); 885 ierr = ISGetIndices(points, &pointsArray);CHKERRQ(ierr); 886 ierr = PetscHMapICreate(&ht);CHKERRQ(ierr); 887 for (v = vStart; v < vEnd; ++v) { 888 PetscInt localIndex = 0; 889 PetscInt dof, off, i, j, k, l; 890 891 ierr = PetscHMapIClear(ht);CHKERRQ(ierr); 892 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 893 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 894 if (dof <= 0) continue; 895 896 /* Calculate the global numbers of the artificial BC dofs here first */ 897 ierr = patch->patchconstructop((void*)patch, dm, v, ownedpts); CHKERRQ(ierr); 898 ierr = PCPatchCompleteCellPatch(pc, ownedpts, seenpts); CHKERRQ(ierr); 899 ierr = PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude); CHKERRQ(ierr); 900 ierr = PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL); CHKERRQ(ierr); 901 ierr = PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs); CHKERRQ(ierr); 902 if (patch->viewPatches) { 903 PetscHSetI globalbcdofs; 904 PetscHashIter hi; 905 MPI_Comm comm = PetscObjectComm((PetscObject)pc); 906 907 ierr = PetscHSetICreate(&globalbcdofs);CHKERRQ(ierr); 908 ierr = PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v); CHKERRQ(ierr); 909 PetscHashIterBegin(owneddofs, hi); 910 while (!PetscHashIterAtEnd(owneddofs, hi)) { 911 PetscInt globalDof; 912 913 PetscHashIterGetKey(owneddofs, hi, globalDof); 914 PetscHashIterNext(owneddofs, hi); 915 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 916 } 917 ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr); 918 ierr = PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v); CHKERRQ(ierr); 919 PetscHashIterBegin(seendofs, hi); 920 while (!PetscHashIterAtEnd(seendofs, hi)) { 921 PetscInt globalDof; 922 PetscBool flg; 923 924 PetscHashIterGetKey(seendofs, hi, globalDof); 925 PetscHashIterNext(seendofs, hi); 926 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 927 928 ierr = PetscHSetIHas(globalBcs, globalDof, &flg);CHKERRQ(ierr); 929 if (flg) {ierr = PetscHSetIAdd(globalbcdofs, globalDof);CHKERRQ(ierr);} 930 } 931 ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr); 932 ierr = PetscSynchronizedPrintf(comm, "Patch %d: global BCs:\n", v); CHKERRQ(ierr); 933 ierr = PetscHSetIGetSize(globalbcdofs, &numBcs);CHKERRQ(ierr); 934 if (numBcs > 0) { 935 PetscHashIterBegin(globalbcdofs, hi); 936 while (!PetscHashIterAtEnd(globalbcdofs, hi)) { 937 PetscInt globalDof; 938 PetscHashIterGetKey(globalbcdofs, hi, globalDof); 939 PetscHashIterNext(globalbcdofs, hi); 940 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof);CHKERRQ(ierr); 941 } 942 } 943 ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr); 944 ierr = PetscSynchronizedPrintf(comm, "Patch %d: artificial BCs:\n", v);CHKERRQ(ierr); 945 ierr = PetscHSetIGetSize(artificialbcs, &numBcs);CHKERRQ(ierr); 946 if (numBcs > 0) { 947 PetscHashIterBegin(artificialbcs, hi); 948 while (!PetscHashIterAtEnd(artificialbcs, hi)) { 949 PetscInt globalDof; 950 PetscHashIterGetKey(artificialbcs, hi, globalDof); 951 PetscHashIterNext(artificialbcs, hi); 952 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 953 } 954 } 955 ierr = PetscSynchronizedPrintf(comm, "\n\n"); CHKERRQ(ierr); 956 ierr = PetscHSetIDestroy(&globalbcdofs);CHKERRQ(ierr); 957 } 958 for (k = 0; k < patch->nsubspaces; ++k) { 959 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 960 PetscInt nodesPerCell = patch->nodesPerCell[k]; 961 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 962 PetscInt bs = patch->bs[k]; 963 964 for (i = off; i < off + dof; ++i) { 965 /* Walk over the cells in this patch. */ 966 const PetscInt c = cellsArray[i]; 967 PetscInt cell = c; 968 969 /* TODO Change this to an IS */ 970 if (cellNumbering) { 971 ierr = PetscSectionGetDof(cellNumbering, c, &cell);CHKERRQ(ierr); 972 if (cell <= 0) SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %D doesn't appear in cell numbering map", c); 973 ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr); 974 } 975 newCellsArray[i] = cell; 976 for (j = 0; j < nodesPerCell; ++j) { 977 /* For each global dof, map it into contiguous local storage. */ 978 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + subspaceOffset; 979 /* finally, loop over block size */ 980 for (l = 0; l < bs; ++l) { 981 PetscInt localDof; 982 PetscBool isGlobalBcDof, isArtificialBcDof; 983 984 /* first, check if this is either a globally enforced or locally enforced BC dof */ 985 ierr = PetscHSetIHas(globalBcs, globalDof + l, &isGlobalBcDof);CHKERRQ(ierr); 986 ierr = PetscHSetIHas(artificialbcs, globalDof + l, &isArtificialBcDof);CHKERRQ(ierr); 987 988 /* if it's either, don't ever give it a local dof number */ 989 if (isGlobalBcDof || isArtificialBcDof) { 990 dofsArray[globalIndex++] = -1; /* don't use this in assembly in this patch */ 991 } else { 992 ierr = PetscHMapIGet(ht, globalDof + l, &localDof);CHKERRQ(ierr); 993 if (localDof == -1) { 994 localDof = localIndex++; 995 ierr = PetscHMapISet(ht, globalDof + l, localDof);CHKERRQ(ierr); 996 } 997 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 998 /* And store. */ 999 dofsArray[globalIndex++] = localDof; 1000 } 1001 } 1002 } 1003 } 1004 } 1005 /* How many local dofs in this patch? */ 1006 ierr = PetscHMapIGetSize(ht, &dof);CHKERRQ(ierr); 1007 ierr = PetscSectionSetDof(gtolCounts, v, dof);CHKERRQ(ierr); 1008 } 1009 if (globalIndex != numDofs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%d) doesn't match found number (%d)", numDofs, globalIndex); 1010 ierr = PetscSectionSetUp(gtolCounts);CHKERRQ(ierr); 1011 ierr = PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);CHKERRQ(ierr); 1012 ierr = PetscMalloc1(numGlobalDofs, &globalDofsArray);CHKERRQ(ierr); 1013 1014 /* Now populate the global to local map. This could be merged into the above loop if we were willing to deal with reallocs. */ 1015 for (v = vStart; v < vEnd; ++v) { 1016 PetscHashIter hi; 1017 PetscInt dof, off, Np, ooff, i, j, k, l; 1018 1019 ierr = PetscHMapIClear(ht);CHKERRQ(ierr); 1020 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 1021 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 1022 ierr = PetscSectionGetDof(pointCounts, v, &Np);CHKERRQ(ierr); 1023 ierr = PetscSectionGetOffset(pointCounts, v, &ooff);CHKERRQ(ierr); 1024 if (dof <= 0) continue; 1025 1026 for (k = 0; k < patch->nsubspaces; ++k) { 1027 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1028 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1029 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1030 PetscInt bs = patch->bs[k]; 1031 PetscInt goff; 1032 1033 for (i = off; i < off + dof; ++i) { 1034 /* Reconstruct mapping of global-to-local on this patch. */ 1035 const PetscInt c = cellsArray[i]; 1036 PetscInt cell = c; 1037 1038 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 1039 for (j = 0; j < nodesPerCell; ++j) { 1040 for (l = 0; l < bs; ++l) { 1041 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1042 const PetscInt localDof = dofsArray[key++]; 1043 1044 if (localDof >= 0) {ierr = PetscHMapISet(ht, globalDof, localDof);CHKERRQ(ierr);} 1045 } 1046 } 1047 } 1048 1049 /* Shove it in the output data structure. */ 1050 ierr = PetscSectionGetOffset(gtolCounts, v, &goff);CHKERRQ(ierr); 1051 PetscHashIterBegin(ht, hi); 1052 while (!PetscHashIterAtEnd(ht, hi)) { 1053 PetscInt globalDof, localDof; 1054 1055 PetscHashIterGetKey(ht, hi, globalDof); 1056 PetscHashIterGetVal(ht, hi, localDof); 1057 if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof; 1058 PetscHashIterNext(ht, hi); 1059 } 1060 1061 for (p = 0; p < Np; ++p) { 1062 const PetscInt point = pointsArray[ooff + p]; 1063 PetscInt globalDof, localDof; 1064 1065 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);CHKERRQ(ierr); 1066 ierr = PetscHMapIGet(ht, globalDof, &localDof);CHKERRQ(ierr); 1067 offsArray[(ooff + p)*Nf + k] = localDof; 1068 } 1069 } 1070 1071 ierr = PetscHSetIDestroy(&globalBcs);CHKERRQ(ierr); 1072 ierr = PetscHSetIDestroy(&ownedpts);CHKERRQ(ierr); 1073 ierr = PetscHSetIDestroy(&seenpts);CHKERRQ(ierr); 1074 ierr = PetscHSetIDestroy(&owneddofs);CHKERRQ(ierr); 1075 ierr = PetscHSetIDestroy(&seendofs);CHKERRQ(ierr); 1076 ierr = PetscHSetIDestroy(&artificialbcs);CHKERRQ(ierr); 1077 1078 /* At this point, we have a hash table ht built that maps globalDof -> localDof. 1079 We need to create the dof table laid out cellwise first, then by subspace, 1080 as the assembler assembles cell-wise and we need to stuff the different 1081 contributions of the different function spaces to the right places. So we loop 1082 over cells, then over subspaces. */ 1083 if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */ 1084 for (i = off; i < off + dof; ++i) { 1085 const PetscInt c = cellsArray[i]; 1086 PetscInt cell = c; 1087 1088 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 1089 for (k = 0; k < patch->nsubspaces; ++k) { 1090 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1091 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1092 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1093 PetscInt bs = patch->bs[k]; 1094 1095 for (j = 0; j < nodesPerCell; ++j) { 1096 for (l = 0; l < bs; ++l) { 1097 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1098 PetscInt localDof; 1099 1100 ierr = PetscHMapIGet(ht, globalDof, &localDof);CHKERRQ(ierr); 1101 /* If it's not in the hash table, i.e. is a BC dof, 1102 then the PetscHSetIMap above gives -1, which matches 1103 exactly the convention for PETSc's matrix assembly to 1104 ignore the dof. So we don't need to do anything here */ 1105 asmArray[asmKey++] = localDof; 1106 } 1107 } 1108 } 1109 } 1110 } 1111 } 1112 if (1 == patch->nsubspaces) {ierr = PetscMemcpy(asmArray, dofsArray, numDofs * sizeof(PetscInt));CHKERRQ(ierr);} 1113 1114 ierr = PetscHMapIDestroy(&ht);CHKERRQ(ierr); 1115 ierr = ISRestoreIndices(cells, &cellsArray);CHKERRQ(ierr); 1116 ierr = ISRestoreIndices(points, &pointsArray);CHKERRQ(ierr); 1117 ierr = PetscFree(dofsArray);CHKERRQ(ierr); 1118 /* Create placeholder section for map from points to patch dofs */ 1119 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);CHKERRQ(ierr); 1120 ierr = PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);CHKERRQ(ierr); 1121 if (patch->combined) { 1122 PetscInt numFields; 1123 ierr = PetscSectionGetNumFields(patch->dofSection[0], &numFields);CHKERRQ(ierr); 1124 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); 1125 ierr = PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);CHKERRQ(ierr); 1126 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1127 for (p = pStart; p < pEnd; ++p) { 1128 PetscInt dof, fdof, f; 1129 1130 ierr = PetscSectionGetDof(patch->dofSection[0], p, &dof);CHKERRQ(ierr); 1131 ierr = PetscSectionSetDof(patch->patchSection, p, dof);CHKERRQ(ierr); 1132 for (f = 0; f < patch->nsubspaces; ++f) { 1133 ierr = PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof);CHKERRQ(ierr); 1134 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1135 } 1136 } 1137 } else { 1138 PetscInt pStartf, pEndf, f; 1139 pStart = PETSC_MAX_INT; 1140 pEnd = PETSC_MIN_INT; 1141 for (f = 0; f < patch->nsubspaces; ++f) { 1142 ierr = PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);CHKERRQ(ierr); 1143 pStart = PetscMin(pStart, pStartf); 1144 pEnd = PetscMax(pEnd, pEndf); 1145 } 1146 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1147 for (f = 0; f < patch->nsubspaces; ++f) { 1148 ierr = PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);CHKERRQ(ierr); 1149 for (p = pStartf; p < pEndf; ++p) { 1150 PetscInt fdof; 1151 ierr = PetscSectionGetDof(patch->dofSection[f], p, &fdof);CHKERRQ(ierr); 1152 ierr = PetscSectionAddDof(patch->patchSection, p, fdof);CHKERRQ(ierr); 1153 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1154 } 1155 } 1156 } 1157 ierr = PetscSectionSetUp(patch->patchSection);CHKERRQ(ierr); 1158 ierr = PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE);CHKERRQ(ierr); 1159 /* Replace cell indices with firedrake-numbered ones. */ 1160 ierr = ISGeneralSetIndices(cells, numCells, (const PetscInt *) newCellsArray, PETSC_OWN_POINTER);CHKERRQ(ierr); 1161 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol);CHKERRQ(ierr); 1162 ierr = PetscObjectSetName((PetscObject) patch->gtol, "Global Indices");CHKERRQ(ierr); 1163 ierr = PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject) pc, "-pc_patch_g2l_view");CHKERRQ(ierr); 1164 ierr = ISViewFromOptions(patch->gtol, (PetscObject) pc, "-pc_patch_g2l_view");CHKERRQ(ierr); 1165 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs);CHKERRQ(ierr); 1166 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArray, PETSC_OWN_POINTER, &patch->offs);CHKERRQ(ierr); 1167 PetscFunctionReturn(0); 1168 } 1169 1170 static PetscErrorCode PCPatchZeroFillMatrix_Private(Mat mat, const PetscInt ncell, const PetscInt ndof, const PetscInt *dof) 1171 { 1172 PetscScalar *values = NULL; 1173 PetscInt rows, c, i; 1174 PetscErrorCode ierr; 1175 1176 PetscFunctionBegin; 1177 ierr = PetscCalloc1(ndof*ndof, &values);CHKERRQ(ierr); 1178 for (c = 0; c < ncell; ++c) { 1179 const PetscInt *idx = &dof[ndof*c]; 1180 ierr = MatSetValues(mat, ndof, idx, ndof, idx, values, INSERT_VALUES);CHKERRQ(ierr); 1181 } 1182 ierr = MatGetLocalSize(mat, &rows, NULL);CHKERRQ(ierr); 1183 for (i = 0; i < rows; ++i) { 1184 ierr = MatSetValues(mat, 1, &i, 1, &i, values, INSERT_VALUES);CHKERRQ(ierr); 1185 } 1186 ierr = MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1187 ierr = MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1188 ierr = PetscFree(values);CHKERRQ(ierr); 1189 PetscFunctionReturn(0); 1190 } 1191 1192 static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat) 1193 { 1194 PC_PATCH *patch = (PC_PATCH *) pc->data; 1195 Vec x, y; 1196 PetscBool flg; 1197 PetscInt csize, rsize; 1198 const char *prefix = NULL; 1199 PetscErrorCode ierr; 1200 1201 PetscFunctionBegin; 1202 x = patch->patchX[point]; 1203 y = patch->patchY[point]; 1204 ierr = VecGetSize(x, &csize);CHKERRQ(ierr); 1205 ierr = VecGetSize(y, &rsize);CHKERRQ(ierr); 1206 ierr = MatCreate(PETSC_COMM_SELF, mat);CHKERRQ(ierr); 1207 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1208 ierr = MatSetOptionsPrefix(*mat, prefix);CHKERRQ(ierr); 1209 ierr = MatAppendOptionsPrefix(*mat, "pc_patch_sub_");CHKERRQ(ierr); 1210 if (patch->sub_mat_type) {ierr = MatSetType(*mat, patch->sub_mat_type);CHKERRQ(ierr);} 1211 else if (!patch->sub_mat_type) {ierr = MatSetType(*mat, MATDENSE);CHKERRQ(ierr);} 1212 ierr = MatSetSizes(*mat, rsize, csize, rsize, csize);CHKERRQ(ierr); 1213 ierr = PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);CHKERRQ(ierr); 1214 if (!flg) {ierr = PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);CHKERRQ(ierr);} 1215 /* Sparse patch matrices */ 1216 if (!flg) { 1217 PetscBT bt; 1218 PetscInt *dnnz = NULL; 1219 const PetscInt *dofsArray = NULL; 1220 PetscInt pStart, pEnd, ncell, offset, c, i, j; 1221 1222 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1223 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1224 point += pStart; 1225 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1226 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1227 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1228 ierr = PetscCalloc1(rsize, &dnnz);CHKERRQ(ierr); 1229 ierr = PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1230 /* XXX: This uses N^2 bits to store the sparsity pattern on a 1231 * patch. This is probably OK if the patches are not too big, 1232 * but could use quite a bit of memory for planes in 3D. 1233 * Should we switch based on the value of rsize to a 1234 * hash-table (slower, but more memory efficient) approach? */ 1235 ierr = PetscBTCreate(rsize*rsize, &bt);CHKERRQ(ierr); 1236 for (c = 0; c < ncell; ++c) { 1237 const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell; 1238 for (i = 0; i < patch->totalDofsPerCell; ++i) { 1239 const PetscInt row = idx[i]; 1240 if (row < 0) continue; 1241 for (j = 0; j < patch->totalDofsPerCell; ++j) { 1242 const PetscInt col = idx[j]; 1243 const PetscInt key = row*rsize + col; 1244 if (col < 0) continue; 1245 if (!PetscBTLookupSet(bt, key)) ++dnnz[row]; 1246 } 1247 } 1248 } 1249 ierr = PetscBTDestroy(&bt);CHKERRQ(ierr); 1250 ierr = MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);CHKERRQ(ierr); 1251 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1252 ierr = PCPatchZeroFillMatrix_Private(*mat, ncell, patch->totalDofsPerCell, &dofsArray[offset*patch->totalDofsPerCell]);CHKERRQ(ierr); 1253 ierr = PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1254 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1255 } 1256 ierr = MatSetUp(*mat);CHKERRQ(ierr); 1257 PetscFunctionReturn(0); 1258 } 1259 1260 static PetscErrorCode PCPatchComputeOperator_DMPlex_Private(PC pc, PetscInt patchNum, Mat J, IS cellIS, PetscInt n, const PetscInt *l2p, void *ctx) 1261 { 1262 PC_PATCH *patch = (PC_PATCH *) pc->data; 1263 DM dm; 1264 PetscSection s; 1265 const PetscInt *parray, *oarray; 1266 PetscInt Nf = patch->nsubspaces, Np, poff, p, f; 1267 PetscErrorCode ierr; 1268 1269 PetscFunctionBegin; 1270 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1271 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1272 /* Set offset into patch */ 1273 ierr = PetscSectionGetDof(patch->pointCounts, patchNum, &Np);CHKERRQ(ierr); 1274 ierr = PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);CHKERRQ(ierr); 1275 ierr = ISGetIndices(patch->points, &parray);CHKERRQ(ierr); 1276 ierr = ISGetIndices(patch->offs, &oarray);CHKERRQ(ierr); 1277 for (f = 0; f < Nf; ++f) { 1278 for (p = 0; p < Np; ++p) { 1279 const PetscInt point = parray[poff+p]; 1280 PetscInt dof; 1281 1282 ierr = PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);CHKERRQ(ierr); 1283 ierr = PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr); 1284 if (patch->nsubspaces == 1) {ierr = PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr);} 1285 else {ierr = PetscSectionSetOffset(patch->patchSection, point, -1);CHKERRQ(ierr);} 1286 } 1287 } 1288 ierr = ISRestoreIndices(patch->points, &parray);CHKERRQ(ierr); 1289 ierr = ISRestoreIndices(patch->offs, &oarray);CHKERRQ(ierr); 1290 if (patch->viewSection) {ierr = ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);CHKERRQ(ierr);} 1291 /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */ 1292 ierr = DMPlexComputeJacobian_Patch_Internal(pc->dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, NULL, NULL, J, J, ctx);CHKERRQ(ierr); 1293 PetscFunctionReturn(0); 1294 } 1295 1296 static PetscErrorCode PCPatchComputeOperator_Private(PC pc, Mat mat, PetscInt point) 1297 { 1298 PC_PATCH *patch = (PC_PATCH *) pc->data; 1299 const PetscInt *dofsArray; 1300 const PetscInt *cellsArray; 1301 PetscInt ncell, offset, pStart, pEnd; 1302 PetscErrorCode ierr; 1303 1304 PetscFunctionBegin; 1305 ierr = PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1306 if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n"); 1307 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1308 ierr = ISGetIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1309 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1310 1311 point += pStart; 1312 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1313 1314 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1315 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1316 if (ncell <= 0) { 1317 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1318 PetscFunctionReturn(0); 1319 } 1320 PetscStackPush("PCPatch user callback"); 1321 /* Cannot reuse the same IS because the geometry info is being cached in it */ 1322 ierr = ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);CHKERRQ(ierr); 1323 ierr = patch->usercomputeop(pc, point, mat, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, patch->usercomputectx);CHKERRQ(ierr); 1324 PetscStackPop; 1325 ierr = ISDestroy(&patch->cellIS);CHKERRQ(ierr); 1326 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1327 ierr = ISRestoreIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1328 if (patch->viewMatrix) { 1329 char name[PETSC_MAX_PATH_LEN]; 1330 1331 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);CHKERRQ(ierr); 1332 ierr = PetscObjectSetName((PetscObject) mat, name);CHKERRQ(ierr); 1333 ierr = ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);CHKERRQ(ierr); 1334 } 1335 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1336 PetscFunctionReturn(0); 1337 } 1338 1339 static PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat) 1340 { 1341 PC_PATCH *patch = (PC_PATCH *) pc->data; 1342 const PetscScalar *xArray = NULL; 1343 PetscScalar *yArray = NULL; 1344 const PetscInt *gtolArray = NULL; 1345 PetscInt dof, offset, lidx; 1346 PetscErrorCode ierr; 1347 1348 PetscFunctionBeginHot; 1349 ierr = PetscLogEventBegin(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 1350 ierr = VecGetArrayRead(x, &xArray);CHKERRQ(ierr); 1351 ierr = VecGetArray(y, &yArray);CHKERRQ(ierr); 1352 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 1353 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 1354 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1355 if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n"); 1356 if (mode == ADD_VALUES && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n"); 1357 for (lidx = 0; lidx < dof; ++lidx) { 1358 const PetscInt gidx = gtolArray[offset+lidx]; 1359 1360 if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */ 1361 else yArray[gidx] += xArray[lidx]; /* Reverse */ 1362 } 1363 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1364 ierr = VecRestoreArrayRead(x, &xArray);CHKERRQ(ierr); 1365 ierr = VecRestoreArray(y, &yArray);CHKERRQ(ierr); 1366 ierr = PetscLogEventEnd(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 1367 PetscFunctionReturn(0); 1368 } 1369 1370 static PetscErrorCode PCSetUp_PATCH(PC pc) 1371 { 1372 PC_PATCH *patch = (PC_PATCH *) pc->data; 1373 PetscInt i; 1374 const char *prefix; 1375 PetscErrorCode ierr; 1376 1377 PetscFunctionBegin; 1378 if (!pc->setupcalled) { 1379 PetscInt pStart, pEnd, p; 1380 PetscInt localSize; 1381 1382 ierr = PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 1383 1384 if (!patch->nsubspaces) { 1385 DM dm; 1386 PetscDS prob; 1387 PetscSection s; 1388 PetscInt cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, totNb = 0, **cellDofs; 1389 1390 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1391 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()"); 1392 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1393 ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr); 1394 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr); 1395 for (p = pStart; p < pEnd; ++p) { 1396 PetscInt cdof; 1397 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 1398 numGlobalBcs += cdof; 1399 } 1400 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1401 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1402 ierr = PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);CHKERRQ(ierr); 1403 for (f = 0; f < Nf; ++f) { 1404 PetscFE fe; 1405 PetscDualSpace sp; 1406 PetscInt cdoff = 0; 1407 1408 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1409 /* ierr = PetscFEGetNumComponents(fe, &Nc[f]);CHKERRQ(ierr); */ 1410 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 1411 ierr = PetscDualSpaceGetDimension(sp, &Nb[f]);CHKERRQ(ierr); 1412 totNb += Nb[f]; 1413 1414 ierr = PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);CHKERRQ(ierr); 1415 for (c = cStart; c < cEnd; ++c) { 1416 PetscInt *closure = NULL; 1417 PetscInt clSize = 0, cl; 1418 1419 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1420 for (cl = 0; cl < clSize*2; cl += 2) { 1421 const PetscInt p = closure[cl]; 1422 PetscInt fdof, d, foff; 1423 1424 ierr = PetscSectionGetFieldDof(s, p, f, &fdof);CHKERRQ(ierr); 1425 ierr = PetscSectionGetFieldOffset(s, p, f, &foff);CHKERRQ(ierr); 1426 for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d; 1427 } 1428 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1429 } 1430 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]); 1431 } 1432 numGlobalBcs = 0; 1433 for (p = pStart; p < pEnd; ++p) { 1434 const PetscInt *ind; 1435 PetscInt off, cdof, d; 1436 1437 ierr = PetscSectionGetOffset(s, p, &off);CHKERRQ(ierr); 1438 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 1439 ierr = PetscSectionGetConstraintIndices(s, p, &ind);CHKERRQ(ierr); 1440 for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d]; 1441 } 1442 1443 ierr = PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);CHKERRQ(ierr); 1444 for (f = 0; f < Nf; ++f) { 1445 ierr = PetscFree(cellDofs[f]);CHKERRQ(ierr); 1446 } 1447 ierr = PetscFree3(Nb, cellDofs, globalBcs);CHKERRQ(ierr); 1448 ierr = PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);CHKERRQ(ierr); 1449 } 1450 1451 localSize = patch->subspaceOffsets[patch->nsubspaces]; 1452 ierr = VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localX);CHKERRQ(ierr); 1453 ierr = VecSetUp(patch->localX);CHKERRQ(ierr); 1454 ierr = VecDuplicate(patch->localX, &patch->localY);CHKERRQ(ierr); 1455 ierr = PCPatchCreateCellPatches(pc);CHKERRQ(ierr); 1456 ierr = PCPatchCreateCellPatchDiscretisationInfo(pc);CHKERRQ(ierr); 1457 1458 /* OK, now build the work vectors */ 1459 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd);CHKERRQ(ierr); 1460 ierr = PetscMalloc1(patch->npatch, &patch->patchX);CHKERRQ(ierr); 1461 ierr = PetscMalloc1(patch->npatch, &patch->patchY);CHKERRQ(ierr); 1462 for (p = pStart; p < pEnd; ++p) { 1463 PetscInt dof; 1464 1465 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 1466 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchX[p-pStart]);CHKERRQ(ierr); 1467 ierr = VecSetUp(patch->patchX[p-pStart]);CHKERRQ(ierr); 1468 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchY[p-pStart]);CHKERRQ(ierr); 1469 ierr = VecSetUp(patch->patchY[p-pStart]);CHKERRQ(ierr); 1470 } 1471 ierr = PetscMalloc1(patch->npatch, &patch->ksp);CHKERRQ(ierr); 1472 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1473 for (i = 0; i < patch->npatch; ++i) { 1474 PC subpc; 1475 1476 ierr = KSPCreate(PETSC_COMM_SELF, &patch->ksp[i]);CHKERRQ(ierr); 1477 ierr = KSPSetOptionsPrefix(patch->ksp[i], prefix);CHKERRQ(ierr); 1478 ierr = KSPAppendOptionsPrefix(patch->ksp[i], "sub_");CHKERRQ(ierr); 1479 ierr = PetscObjectIncrementTabLevel((PetscObject) patch->ksp[i], (PetscObject) pc, 1);CHKERRQ(ierr); 1480 ierr = KSPGetPC(patch->ksp[i], &subpc);CHKERRQ(ierr); 1481 ierr = PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);CHKERRQ(ierr); 1482 ierr = PetscLogObjectParent((PetscObject) pc, (PetscObject) patch->ksp[i]);CHKERRQ(ierr); 1483 } 1484 if (patch->save_operators) { 1485 ierr = PetscMalloc1(patch->npatch, &patch->mat);CHKERRQ(ierr); 1486 for (i = 0; i < patch->npatch; ++i) { 1487 ierr = PCPatchCreateMatrix_Private(pc, i, &patch->mat[i]);CHKERRQ(ierr); 1488 } 1489 } 1490 ierr = PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 1491 1492 /* If desired, calculate weights for dof multiplicity */ 1493 if (patch->partition_of_unity) { 1494 ierr = VecDuplicate(patch->localX, &patch->dof_weights);CHKERRQ(ierr); 1495 for (i = 0; i < patch->npatch; ++i) { 1496 PetscInt dof; 1497 1498 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);CHKERRQ(ierr); 1499 if (dof <= 0) continue; 1500 ierr = VecSet(patch->patchX[i], 1.0);CHKERRQ(ierr); 1501 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchX[i], patch->dof_weights, ADD_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 1502 } 1503 1504 PetscScalar *input = NULL; 1505 PetscScalar *output = NULL; 1506 Vec global; 1507 VecDuplicate(patch->dof_weights, &global); 1508 VecSet(global, 0.); 1509 1510 ierr = VecGetArray(patch->dof_weights, &input);CHKERRQ(ierr); 1511 ierr = VecGetArray(global, &output);CHKERRQ(ierr); 1512 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);CHKERRQ(ierr); 1513 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);CHKERRQ(ierr); 1514 ierr = VecRestoreArray(patch->dof_weights, &input);CHKERRQ(ierr); 1515 ierr = VecRestoreArray(global, &output);CHKERRQ(ierr); 1516 1517 ierr = VecReciprocal(patch->dof_weights);CHKERRQ(ierr); 1518 1519 ierr = VecGetArray(patch->dof_weights, &output);CHKERRQ(ierr); 1520 ierr = VecGetArray(global, &input);CHKERRQ(ierr); 1521 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, input, output);CHKERRQ(ierr); 1522 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, input, output);CHKERRQ(ierr); 1523 ierr = VecRestoreArray(patch->dof_weights, &output);CHKERRQ(ierr); 1524 ierr = VecRestoreArray(global, &input);CHKERRQ(ierr); 1525 ierr = VecDestroy(&global);CHKERRQ(ierr); 1526 } 1527 } 1528 if (patch->save_operators) { 1529 for (i = 0; i < patch->npatch; ++i) { 1530 ierr = MatZeroEntries(patch->mat[i]);CHKERRQ(ierr); 1531 ierr = PCPatchComputeOperator_Private(pc, patch->mat[i], i);CHKERRQ(ierr); 1532 ierr = KSPSetOperators(patch->ksp[i], patch->mat[i], patch->mat[i]);CHKERRQ(ierr); 1533 } 1534 } 1535 PetscFunctionReturn(0); 1536 } 1537 1538 static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y) 1539 { 1540 PC_PATCH *patch = (PC_PATCH *) pc->data; 1541 const PetscScalar *globalX = NULL; 1542 PetscScalar *localX = NULL; 1543 PetscScalar *globalY = NULL; 1544 const PetscInt *bcNodes = NULL; 1545 PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1; 1546 PetscInt start[2] = {0, 0}; 1547 PetscInt end[2] = {-1, -1}; 1548 const PetscInt inc[2] = {1, -1}; 1549 const PetscScalar *localY; 1550 const PetscInt *iterationSet; 1551 PetscInt pStart, numBcs, n, sweep, bc, j; 1552 PetscErrorCode ierr; 1553 1554 PetscFunctionBegin; 1555 ierr = PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 1556 ierr = PetscOptionsPushGetViewerOff(PETSC_TRUE);CHKERRQ(ierr); 1557 end[0] = patch->npatch; 1558 start[1] = patch->npatch-1; 1559 if (patch->user_patches) { 1560 ierr = ISGetLocalSize(patch->iterationSet, &end[0]);CHKERRQ(ierr); 1561 start[1] = end[0] - 1; 1562 ierr = ISGetIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr); 1563 } 1564 /* Scatter from global space into overlapped local spaces */ 1565 ierr = VecGetArrayRead(x, &globalX);CHKERRQ(ierr); 1566 ierr = VecGetArray(patch->localX, &localX);CHKERRQ(ierr); 1567 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, globalX, localX);CHKERRQ(ierr); 1568 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, globalX, localX);CHKERRQ(ierr); 1569 ierr = VecRestoreArrayRead(x, &globalX);CHKERRQ(ierr); 1570 ierr = VecRestoreArray(patch->localX, &localX);CHKERRQ(ierr); 1571 1572 ierr = VecSet(patch->localY, 0.0);CHKERRQ(ierr); 1573 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);CHKERRQ(ierr); 1574 for (sweep = 0; sweep < nsweep; sweep++) { 1575 for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) { 1576 PetscInt i = patch->user_patches ? iterationSet[j] : j; 1577 PetscInt start, len; 1578 1579 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);CHKERRQ(ierr); 1580 ierr = PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);CHKERRQ(ierr); 1581 /* TODO: Squash out these guys in the setup as well. */ 1582 if (len <= 0) continue; 1583 /* TODO: Do we need different scatters for X and Y? */ 1584 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localX, patch->patchX[i], INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1585 if (!patch->save_operators) { 1586 Mat mat; 1587 1588 ierr = PCPatchCreateMatrix_Private(pc, i, &mat);CHKERRQ(ierr); 1589 /* Populate operator here. */ 1590 ierr = PCPatchComputeOperator_Private(pc, mat, i);CHKERRQ(ierr); 1591 ierr = KSPSetOperators(patch->ksp[i], mat, mat);CHKERRQ(ierr); 1592 /* Drop reference so the KSPSetOperators below will blow it away. */ 1593 ierr = MatDestroy(&mat);CHKERRQ(ierr); 1594 } 1595 ierr = PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 1596 if (!patch->ksp[i]->setfromoptionscalled) { 1597 ierr = KSPSetFromOptions(patch->ksp[i]);CHKERRQ(ierr); 1598 } 1599 ierr = KSPSolve(patch->ksp[i], patch->patchX[i], patch->patchY[i]);CHKERRQ(ierr); 1600 ierr = PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 1601 1602 if (!patch->save_operators) { 1603 PC pc; 1604 ierr = KSPSetOperators(patch->ksp[i], NULL, NULL);CHKERRQ(ierr); 1605 ierr = KSPGetPC(patch->ksp[i], &pc);CHKERRQ(ierr); 1606 /* Destroy PC context too, otherwise the factored matrix hangs around. */ 1607 ierr = PCReset(pc);CHKERRQ(ierr); 1608 } 1609 1610 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchY[i], patch->localY, ADD_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 1611 } 1612 } 1613 if (patch->user_patches) {ierr = ISRestoreIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr);} 1614 /* XXX: should we do this on the global vector? */ 1615 if (patch->partition_of_unity) { 1616 ierr = VecPointwiseMult(patch->localY, patch->localY, patch->dof_weights);CHKERRQ(ierr); 1617 } 1618 /* Now patch->localY contains the solution of the patch solves, so we need to combine them all. */ 1619 ierr = VecSet(y, 0.0);CHKERRQ(ierr); 1620 ierr = VecGetArray(y, &globalY);CHKERRQ(ierr); 1621 ierr = VecGetArrayRead(patch->localY, &localY);CHKERRQ(ierr); 1622 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, localY, globalY, MPI_SUM);CHKERRQ(ierr); 1623 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, localY, globalY, MPI_SUM);CHKERRQ(ierr); 1624 ierr = VecRestoreArrayRead(patch->localY, &localY);CHKERRQ(ierr); 1625 1626 /* Now we need to send the global BC values through */ 1627 ierr = VecGetArrayRead(x, &globalX);CHKERRQ(ierr); 1628 ierr = ISGetSize(patch->globalBcNodes, &numBcs);CHKERRQ(ierr); 1629 ierr = ISGetIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 1630 ierr = VecGetLocalSize(x, &n);CHKERRQ(ierr); 1631 for (bc = 0; bc < numBcs; ++bc) { 1632 const PetscInt idx = bcNodes[bc]; 1633 if (idx < n) globalY[idx] = globalX[idx]; 1634 } 1635 1636 ierr = ISRestoreIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 1637 ierr = VecRestoreArrayRead(x, &globalX);CHKERRQ(ierr); 1638 ierr = VecRestoreArray(y, &globalY);CHKERRQ(ierr); 1639 1640 ierr = PetscOptionsPopGetViewerOff();CHKERRQ(ierr); 1641 ierr = PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 1642 PetscFunctionReturn(0); 1643 } 1644 1645 static PetscErrorCode PCReset_PATCH(PC pc) 1646 { 1647 PC_PATCH *patch = (PC_PATCH *) pc->data; 1648 PetscInt i; 1649 PetscErrorCode ierr; 1650 1651 PetscFunctionBegin; 1652 /* TODO: Get rid of all these ifs */ 1653 ierr = PetscSFDestroy(&patch->defaultSF);CHKERRQ(ierr); 1654 ierr = PetscSectionDestroy(&patch->cellCounts);CHKERRQ(ierr); 1655 ierr = PetscSectionDestroy(&patch->pointCounts);CHKERRQ(ierr); 1656 ierr = PetscSectionDestroy(&patch->cellNumbering);CHKERRQ(ierr); 1657 ierr = PetscSectionDestroy(&patch->gtolCounts);CHKERRQ(ierr); 1658 ierr = ISDestroy(&patch->gtol);CHKERRQ(ierr); 1659 ierr = ISDestroy(&patch->cells);CHKERRQ(ierr); 1660 ierr = ISDestroy(&patch->points);CHKERRQ(ierr); 1661 ierr = ISDestroy(&patch->dofs);CHKERRQ(ierr); 1662 ierr = ISDestroy(&patch->offs);CHKERRQ(ierr); 1663 ierr = PetscSectionDestroy(&patch->patchSection);CHKERRQ(ierr); 1664 ierr = ISDestroy(&patch->ghostBcNodes);CHKERRQ(ierr); 1665 ierr = ISDestroy(&patch->globalBcNodes);CHKERRQ(ierr); 1666 1667 if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscSectionDestroy(&patch->dofSection[i]);CHKERRQ(ierr);} 1668 ierr = PetscFree(patch->dofSection);CHKERRQ(ierr); 1669 ierr = PetscFree(patch->bs);CHKERRQ(ierr); 1670 ierr = PetscFree(patch->nodesPerCell);CHKERRQ(ierr); 1671 if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscFree(patch->cellNodeMap[i]);CHKERRQ(ierr);} 1672 ierr = PetscFree(patch->cellNodeMap);CHKERRQ(ierr); 1673 ierr = PetscFree(patch->subspaceOffsets);CHKERRQ(ierr); 1674 1675 if (patch->ksp) { 1676 for (i = 0; i < patch->npatch; ++i) {ierr = KSPReset(patch->ksp[i]);CHKERRQ(ierr);} 1677 } 1678 1679 if (patch->subspaces_to_exclude) { 1680 PetscHSetIDestroy(&patch->subspaces_to_exclude); 1681 } 1682 1683 ierr = VecDestroy(&patch->localX);CHKERRQ(ierr); 1684 ierr = VecDestroy(&patch->localY);CHKERRQ(ierr); 1685 if (patch->patchX) { 1686 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchX[i]);CHKERRQ(ierr);} 1687 ierr = PetscFree(patch->patchX);CHKERRQ(ierr); 1688 } 1689 if (patch->patchY) { 1690 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchY[i]);CHKERRQ(ierr);} 1691 ierr = PetscFree(patch->patchY);CHKERRQ(ierr); 1692 } 1693 ierr = VecDestroy(&patch->dof_weights);CHKERRQ(ierr); 1694 if (patch->patch_dof_weights) { 1695 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patch_dof_weights[i]);CHKERRQ(ierr);} 1696 ierr = PetscFree(patch->patch_dof_weights);CHKERRQ(ierr); 1697 } 1698 if (patch->mat) { 1699 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->mat[i]);CHKERRQ(ierr);} 1700 ierr = PetscFree(patch->mat);CHKERRQ(ierr); 1701 } 1702 ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr); 1703 if (patch->userIS) { 1704 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->userIS[i]);CHKERRQ(ierr);} 1705 ierr = PetscFree(patch->userIS);CHKERRQ(ierr); 1706 } 1707 patch->bs = 0; 1708 patch->cellNodeMap = NULL; 1709 patch->nsubspaces = 0; 1710 ierr = ISDestroy(&patch->iterationSet);CHKERRQ(ierr); 1711 1712 ierr = PetscViewerDestroy(&patch->viewerSection);CHKERRQ(ierr); 1713 PetscFunctionReturn(0); 1714 } 1715 1716 static PetscErrorCode PCDestroy_PATCH(PC pc) 1717 { 1718 PC_PATCH *patch = (PC_PATCH *) pc->data; 1719 PetscInt i; 1720 PetscErrorCode ierr; 1721 1722 PetscFunctionBegin; 1723 ierr = PCReset_PATCH(pc);CHKERRQ(ierr); 1724 if (patch->ksp) { 1725 for (i = 0; i < patch->npatch; ++i) {ierr = KSPDestroy(&patch->ksp[i]);CHKERRQ(ierr);} 1726 ierr = PetscFree(patch->ksp);CHKERRQ(ierr); 1727 } 1728 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1729 PetscFunctionReturn(0); 1730 } 1731 1732 static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc) 1733 { 1734 PC_PATCH *patch = (PC_PATCH *) pc->data; 1735 PCPatchConstructType patchConstructionType = PC_PATCH_STAR; 1736 char sub_mat_type[PETSC_MAX_PATH_LEN]; 1737 const char *prefix; 1738 PetscBool flg, dimflg, codimflg; 1739 MPI_Comm comm; 1740 PetscInt *ifields, nfields, k; 1741 PetscErrorCode ierr; 1742 1743 PetscFunctionBegin; 1744 ierr = PetscObjectGetComm((PetscObject) pc, &comm);CHKERRQ(ierr); 1745 ierr = PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);CHKERRQ(ierr); 1746 ierr = PetscOptionsHead(PetscOptionsObject, "Vertex-patch Additive Schwarz options");CHKERRQ(ierr); 1747 ierr = PetscOptionsBool("-pc_patch_save_operators", "Store all patch operators for lifetime of PC?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg);CHKERRQ(ierr); 1748 ierr = PetscOptionsBool("-pc_patch_partition_of_unity", "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg);CHKERRQ(ierr); 1749 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); 1750 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); 1751 if (dimflg && codimflg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");CHKERRQ(ierr); 1752 ierr = PetscOptionsEnum("-pc_patch_construct_type", "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum) patchConstructionType, (PetscEnum *) &patchConstructionType, &flg);CHKERRQ(ierr); 1753 if (flg) {ierr = PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL);CHKERRQ(ierr);} 1754 ierr = PetscOptionsInt("-pc_patch_vanka_dim", "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg);CHKERRQ(ierr); 1755 ierr = PetscOptionsInt("-pc_patch_ignore_dim", "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg);CHKERRQ(ierr); 1756 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); 1757 if (flg) {ierr = PCPatchSetSubMatType(pc, sub_mat_type);CHKERRQ(ierr);} 1758 ierr = PetscOptionsBool("-pc_patch_symmetrise_sweep", "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg);CHKERRQ(ierr); 1759 1760 /* If the user has set the number of subspaces, use that for the buffer size, 1761 otherwise use a large number */ 1762 if (patch->nsubspaces <= 0) { 1763 nfields = 128; 1764 } else { 1765 nfields = patch->nsubspaces; 1766 } 1767 ierr = PetscMalloc1(nfields, &ifields);CHKERRQ(ierr); 1768 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_patch_exclude_subspaces",ifields,&nfields,&flg);CHKERRQ(ierr); 1769 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"); 1770 if (flg) { 1771 PetscHSetIClear(patch->subspaces_to_exclude); 1772 for (k = 0; k < nfields; k++) { 1773 PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]); 1774 } 1775 } 1776 ierr = PetscFree(ifields);CHKERRQ(ierr); 1777 1778 ierr = PetscOptionsBool("-pc_patch_patches_view", "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);CHKERRQ(ierr); 1779 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, "-pc_patch_cells_view", &patch->viewerCells, &patch->formatCells, &patch->viewCells);CHKERRQ(ierr); 1780 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, "-pc_patch_points_view", &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);CHKERRQ(ierr); 1781 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, "-pc_patch_section_view", &patch->viewerSection, &patch->formatSection, &patch->viewSection);CHKERRQ(ierr); 1782 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, "-pc_patch_mat_view", &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);CHKERRQ(ierr); 1783 ierr = PetscOptionsTail();CHKERRQ(ierr); 1784 patch->optionsSet = PETSC_TRUE; 1785 PetscFunctionReturn(0); 1786 } 1787 1788 static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc) 1789 { 1790 PC_PATCH *patch = (PC_PATCH*) pc->data; 1791 KSPConvergedReason reason; 1792 PetscInt i; 1793 PetscErrorCode ierr; 1794 1795 PetscFunctionBegin; 1796 if (!patch->save_operators) { 1797 /* Can't do this here because the sub KSPs don't have an operator attached yet. */ 1798 PetscFunctionReturn(0); 1799 } 1800 for (i = 0; i < patch->npatch; ++i) { 1801 if (!patch->ksp[i]->setfromoptionscalled) { 1802 ierr = KSPSetFromOptions(patch->ksp[i]);CHKERRQ(ierr); 1803 } 1804 ierr = KSPSetUp(patch->ksp[i]);CHKERRQ(ierr); 1805 ierr = KSPGetConvergedReason(patch->ksp[i], &reason);CHKERRQ(ierr); 1806 if (reason == KSP_DIVERGED_PCSETUP_FAILED) pc->failedreason = PC_SUBPC_ERROR; 1807 } 1808 PetscFunctionReturn(0); 1809 } 1810 1811 static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer) 1812 { 1813 PC_PATCH *patch = (PC_PATCH *) pc->data; 1814 PetscViewer sviewer; 1815 PetscBool isascii; 1816 PetscMPIInt rank; 1817 PetscErrorCode ierr; 1818 1819 PetscFunctionBegin; 1820 /* TODO Redo tabbing with set tbas in new style */ 1821 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 1822 if (!isascii) PetscFunctionReturn(0); 1823 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);CHKERRQ(ierr); 1824 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1825 ierr = PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);CHKERRQ(ierr); 1826 ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");CHKERRQ(ierr); 1827 if (patch->partition_of_unity) {ierr = PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");CHKERRQ(ierr);} 1828 else {ierr = PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");CHKERRQ(ierr);} 1829 if (patch->symmetrise_sweep) {ierr = PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");CHKERRQ(ierr);} 1830 else {ierr = PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");CHKERRQ(ierr);} 1831 if (!patch->save_operators) {ierr = PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");CHKERRQ(ierr);} 1832 else {ierr = PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");CHKERRQ(ierr);} 1833 if (patch->patchconstructop == PCPatchConstruct_Star) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");CHKERRQ(ierr);} 1834 else if (patch->patchconstructop == PCPatchConstruct_Vanka) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");CHKERRQ(ierr);} 1835 else if (patch->patchconstructop == PCPatchConstruct_User) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");CHKERRQ(ierr);} 1836 else {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");CHKERRQ(ierr);} 1837 ierr = PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");CHKERRQ(ierr); 1838 if (patch->ksp) { 1839 ierr = PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 1840 if (!rank) { 1841 ierr = PetscViewerASCIIPushTab(sviewer);CHKERRQ(ierr); 1842 ierr = KSPView(patch->ksp[0], sviewer);CHKERRQ(ierr); 1843 ierr = PetscViewerASCIIPopTab(sviewer);CHKERRQ(ierr); 1844 } 1845 ierr = PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 1846 } else { 1847 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1848 ierr = PetscViewerASCIIPrintf(viewer, "KSP not yet set.\n");CHKERRQ(ierr); 1849 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1850 } 1851 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1852 PetscFunctionReturn(0); 1853 } 1854 1855 /*MC 1856 PCPATCH - A PC object that encapsulates flexible definition of blocks for overlapping and non-overlapping 1857 small block additive preconditioners. Block definition is based on topology from 1858 a DM and equation numbering from a PetscSection. 1859 1860 Options Database Keys: 1861 + -pc_patch_cells_view - Views the process local cell numbers for each patch 1862 . -pc_patch_points_view - Views the process local mesh point numbers for each patch 1863 . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch 1864 . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary 1865 - -pc_patch_sub_mat_view - Views the matrix associated with each patch 1866 1867 Level: intermediate 1868 1869 .seealso: PCType, PCCreate(), PCSetType() 1870 M*/ 1871 PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc) 1872 { 1873 PC_PATCH *patch; 1874 PetscErrorCode ierr; 1875 1876 PetscFunctionBegin; 1877 ierr = PetscNewLog(pc, &patch);CHKERRQ(ierr); 1878 1879 if (patch->subspaces_to_exclude) { 1880 PetscHSetIDestroy(&patch->subspaces_to_exclude); 1881 } 1882 PetscHSetICreate(&patch->subspaces_to_exclude); 1883 1884 /* Set some defaults */ 1885 patch->combined = PETSC_FALSE; 1886 patch->save_operators = PETSC_TRUE; 1887 patch->partition_of_unity = PETSC_FALSE; 1888 patch->codim = -1; 1889 patch->dim = -1; 1890 patch->vankadim = -1; 1891 patch->ignoredim = -1; 1892 patch->patchconstructop = PCPatchConstruct_Star; 1893 patch->symmetrise_sweep = PETSC_FALSE; 1894 patch->npatch = 0; 1895 patch->userIS = NULL; 1896 patch->optionsSet = PETSC_FALSE; 1897 patch->iterationSet = NULL; 1898 patch->user_patches = PETSC_FALSE; 1899 ierr = PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 1900 patch->viewPatches = PETSC_FALSE; 1901 patch->viewCells = PETSC_FALSE; 1902 patch->viewPoints = PETSC_FALSE; 1903 patch->viewSection = PETSC_FALSE; 1904 patch->viewMatrix = PETSC_FALSE; 1905 1906 pc->data = (void *) patch; 1907 pc->ops->apply = PCApply_PATCH; 1908 pc->ops->applytranspose = 0; /* PCApplyTranspose_PATCH; */ 1909 pc->ops->setup = PCSetUp_PATCH; 1910 pc->ops->reset = PCReset_PATCH; 1911 pc->ops->destroy = PCDestroy_PATCH; 1912 pc->ops->setfromoptions = PCSetFromOptions_PATCH; 1913 pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH; 1914 pc->ops->view = PCView_PATCH; 1915 pc->ops->applyrichardson = 0; 1916 1917 PetscFunctionReturn(0); 1918 } 1919