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], &patch->cellNodeMap[i]);CHKERRQ(ierr); 410 for (j = 0; j < (cEnd-cStart)*nodesPerCell[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 ierr = VecReciprocal(patch->dof_weights);CHKERRQ(ierr); 1504 } 1505 } 1506 if (patch->save_operators) { 1507 for (i = 0; i < patch->npatch; ++i) { 1508 ierr = MatZeroEntries(patch->mat[i]);CHKERRQ(ierr); 1509 ierr = PCPatchComputeOperator_Private(pc, patch->mat[i], i);CHKERRQ(ierr); 1510 ierr = KSPSetOperators(patch->ksp[i], patch->mat[i], patch->mat[i]);CHKERRQ(ierr); 1511 } 1512 } 1513 PetscFunctionReturn(0); 1514 } 1515 1516 static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y) 1517 { 1518 PC_PATCH *patch = (PC_PATCH *) pc->data; 1519 const PetscScalar *globalX = NULL; 1520 PetscScalar *localX = NULL; 1521 PetscScalar *globalY = NULL; 1522 const PetscInt *bcNodes = NULL; 1523 PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1; 1524 PetscInt start[2] = {0, 0}; 1525 PetscInt end[2] = {-1, -1}; 1526 const PetscInt inc[2] = {1, -1}; 1527 const PetscScalar *localY; 1528 const PetscInt *iterationSet; 1529 PetscInt pStart, numBcs, n, sweep, bc, j; 1530 PetscErrorCode ierr; 1531 1532 PetscFunctionBegin; 1533 ierr = PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 1534 ierr = PetscOptionsPushGetViewerOff(PETSC_TRUE);CHKERRQ(ierr); 1535 end[0] = patch->npatch; 1536 start[1] = patch->npatch-1; 1537 if (patch->user_patches) { 1538 ierr = ISGetLocalSize(patch->iterationSet, &end[0]);CHKERRQ(ierr); 1539 start[1] = end[0] - 1; 1540 ierr = ISGetIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr); 1541 } 1542 /* Scatter from global space into overlapped local spaces */ 1543 ierr = VecGetArrayRead(x, &globalX);CHKERRQ(ierr); 1544 ierr = VecGetArray(patch->localX, &localX);CHKERRQ(ierr); 1545 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, globalX, localX);CHKERRQ(ierr); 1546 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, globalX, localX);CHKERRQ(ierr); 1547 ierr = VecRestoreArrayRead(x, &globalX);CHKERRQ(ierr); 1548 ierr = VecRestoreArray(patch->localX, &localX);CHKERRQ(ierr); 1549 1550 ierr = VecSet(patch->localY, 0.0);CHKERRQ(ierr); 1551 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);CHKERRQ(ierr); 1552 for (sweep = 0; sweep < nsweep; sweep++) { 1553 for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) { 1554 PetscInt i = patch->user_patches ? iterationSet[j] : j; 1555 PetscInt start, len; 1556 1557 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);CHKERRQ(ierr); 1558 ierr = PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);CHKERRQ(ierr); 1559 /* TODO: Squash out these guys in the setup as well. */ 1560 if (len <= 0) continue; 1561 /* TODO: Do we need different scatters for X and Y? */ 1562 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localX, patch->patchX[i], INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 1563 if (!patch->save_operators) { 1564 Mat mat; 1565 1566 ierr = PCPatchCreateMatrix_Private(pc, i, &mat);CHKERRQ(ierr); 1567 /* Populate operator here. */ 1568 ierr = PCPatchComputeOperator_Private(pc, mat, i);CHKERRQ(ierr); 1569 ierr = KSPSetOperators(patch->ksp[i], mat, mat);CHKERRQ(ierr); 1570 /* Drop reference so the KSPSetOperators below will blow it away. */ 1571 ierr = MatDestroy(&mat);CHKERRQ(ierr); 1572 } 1573 ierr = PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 1574 if (!patch->ksp[i]->setfromoptionscalled) { 1575 ierr = KSPSetFromOptions(patch->ksp[i]);CHKERRQ(ierr); 1576 } 1577 ierr = KSPSolve(patch->ksp[i], patch->patchX[i], patch->patchY[i]);CHKERRQ(ierr); 1578 ierr = PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 1579 1580 if (!patch->save_operators) { 1581 PC pc; 1582 ierr = KSPSetOperators(patch->ksp[i], NULL, NULL);CHKERRQ(ierr); 1583 ierr = KSPGetPC(patch->ksp[i], &pc);CHKERRQ(ierr); 1584 /* Destroy PC context too, otherwise the factored matrix hangs around. */ 1585 ierr = PCReset(pc);CHKERRQ(ierr); 1586 } 1587 1588 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchY[i], patch->localY, ADD_VALUES, SCATTER_REVERSE);CHKERRQ(ierr); 1589 } 1590 } 1591 if (patch->user_patches) {ierr = ISRestoreIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr);} 1592 /* XXX: should we do this on the global vector? */ 1593 if (patch->partition_of_unity) { 1594 ierr = VecPointwiseMult(patch->localY, patch->localY, patch->dof_weights);CHKERRQ(ierr); 1595 } 1596 /* Now patch->localY contains the solution of the patch solves, so we need to combine them all. */ 1597 ierr = VecSet(y, 0.0);CHKERRQ(ierr); 1598 ierr = VecGetArray(y, &globalY);CHKERRQ(ierr); 1599 ierr = VecGetArrayRead(patch->localY, &localY);CHKERRQ(ierr); 1600 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, localY, globalY, MPI_SUM);CHKERRQ(ierr); 1601 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, localY, globalY, MPI_SUM);CHKERRQ(ierr); 1602 ierr = VecRestoreArrayRead(patch->localY, &localY);CHKERRQ(ierr); 1603 1604 /* Now we need to send the global BC values through */ 1605 ierr = VecGetArrayRead(x, &globalX);CHKERRQ(ierr); 1606 ierr = ISGetSize(patch->globalBcNodes, &numBcs);CHKERRQ(ierr); 1607 ierr = ISGetIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 1608 ierr = VecGetLocalSize(x, &n);CHKERRQ(ierr); 1609 for (bc = 0; bc < numBcs; ++bc) { 1610 const PetscInt idx = bcNodes[bc]; 1611 if (idx < n) globalY[idx] = globalX[idx]; 1612 } 1613 1614 ierr = ISRestoreIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 1615 ierr = VecRestoreArrayRead(x, &globalX);CHKERRQ(ierr); 1616 ierr = VecRestoreArray(y, &globalY);CHKERRQ(ierr); 1617 1618 ierr = PetscOptionsPopGetViewerOff();CHKERRQ(ierr); 1619 ierr = PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 1620 PetscFunctionReturn(0); 1621 } 1622 1623 static PetscErrorCode PCReset_PATCH(PC pc) 1624 { 1625 PC_PATCH *patch = (PC_PATCH *) pc->data; 1626 PetscInt i; 1627 PetscErrorCode ierr; 1628 1629 PetscFunctionBegin; 1630 /* TODO: Get rid of all these ifs */ 1631 ierr = PetscSFDestroy(&patch->defaultSF);CHKERRQ(ierr); 1632 ierr = PetscSectionDestroy(&patch->cellCounts);CHKERRQ(ierr); 1633 ierr = PetscSectionDestroy(&patch->pointCounts);CHKERRQ(ierr); 1634 ierr = PetscSectionDestroy(&patch->cellNumbering);CHKERRQ(ierr); 1635 ierr = PetscSectionDestroy(&patch->gtolCounts);CHKERRQ(ierr); 1636 ierr = ISDestroy(&patch->gtol);CHKERRQ(ierr); 1637 ierr = ISDestroy(&patch->cells);CHKERRQ(ierr); 1638 ierr = ISDestroy(&patch->points);CHKERRQ(ierr); 1639 ierr = ISDestroy(&patch->dofs);CHKERRQ(ierr); 1640 ierr = ISDestroy(&patch->offs);CHKERRQ(ierr); 1641 ierr = PetscSectionDestroy(&patch->patchSection);CHKERRQ(ierr); 1642 ierr = ISDestroy(&patch->ghostBcNodes);CHKERRQ(ierr); 1643 ierr = ISDestroy(&patch->globalBcNodes);CHKERRQ(ierr); 1644 1645 if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscSectionDestroy(&patch->dofSection[i]);CHKERRQ(ierr);} 1646 ierr = PetscFree(patch->dofSection);CHKERRQ(ierr); 1647 ierr = PetscFree(patch->bs);CHKERRQ(ierr); 1648 ierr = PetscFree(patch->nodesPerCell);CHKERRQ(ierr); 1649 if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscFree(patch->cellNodeMap[i]);CHKERRQ(ierr);} 1650 ierr = PetscFree(patch->cellNodeMap);CHKERRQ(ierr); 1651 ierr = PetscFree(patch->subspaceOffsets);CHKERRQ(ierr); 1652 1653 if (patch->ksp) { 1654 for (i = 0; i < patch->npatch; ++i) {ierr = KSPReset(patch->ksp[i]);CHKERRQ(ierr);} 1655 } 1656 1657 if (patch->subspaces_to_exclude) { 1658 PetscHSetIDestroy(&patch->subspaces_to_exclude); 1659 } 1660 1661 ierr = VecDestroy(&patch->localX);CHKERRQ(ierr); 1662 ierr = VecDestroy(&patch->localY);CHKERRQ(ierr); 1663 if (patch->patchX) { 1664 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchX[i]);CHKERRQ(ierr);} 1665 ierr = PetscFree(patch->patchX);CHKERRQ(ierr); 1666 } 1667 if (patch->patchY) { 1668 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchY[i]);CHKERRQ(ierr);} 1669 ierr = PetscFree(patch->patchY);CHKERRQ(ierr); 1670 } 1671 ierr = VecDestroy(&patch->dof_weights);CHKERRQ(ierr); 1672 if (patch->patch_dof_weights) { 1673 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patch_dof_weights[i]);CHKERRQ(ierr);} 1674 ierr = PetscFree(patch->patch_dof_weights);CHKERRQ(ierr); 1675 } 1676 if (patch->mat) { 1677 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->mat[i]);CHKERRQ(ierr);} 1678 ierr = PetscFree(patch->mat);CHKERRQ(ierr); 1679 } 1680 ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr); 1681 if (patch->userIS) { 1682 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->userIS[i]);CHKERRQ(ierr);} 1683 ierr = PetscFree(patch->userIS);CHKERRQ(ierr); 1684 } 1685 patch->bs = 0; 1686 patch->cellNodeMap = NULL; 1687 patch->nsubspaces = 0; 1688 ierr = ISDestroy(&patch->iterationSet);CHKERRQ(ierr); 1689 1690 ierr = PetscViewerDestroy(&patch->viewerSection);CHKERRQ(ierr); 1691 PetscFunctionReturn(0); 1692 } 1693 1694 static PetscErrorCode PCDestroy_PATCH(PC pc) 1695 { 1696 PC_PATCH *patch = (PC_PATCH *) pc->data; 1697 PetscInt i; 1698 PetscErrorCode ierr; 1699 1700 PetscFunctionBegin; 1701 ierr = PCReset_PATCH(pc);CHKERRQ(ierr); 1702 if (patch->ksp) { 1703 for (i = 0; i < patch->npatch; ++i) {ierr = KSPDestroy(&patch->ksp[i]);CHKERRQ(ierr);} 1704 ierr = PetscFree(patch->ksp);CHKERRQ(ierr); 1705 } 1706 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1707 PetscFunctionReturn(0); 1708 } 1709 1710 static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc) 1711 { 1712 PC_PATCH *patch = (PC_PATCH *) pc->data; 1713 PCPatchConstructType patchConstructionType = PC_PATCH_STAR; 1714 char sub_mat_type[PETSC_MAX_PATH_LEN]; 1715 const char *prefix; 1716 PetscBool flg, dimflg, codimflg; 1717 MPI_Comm comm; 1718 PetscInt *ifields, nfields, k; 1719 PetscErrorCode ierr; 1720 1721 PetscFunctionBegin; 1722 ierr = PetscObjectGetComm((PetscObject) pc, &comm);CHKERRQ(ierr); 1723 ierr = PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);CHKERRQ(ierr); 1724 ierr = PetscOptionsHead(PetscOptionsObject, "Vertex-patch Additive Schwarz options");CHKERRQ(ierr); 1725 ierr = PetscOptionsBool("-pc_patch_save_operators", "Store all patch operators for lifetime of PC?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg);CHKERRQ(ierr); 1726 ierr = PetscOptionsBool("-pc_patch_partition_of_unity", "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg);CHKERRQ(ierr); 1727 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); 1728 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); 1729 if (dimflg && codimflg) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");CHKERRQ(ierr); 1730 ierr = PetscOptionsEnum("-pc_patch_construct_type", "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum) patchConstructionType, (PetscEnum *) &patchConstructionType, &flg);CHKERRQ(ierr); 1731 if (flg) {ierr = PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL);CHKERRQ(ierr);} 1732 ierr = PetscOptionsInt("-pc_patch_vanka_dim", "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg);CHKERRQ(ierr); 1733 ierr = PetscOptionsInt("-pc_patch_ignore_dim", "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg);CHKERRQ(ierr); 1734 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); 1735 if (flg) {ierr = PCPatchSetSubMatType(pc, sub_mat_type);CHKERRQ(ierr);} 1736 ierr = PetscOptionsBool("-pc_patch_symmetrise_sweep", "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg);CHKERRQ(ierr); 1737 1738 /* If the user has set the number of subspaces, use that for the buffer size, 1739 otherwise use a large number */ 1740 if (patch->nsubspaces <= 0) { 1741 nfields = 128; 1742 } else { 1743 nfields = patch->nsubspaces; 1744 } 1745 ierr = PetscMalloc1(nfields, &ifields);CHKERRQ(ierr); 1746 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_patch_exclude_subspaces",ifields,&nfields,&flg);CHKERRQ(ierr); 1747 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"); 1748 if (flg) { 1749 PetscHSetIClear(patch->subspaces_to_exclude); 1750 for (k = 0; k < nfields; k++) { 1751 PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]); 1752 } 1753 } 1754 ierr = PetscFree(ifields);CHKERRQ(ierr); 1755 1756 ierr = PetscOptionsBool("-pc_patch_patches_view", "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);CHKERRQ(ierr); 1757 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, "-pc_patch_cells_view", &patch->viewerCells, &patch->formatCells, &patch->viewCells);CHKERRQ(ierr); 1758 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, "-pc_patch_points_view", &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);CHKERRQ(ierr); 1759 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, "-pc_patch_section_view", &patch->viewerSection, &patch->formatSection, &patch->viewSection);CHKERRQ(ierr); 1760 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, "-pc_patch_mat_view", &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);CHKERRQ(ierr); 1761 ierr = PetscOptionsTail();CHKERRQ(ierr); 1762 patch->optionsSet = PETSC_TRUE; 1763 PetscFunctionReturn(0); 1764 } 1765 1766 static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc) 1767 { 1768 PC_PATCH *patch = (PC_PATCH*) pc->data; 1769 KSPConvergedReason reason; 1770 PetscInt i; 1771 PetscErrorCode ierr; 1772 1773 PetscFunctionBegin; 1774 if (!patch->save_operators) { 1775 /* Can't do this here because the sub KSPs don't have an operator attached yet. */ 1776 PetscFunctionReturn(0); 1777 } 1778 for (i = 0; i < patch->npatch; ++i) { 1779 if (!patch->ksp[i]->setfromoptionscalled) { 1780 ierr = KSPSetFromOptions(patch->ksp[i]);CHKERRQ(ierr); 1781 } 1782 ierr = KSPSetUp(patch->ksp[i]);CHKERRQ(ierr); 1783 ierr = KSPGetConvergedReason(patch->ksp[i], &reason);CHKERRQ(ierr); 1784 if (reason == KSP_DIVERGED_PCSETUP_FAILED) pc->failedreason = PC_SUBPC_ERROR; 1785 } 1786 PetscFunctionReturn(0); 1787 } 1788 1789 static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer) 1790 { 1791 PC_PATCH *patch = (PC_PATCH *) pc->data; 1792 PetscViewer sviewer; 1793 PetscBool isascii; 1794 PetscMPIInt rank; 1795 PetscErrorCode ierr; 1796 1797 PetscFunctionBegin; 1798 /* TODO Redo tabbing with set tbas in new style */ 1799 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 1800 if (!isascii) PetscFunctionReturn(0); 1801 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);CHKERRQ(ierr); 1802 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1803 ierr = PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);CHKERRQ(ierr); 1804 ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");CHKERRQ(ierr); 1805 if (patch->partition_of_unity) {ierr = PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");CHKERRQ(ierr);} 1806 else {ierr = PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");CHKERRQ(ierr);} 1807 if (patch->symmetrise_sweep) {ierr = PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");CHKERRQ(ierr);} 1808 else {ierr = PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");CHKERRQ(ierr);} 1809 if (!patch->save_operators) {ierr = PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");CHKERRQ(ierr);} 1810 else {ierr = PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");CHKERRQ(ierr);} 1811 if (patch->patchconstructop == PCPatchConstruct_Star) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");CHKERRQ(ierr);} 1812 else if (patch->patchconstructop == PCPatchConstruct_Vanka) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");CHKERRQ(ierr);} 1813 else if (patch->patchconstructop == PCPatchConstruct_User) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");CHKERRQ(ierr);} 1814 else {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");CHKERRQ(ierr);} 1815 ierr = PetscViewerASCIIPrintf(viewer, "KSP on patches (all same):\n");CHKERRQ(ierr); 1816 if (patch->ksp) { 1817 ierr = PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 1818 if (!rank) { 1819 ierr = PetscViewerASCIIPushTab(sviewer);CHKERRQ(ierr); 1820 ierr = KSPView(patch->ksp[0], sviewer);CHKERRQ(ierr); 1821 ierr = PetscViewerASCIIPopTab(sviewer);CHKERRQ(ierr); 1822 } 1823 ierr = PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 1824 } else { 1825 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 1826 ierr = PetscViewerASCIIPrintf(viewer, "KSP not yet set.\n");CHKERRQ(ierr); 1827 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1828 } 1829 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 1830 PetscFunctionReturn(0); 1831 } 1832 1833 /*MC 1834 PCPATCH - A PC object that encapsulates flexible definition of blocks for overlapping and non-overlapping 1835 small block additive preconditioners. Block definition is based on topology from 1836 a DM and equation numbering from a PetscSection. 1837 1838 Options Database Keys: 1839 + -pc_patch_cells_view - Views the process local cell numbers for each patch 1840 . -pc_patch_points_view - Views the process local mesh point numbers for each patch 1841 . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch 1842 . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary 1843 - -pc_patch_sub_mat_view - Views the matrix associated with each patch 1844 1845 Level: intermediate 1846 1847 .seealso: PCType, PCCreate(), PCSetType() 1848 M*/ 1849 PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc) 1850 { 1851 PC_PATCH *patch; 1852 PetscErrorCode ierr; 1853 1854 PetscFunctionBegin; 1855 ierr = PetscNewLog(pc, &patch);CHKERRQ(ierr); 1856 1857 if (patch->subspaces_to_exclude) { 1858 PetscHSetIDestroy(&patch->subspaces_to_exclude); 1859 } 1860 PetscHSetICreate(&patch->subspaces_to_exclude); 1861 1862 /* Set some defaults */ 1863 patch->combined = PETSC_FALSE; 1864 patch->save_operators = PETSC_TRUE; 1865 patch->partition_of_unity = PETSC_FALSE; 1866 patch->codim = -1; 1867 patch->dim = -1; 1868 patch->vankadim = -1; 1869 patch->ignoredim = -1; 1870 patch->patchconstructop = PCPatchConstruct_Star; 1871 patch->symmetrise_sweep = PETSC_FALSE; 1872 patch->npatch = 0; 1873 patch->userIS = NULL; 1874 patch->optionsSet = PETSC_FALSE; 1875 patch->iterationSet = NULL; 1876 patch->user_patches = PETSC_FALSE; 1877 ierr = PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 1878 patch->viewPatches = PETSC_FALSE; 1879 patch->viewCells = PETSC_FALSE; 1880 patch->viewPoints = PETSC_FALSE; 1881 patch->viewSection = PETSC_FALSE; 1882 patch->viewMatrix = PETSC_FALSE; 1883 1884 pc->data = (void *) patch; 1885 pc->ops->apply = PCApply_PATCH; 1886 pc->ops->applytranspose = 0; /* PCApplyTranspose_PATCH; */ 1887 pc->ops->setup = PCSetUp_PATCH; 1888 pc->ops->reset = PCReset_PATCH; 1889 pc->ops->destroy = PCDestroy_PATCH; 1890 pc->ops->setfromoptions = PCSetFromOptions_PATCH; 1891 pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH; 1892 pc->ops->view = PCView_PATCH; 1893 pc->ops->applyrichardson = 0; 1894 1895 PetscFunctionReturn(0); 1896 } 1897