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 PCPatchSetLocalComposition(PC pc, PCCompositeType type) 287 { 288 PC_PATCH *patch = (PC_PATCH *) pc->data; 289 PetscFunctionBegin; 290 if (type != PC_COMPOSITE_ADDITIVE && type != PC_COMPOSITE_MULTIPLICATIVE) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Only supports additive or multiplicative as the local type"); 291 patch->local_composition_type = type; 292 PetscFunctionReturn(0); 293 } 294 295 /* TODO: Docs */ 296 PetscErrorCode PCPatchGetLocalComposition(PC pc, PCCompositeType *type) 297 { 298 PC_PATCH *patch = (PC_PATCH *) pc->data; 299 PetscFunctionBegin; 300 *type = patch->local_composition_type; 301 PetscFunctionReturn(0); 302 } 303 304 /* TODO: Docs */ 305 PetscErrorCode PCPatchSetSubMatType(PC pc, MatType sub_mat_type) 306 { 307 PC_PATCH *patch = (PC_PATCH *) pc->data; 308 PetscErrorCode ierr; 309 310 PetscFunctionBegin; 311 if (patch->sub_mat_type) {ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr);} 312 ierr = PetscStrallocpy(sub_mat_type, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 /* TODO: Docs */ 317 PetscErrorCode PCPatchGetSubMatType(PC pc, MatType *sub_mat_type) 318 { 319 PC_PATCH *patch = (PC_PATCH *) pc->data; 320 PetscFunctionBegin; 321 *sub_mat_type = patch->sub_mat_type; 322 PetscFunctionReturn(0); 323 } 324 325 /* TODO: Docs */ 326 PetscErrorCode PCPatchSetCellNumbering(PC pc, PetscSection cellNumbering) 327 { 328 PC_PATCH *patch = (PC_PATCH *) pc->data; 329 PetscErrorCode ierr; 330 331 PetscFunctionBegin; 332 patch->cellNumbering = cellNumbering; 333 ierr = PetscObjectReference((PetscObject) cellNumbering);CHKERRQ(ierr); 334 PetscFunctionReturn(0); 335 } 336 337 /* TODO: Docs */ 338 PetscErrorCode PCPatchGetCellNumbering(PC pc, PetscSection *cellNumbering) 339 { 340 PC_PATCH *patch = (PC_PATCH *) pc->data; 341 PetscFunctionBegin; 342 *cellNumbering = patch->cellNumbering; 343 PetscFunctionReturn(0); 344 } 345 346 /* TODO: Docs */ 347 PetscErrorCode PCPatchSetConstructType(PC pc, PCPatchConstructType ctype, PetscErrorCode (*func)(PC, PetscInt *, IS **, IS *, void *), void *ctx) 348 { 349 PC_PATCH *patch = (PC_PATCH *) pc->data; 350 351 PetscFunctionBegin; 352 patch->ctype = ctype; 353 switch (ctype) { 354 case PC_PATCH_STAR: 355 patch->user_patches = PETSC_FALSE; 356 patch->patchconstructop = PCPatchConstruct_Star; 357 break; 358 case PC_PATCH_VANKA: 359 patch->user_patches = PETSC_FALSE; 360 patch->patchconstructop = PCPatchConstruct_Vanka; 361 break; 362 case PC_PATCH_USER: 363 case PC_PATCH_PYTHON: 364 patch->user_patches = PETSC_TRUE; 365 patch->patchconstructop = PCPatchConstruct_User; 366 if (func) { 367 patch->userpatchconstructionop = func; 368 patch->userpatchconstructctx = ctx; 369 } 370 break; 371 default: 372 SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype); 373 } 374 PetscFunctionReturn(0); 375 } 376 377 /* TODO: Docs */ 378 PetscErrorCode PCPatchGetConstructType(PC pc, PCPatchConstructType *ctype, PetscErrorCode (**func)(PC, PetscInt *, IS **, IS *, void *), void **ctx) 379 { 380 PC_PATCH *patch = (PC_PATCH *) pc->data; 381 382 PetscFunctionBegin; 383 *ctype = patch->ctype; 384 switch (patch->ctype) { 385 case PC_PATCH_STAR: 386 case PC_PATCH_VANKA: 387 break; 388 case PC_PATCH_USER: 389 case PC_PATCH_PYTHON: 390 *func = patch->userpatchconstructionop; 391 *ctx = patch->userpatchconstructctx; 392 break; 393 default: 394 SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_USER, "Unknown patch construction type %D", (PetscInt) patch->ctype); 395 } 396 PetscFunctionReturn(0); 397 } 398 399 /* TODO: Docs */ 400 PetscErrorCode PCPatchSetDiscretisationInfo(PC pc, PetscInt nsubspaces, DM *dms, PetscInt *bs, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, 401 const PetscInt *subspaceOffsets, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 402 { 403 PC_PATCH *patch = (PC_PATCH *) pc->data; 404 DM dm; 405 PetscSF *sfs; 406 PetscInt cStart, cEnd, i, j; 407 PetscErrorCode ierr; 408 409 PetscFunctionBegin; 410 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 411 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 412 ierr = PetscMalloc1(nsubspaces, &sfs);CHKERRQ(ierr); 413 ierr = PetscMalloc1(nsubspaces, &patch->dofSection);CHKERRQ(ierr); 414 ierr = PetscMalloc1(nsubspaces, &patch->bs);CHKERRQ(ierr); 415 ierr = PetscMalloc1(nsubspaces, &patch->nodesPerCell);CHKERRQ(ierr); 416 ierr = PetscMalloc1(nsubspaces, &patch->cellNodeMap);CHKERRQ(ierr); 417 ierr = PetscMalloc1(nsubspaces+1, &patch->subspaceOffsets);CHKERRQ(ierr); 418 419 patch->nsubspaces = nsubspaces; 420 patch->totalDofsPerCell = 0; 421 for (i = 0; i < nsubspaces; ++i) { 422 ierr = DMGetDefaultSection(dms[i], &patch->dofSection[i]);CHKERRQ(ierr); 423 ierr = PetscObjectReference((PetscObject) patch->dofSection[i]);CHKERRQ(ierr); 424 ierr = DMGetDefaultSF(dms[i], &sfs[i]);CHKERRQ(ierr); 425 patch->bs[i] = bs[i]; 426 patch->nodesPerCell[i] = nodesPerCell[i]; 427 patch->totalDofsPerCell += nodesPerCell[i]*bs[i]; 428 ierr = PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);CHKERRQ(ierr); 429 for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j]; 430 patch->subspaceOffsets[i] = subspaceOffsets[i]; 431 } 432 ierr = PCPatchCreateDefaultSF_Private(pc, nsubspaces, sfs, patch->bs);CHKERRQ(ierr); 433 ierr = PetscFree(sfs);CHKERRQ(ierr); 434 435 patch->subspaceOffsets[nsubspaces] = subspaceOffsets[nsubspaces]; 436 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);CHKERRQ(ierr); 437 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 441 /* TODO: Docs */ 442 PetscErrorCode PCPatchSetDiscretisationInfoCombined(PC pc, DM dm, PetscInt *nodesPerCell, const PetscInt **cellNodeMap, PetscInt numGhostBcs, const PetscInt *ghostBcNodes, PetscInt numGlobalBcs, const PetscInt *globalBcNodes) 443 { 444 PC_PATCH *patch = (PC_PATCH *) pc->data; 445 PetscInt cStart, cEnd, i, j; 446 PetscErrorCode ierr; 447 448 PetscFunctionBegin; 449 patch->combined = PETSC_TRUE; 450 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 451 ierr = DMGetNumFields(dm, &patch->nsubspaces);CHKERRQ(ierr); 452 ierr = PetscCalloc1(patch->nsubspaces, &patch->dofSection);CHKERRQ(ierr); 453 ierr = PetscMalloc1(patch->nsubspaces, &patch->bs);CHKERRQ(ierr); 454 ierr = PetscMalloc1(patch->nsubspaces, &patch->nodesPerCell);CHKERRQ(ierr); 455 ierr = PetscMalloc1(patch->nsubspaces, &patch->cellNodeMap);CHKERRQ(ierr); 456 ierr = PetscCalloc1(patch->nsubspaces+1, &patch->subspaceOffsets);CHKERRQ(ierr); 457 ierr = DMGetDefaultSection(dm, &patch->dofSection[0]);CHKERRQ(ierr); 458 ierr = PetscObjectReference((PetscObject) patch->dofSection[0]);CHKERRQ(ierr); 459 ierr = PetscSectionGetStorageSize(patch->dofSection[0], &patch->subspaceOffsets[patch->nsubspaces]);CHKERRQ(ierr); 460 patch->totalDofsPerCell = 0; 461 for (i = 0; i < patch->nsubspaces; ++i) { 462 patch->bs[i] = 1; 463 patch->nodesPerCell[i] = nodesPerCell[i]; 464 patch->totalDofsPerCell += nodesPerCell[i]; 465 ierr = PetscMalloc1((cEnd-cStart)*nodesPerCell[i], &patch->cellNodeMap[i]);CHKERRQ(ierr); 466 for (j = 0; j < (cEnd-cStart)*nodesPerCell[i]; ++j) patch->cellNodeMap[i][j] = cellNodeMap[i][j]; 467 } 468 ierr = DMGetDefaultSF(dm, &patch->defaultSF);CHKERRQ(ierr); 469 ierr = PetscObjectReference((PetscObject) patch->defaultSF);CHKERRQ(ierr); 470 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGhostBcs, ghostBcNodes, PETSC_COPY_VALUES, &patch->ghostBcNodes);CHKERRQ(ierr); 471 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalBcs, globalBcNodes, PETSC_COPY_VALUES, &patch->globalBcNodes);CHKERRQ(ierr); 472 PetscFunctionReturn(0); 473 } 474 475 /*@C 476 477 PCPatchSetComputeFunction - Set the callback used to compute patch residuals 478 479 Input Parameters: 480 + pc - The PC 481 . func - The callback 482 - ctx - The user context 483 484 Level: advanced 485 486 Note: 487 The callback has signature: 488 + usercomputef(pc, point, x, f, cellIS, n, u, ctx) 489 + pc - The PC 490 + point - The point 491 + x - The input solution (not used in linear problems) 492 + f - The patch residual vector 493 + cellIS - An array of the cell numbers 494 + n - The size of g2l 495 + g2l - The global to local dof translation table 496 + ctx - The user context 497 and can assume that the matrix entries have been set to zero before the call. 498 499 .seealso: PCPatchSetComputeOperator(), PCPatchGetComputeOperator(), PCPatchSetDiscretisationInfo() 500 @*/ 501 PetscErrorCode PCPatchSetComputeFunction(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Vec, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx) 502 { 503 PC_PATCH *patch = (PC_PATCH *) pc->data; 504 505 PetscFunctionBegin; 506 patch->usercomputef = func; 507 patch->usercomputefctx = ctx; 508 PetscFunctionReturn(0); 509 } 510 511 /*@C 512 513 PCPatchSetComputeOperator - Set the callback used to compute patch matrices 514 515 Input Parameters: 516 + pc - The PC 517 . func - The callback 518 - ctx - The user context 519 520 Level: advanced 521 522 Note: 523 The callback has signature: 524 + usercomputeop(pc, point, x, mat, cellIS, n, u, ctx) 525 + pc - The PC 526 + point - The point 527 + x - The input solution (not used in linear problems) 528 + mat - The patch matrix 529 + cellIS - An array of the cell numbers 530 + n - The size of g2l 531 + g2l - The global to local dof translation table 532 + ctx - The user context 533 and can assume that the matrix entries have been set to zero before the call. 534 535 .seealso: PCPatchGetComputeOperator(), PCPatchSetComputeFunction(), PCPatchSetDiscretisationInfo() 536 @*/ 537 PetscErrorCode PCPatchSetComputeOperator(PC pc, PetscErrorCode (*func)(PC, PetscInt, Vec, Mat, IS, PetscInt, const PetscInt *, const PetscInt *, void *), void *ctx) 538 { 539 PC_PATCH *patch = (PC_PATCH *) pc->data; 540 541 PetscFunctionBegin; 542 patch->usercomputeop = func; 543 patch->usercomputeopctx = ctx; 544 PetscFunctionReturn(0); 545 } 546 547 /* On entry, ht contains the topological entities whose dofs we are responsible for solving for; 548 on exit, cht contains all the topological entities we need to compute their residuals. 549 In full generality this should incorporate knowledge of the sparsity pattern of the matrix; 550 here we assume a standard FE sparsity pattern.*/ 551 /* TODO: Use DMPlexGetAdjacency() */ 552 static PetscErrorCode PCPatchCompleteCellPatch(PC pc, PetscHSetI ht, PetscHSetI cht) 553 { 554 DM dm; 555 PetscHashIter hi; 556 PetscInt point; 557 PetscInt *star = NULL, *closure = NULL; 558 PetscInt ignoredim, iStart = 0, iEnd = -1, starSize, closureSize, si, ci; 559 PetscErrorCode ierr; 560 561 PetscFunctionBegin; 562 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 563 ierr = PCPatchGetIgnoreDim(pc, &ignoredim);CHKERRQ(ierr); 564 if (ignoredim >= 0) {ierr = DMPlexGetDepthStratum(dm, ignoredim, &iStart, &iEnd);CHKERRQ(ierr);} 565 ierr = PetscHSetIClear(cht);CHKERRQ(ierr); 566 PetscHashIterBegin(ht, hi); 567 while (!PetscHashIterAtEnd(ht, hi)) { 568 569 PetscHashIterGetKey(ht, hi, point); 570 PetscHashIterNext(ht, hi); 571 572 /* Loop over all the cells that this point connects to */ 573 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_FALSE, &starSize, &star);CHKERRQ(ierr); 574 for (si = 0; si < starSize*2; si += 2) { 575 const PetscInt ownedpoint = star[si]; 576 /* TODO Check for point in cht before running through closure again */ 577 /* now loop over all entities in the closure of that cell */ 578 ierr = DMPlexGetTransitiveClosure(dm, ownedpoint, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 579 for (ci = 0; ci < closureSize*2; ci += 2) { 580 const PetscInt seenpoint = closure[ci]; 581 if (ignoredim >= 0 && seenpoint >= iStart && seenpoint < iEnd) continue; 582 ierr = PetscHSetIAdd(cht, seenpoint);CHKERRQ(ierr); 583 } 584 } 585 } 586 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, NULL, &closure);CHKERRQ(ierr); 587 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_FALSE, NULL, &star);CHKERRQ(ierr); 588 PetscFunctionReturn(0); 589 } 590 591 static PetscErrorCode PCPatchGetGlobalDofs(PC pc, PetscSection dofSection[], PetscInt f, PetscBool combined, PetscInt p, PetscInt *dof, PetscInt *off) 592 { 593 PetscErrorCode ierr; 594 595 PetscFunctionBegin; 596 if (combined) { 597 if (f < 0) { 598 if (dof) {ierr = PetscSectionGetDof(dofSection[0], p, dof);CHKERRQ(ierr);} 599 if (off) {ierr = PetscSectionGetOffset(dofSection[0], p, off);CHKERRQ(ierr);} 600 } else { 601 if (dof) {ierr = PetscSectionGetFieldDof(dofSection[0], p, f, dof);CHKERRQ(ierr);} 602 if (off) {ierr = PetscSectionGetFieldOffset(dofSection[0], p, f, off);CHKERRQ(ierr);} 603 } 604 } else { 605 if (f < 0) { 606 PC_PATCH *patch = (PC_PATCH *) pc->data; 607 PetscInt fdof, g; 608 609 if (dof) { 610 *dof = 0; 611 for (g = 0; g < patch->nsubspaces; ++g) { 612 ierr = PetscSectionGetDof(dofSection[g], p, &fdof);CHKERRQ(ierr); 613 *dof += fdof; 614 } 615 } 616 if (off) { 617 *off = 0; 618 for (g = 0; g < patch->nsubspaces; ++g) { 619 ierr = PetscSectionGetOffset(dofSection[g], p, &fdof);CHKERRQ(ierr); 620 *off += fdof; 621 } 622 } 623 } else { 624 if (dof) {ierr = PetscSectionGetDof(dofSection[f], p, dof);CHKERRQ(ierr);} 625 if (off) {ierr = PetscSectionGetOffset(dofSection[f], p, off);CHKERRQ(ierr);} 626 } 627 } 628 PetscFunctionReturn(0); 629 } 630 631 /* Given a hash table with a set of topological entities (pts), compute the degrees of 632 freedom in global concatenated numbering on those entities. 633 For Vanka smoothing, this needs to do something special: ignore dofs of the 634 constraint subspace on entities that aren't the base entity we're building the patch 635 around. */ 636 static PetscErrorCode PCPatchGetPointDofs(PC pc, PetscHSetI pts, PetscHSetI dofs, PetscInt base, PetscHSetI* subspaces_to_exclude) 637 { 638 PC_PATCH *patch = (PC_PATCH *) pc->data; 639 PetscHashIter hi; 640 PetscInt ldof, loff; 641 PetscInt k, p; 642 PetscErrorCode ierr; 643 644 PetscFunctionBegin; 645 ierr = PetscHSetIClear(dofs);CHKERRQ(ierr); 646 for (k = 0; k < patch->nsubspaces; ++k) { 647 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 648 PetscInt bs = patch->bs[k]; 649 PetscInt j, l; 650 651 if (subspaces_to_exclude != NULL) { 652 PetscBool should_exclude_k = PETSC_FALSE; 653 PetscHSetIHas(*subspaces_to_exclude, k, &should_exclude_k); 654 if (should_exclude_k) { 655 /* only get this subspace dofs at the base entity, not any others */ 656 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, base, &ldof, &loff);CHKERRQ(ierr); 657 if (0 == ldof) continue; 658 for (j = loff; j < ldof + loff; ++j) { 659 for (l = 0; l < bs; ++l) { 660 PetscInt dof = bs*j + l + subspaceOffset; 661 ierr = PetscHSetIAdd(dofs, dof);CHKERRQ(ierr); 662 } 663 } 664 continue; /* skip the other dofs of this subspace */ 665 } 666 } 667 668 PetscHashIterBegin(pts, hi); 669 while (!PetscHashIterAtEnd(pts, hi)) { 670 PetscHashIterGetKey(pts, hi, p); 671 PetscHashIterNext(pts, hi); 672 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, p, &ldof, &loff);CHKERRQ(ierr); 673 if (0 == ldof) continue; 674 for (j = loff; j < ldof + loff; ++j) { 675 for (l = 0; l < bs; ++l) { 676 PetscInt dof = bs*j + l + subspaceOffset; 677 ierr = PetscHSetIAdd(dofs, dof);CHKERRQ(ierr); 678 } 679 } 680 } 681 } 682 PetscFunctionReturn(0); 683 } 684 685 /* Given two hash tables A and B, compute the keys in B that are not in A, and put them in C */ 686 static PetscErrorCode PCPatchComputeSetDifference_Private(PetscHSetI A, PetscHSetI B, PetscHSetI C) 687 { 688 PetscHashIter hi; 689 PetscInt key; 690 PetscBool flg; 691 PetscErrorCode ierr; 692 693 PetscFunctionBegin; 694 ierr = PetscHSetIClear(C);CHKERRQ(ierr); 695 PetscHashIterBegin(B, hi); 696 while (!PetscHashIterAtEnd(B, hi)) { 697 PetscHashIterGetKey(B, hi, key); 698 PetscHashIterNext(B, hi); 699 ierr = PetscHSetIHas(A, key, &flg);CHKERRQ(ierr); 700 if (!flg) {ierr = PetscHSetIAdd(C, key);CHKERRQ(ierr);} 701 } 702 PetscFunctionReturn(0); 703 } 704 705 /* 706 * PCPatchCreateCellPatches - create patches. 707 * 708 * Input Parameters: 709 * + dm - The DMPlex object defining the mesh 710 * 711 * Output Parameters: 712 * + cellCounts - Section with counts of cells around each vertex 713 * . cells - IS of the cell point indices of cells in each patch 714 * . pointCounts - Section with counts of cells around each vertex 715 * - point - IS of the cell point indices of cells in each patch 716 */ 717 static PetscErrorCode PCPatchCreateCellPatches(PC pc) 718 { 719 PC_PATCH *patch = (PC_PATCH *) pc->data; 720 DMLabel ghost = NULL; 721 DM dm, plex; 722 PetscHSetI ht, cht; 723 PetscSection cellCounts, pointCounts; 724 PetscInt *cellsArray, *pointsArray; 725 PetscInt numCells, numPoints; 726 const PetscInt *leaves; 727 PetscInt nleaves, pStart, pEnd, cStart, cEnd, vStart, vEnd, v; 728 PetscBool isFiredrake; 729 PetscErrorCode ierr; 730 731 PetscFunctionBegin; 732 /* Used to keep track of the cells in the patch. */ 733 ierr = PetscHSetICreate(&ht);CHKERRQ(ierr); 734 ierr = PetscHSetICreate(&cht);CHKERRQ(ierr); 735 736 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 737 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONGSTATE, "DM not yet set on patch PC\n"); 738 ierr = DMConvert(dm, DMPLEX, &plex);CHKERRQ(ierr); 739 ierr = DMPlexGetChart(plex, &pStart, &pEnd);CHKERRQ(ierr); 740 ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr); 741 742 if (patch->user_patches) { 743 ierr = patch->userpatchconstructionop(pc, &patch->npatch, &patch->userIS, &patch->iterationSet, patch->userpatchconstructctx);CHKERRQ(ierr); 744 vStart = 0; vEnd = patch->npatch; 745 } else if (patch->codim < 0) { 746 if (patch->dim < 0) {ierr = DMPlexGetDepthStratum(plex, 0, &vStart, &vEnd);CHKERRQ(ierr);} 747 else {ierr = DMPlexGetDepthStratum(plex, patch->dim, &vStart, &vEnd);CHKERRQ(ierr);} 748 } else {ierr = DMPlexGetHeightStratum(plex, patch->codim, &vStart, &vEnd);CHKERRQ(ierr);} 749 patch->npatch = vEnd - vStart; 750 751 /* These labels mark the owned points. We only create patches around points that this process owns. */ 752 ierr = DMHasLabel(dm, "pyop2_ghost", &isFiredrake);CHKERRQ(ierr); 753 if (isFiredrake) { 754 ierr = DMGetLabel(dm, "pyop2_ghost", &ghost);CHKERRQ(ierr); 755 ierr = DMLabelCreateIndex(ghost, pStart, pEnd);CHKERRQ(ierr); 756 } else { 757 PetscSF sf; 758 759 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 760 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 761 nleaves = PetscMax(nleaves, 0); 762 } 763 764 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->cellCounts);CHKERRQ(ierr); 765 ierr = PetscObjectSetName((PetscObject) patch->cellCounts, "Patch Cell Layout");CHKERRQ(ierr); 766 cellCounts = patch->cellCounts; 767 ierr = PetscSectionSetChart(cellCounts, vStart, vEnd);CHKERRQ(ierr); 768 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->pointCounts);CHKERRQ(ierr); 769 ierr = PetscObjectSetName((PetscObject) patch->pointCounts, "Patch Point Layout");CHKERRQ(ierr); 770 pointCounts = patch->pointCounts; 771 ierr = PetscSectionSetChart(pointCounts, vStart, vEnd);CHKERRQ(ierr); 772 /* Count cells and points in the patch surrounding each entity */ 773 for (v = vStart; v < vEnd; ++v) { 774 PetscHashIter hi; 775 PetscInt chtSize, loc = -1; 776 PetscBool flg; 777 778 if (!patch->user_patches) { 779 if (ghost) {ierr = DMLabelHasPoint(ghost, v, &flg);CHKERRQ(ierr);} 780 else {ierr = PetscFindInt(v, nleaves, leaves, &loc);CHKERRQ(ierr); flg = loc >=0 ? PETSC_TRUE : PETSC_FALSE;} 781 /* Not an owned entity, don't make a cell patch. */ 782 if (flg) continue; 783 } 784 785 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 786 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 787 ierr = PetscHSetIGetSize(cht, &chtSize);CHKERRQ(ierr); 788 /* empty patch, continue */ 789 if (chtSize == 0) continue; 790 791 /* safe because size(cht) > 0 from above */ 792 PetscHashIterBegin(cht, hi); 793 while (!PetscHashIterAtEnd(cht, hi)) { 794 PetscInt point, pdof; 795 796 PetscHashIterGetKey(cht, hi, point); 797 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 798 if (pdof) {ierr = PetscSectionAddDof(pointCounts, v, 1);CHKERRQ(ierr);} 799 if (point >= cStart && point < cEnd) {ierr = PetscSectionAddDof(cellCounts, v, 1);CHKERRQ(ierr);} 800 PetscHashIterNext(cht, hi); 801 } 802 } 803 if (isFiredrake) {ierr = DMLabelDestroyIndex(ghost);CHKERRQ(ierr);} 804 805 ierr = PetscSectionSetUp(cellCounts);CHKERRQ(ierr); 806 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 807 ierr = PetscMalloc1(numCells, &cellsArray);CHKERRQ(ierr); 808 ierr = PetscSectionSetUp(pointCounts);CHKERRQ(ierr); 809 ierr = PetscSectionGetStorageSize(pointCounts, &numPoints);CHKERRQ(ierr); 810 ierr = PetscMalloc1(numPoints, &pointsArray);CHKERRQ(ierr); 811 812 /* Now that we know how much space we need, run through again and actually remember the cells. */ 813 for (v = vStart; v < vEnd; v++ ) { 814 PetscHashIter hi; 815 PetscInt dof, off, cdof, coff, pdof, n = 0, cn = 0; 816 817 ierr = PetscSectionGetDof(pointCounts, v, &dof);CHKERRQ(ierr); 818 ierr = PetscSectionGetOffset(pointCounts, v, &off);CHKERRQ(ierr); 819 ierr = PetscSectionGetDof(cellCounts, v, &cdof);CHKERRQ(ierr); 820 ierr = PetscSectionGetOffset(cellCounts, v, &coff);CHKERRQ(ierr); 821 if (dof <= 0) continue; 822 ierr = patch->patchconstructop((void *) patch, dm, v, ht);CHKERRQ(ierr); 823 ierr = PCPatchCompleteCellPatch(pc, ht, cht);CHKERRQ(ierr); 824 PetscHashIterBegin(cht, hi); 825 while (!PetscHashIterAtEnd(cht, hi)) { 826 PetscInt point; 827 828 PetscHashIterGetKey(cht, hi, point); 829 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, -1, patch->combined, point, &pdof, NULL);CHKERRQ(ierr); 830 if (pdof) {pointsArray[off + n++] = point;} 831 if (point >= cStart && point < cEnd) {cellsArray[coff + cn++] = point;} 832 PetscHashIterNext(cht, hi); 833 } 834 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); 835 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); 836 } 837 ierr = PetscHSetIDestroy(&ht);CHKERRQ(ierr); 838 ierr = PetscHSetIDestroy(&cht);CHKERRQ(ierr); 839 ierr = DMDestroy(&plex);CHKERRQ(ierr); 840 841 ierr = ISCreateGeneral(PETSC_COMM_SELF, numCells, cellsArray, PETSC_OWN_POINTER, &patch->cells);CHKERRQ(ierr); 842 ierr = PetscObjectSetName((PetscObject) patch->cells, "Patch Cells");CHKERRQ(ierr); 843 if (patch->viewCells) { 844 ierr = ObjectView((PetscObject) patch->cellCounts, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 845 ierr = ObjectView((PetscObject) patch->cells, patch->viewerCells, patch->formatCells);CHKERRQ(ierr); 846 } 847 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints, pointsArray, PETSC_OWN_POINTER, &patch->points);CHKERRQ(ierr); 848 ierr = PetscObjectSetName((PetscObject) patch->points, "Patch Points");CHKERRQ(ierr); 849 if (patch->viewPoints) { 850 ierr = ObjectView((PetscObject) patch->pointCounts, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 851 ierr = ObjectView((PetscObject) patch->points, patch->viewerPoints, patch->formatPoints);CHKERRQ(ierr); 852 } 853 PetscFunctionReturn(0); 854 } 855 856 /* 857 * PCPatchCreateCellPatchDiscretisationInfo - Build the dof maps for cell patches 858 * 859 * Input Parameters: 860 * + dm - The DMPlex object defining the mesh 861 * . cellCounts - Section with counts of cells around each vertex 862 * . cells - IS of the cell point indices of cells in each patch 863 * . cellNumbering - Section mapping plex cell points to Firedrake cell indices. 864 * . nodesPerCell - number of nodes per cell. 865 * - cellNodeMap - map from cells to node indices (nodesPerCell * numCells) 866 * 867 * Output Parameters: 868 * + dofs - IS of local dof numbers of each cell in the patch, where local is a patch local numbering 869 * . gtolCounts - Section with counts of dofs per cell patch 870 * - gtol - IS mapping from global dofs to local dofs for each patch. 871 */ 872 static PetscErrorCode PCPatchCreateCellPatchDiscretisationInfo(PC pc) 873 { 874 PC_PATCH *patch = (PC_PATCH *) pc->data; 875 PetscSection cellCounts = patch->cellCounts; 876 PetscSection pointCounts = patch->pointCounts; 877 PetscSection gtolCounts, gtolCountsWithArtificial = NULL, gtolCountsWithAll = NULL; 878 IS cells = patch->cells; 879 IS points = patch->points; 880 PetscSection cellNumbering = patch->cellNumbering; 881 PetscInt Nf = patch->nsubspaces; 882 PetscInt numCells, numPoints; 883 PetscInt numDofs; 884 PetscInt numGlobalDofs, numGlobalDofsWithArtificial, numGlobalDofsWithAll; 885 PetscInt totalDofsPerCell = patch->totalDofsPerCell; 886 PetscInt vStart, vEnd, v; 887 const PetscInt *cellsArray, *pointsArray; 888 PetscInt *newCellsArray = NULL; 889 PetscInt *dofsArray = NULL; 890 PetscInt *dofsArrayWithArtificial = NULL; 891 PetscInt *dofsArrayWithAll = NULL; 892 PetscInt *offsArray = NULL; 893 PetscInt *offsArrayWithArtificial = NULL; 894 PetscInt *offsArrayWithAll = NULL; 895 PetscInt *asmArray = NULL; 896 PetscInt *asmArrayWithArtificial = NULL; 897 PetscInt *asmArrayWithAll = NULL; 898 PetscInt *globalDofsArray = NULL; 899 PetscInt *globalDofsArrayWithArtificial = NULL; 900 PetscInt *globalDofsArrayWithAll = NULL; 901 PetscInt globalIndex = 0; 902 PetscInt key = 0; 903 PetscInt asmKey = 0; 904 DM dm = NULL; 905 const PetscInt *bcNodes = NULL; 906 PetscHMapI ht; 907 PetscHMapI htWithArtificial; 908 PetscHMapI htWithAll; 909 PetscHSetI globalBcs; 910 PetscInt numBcs; 911 PetscHSetI ownedpts, seenpts, owneddofs, seendofs, artificialbcs; 912 PetscInt pStart, pEnd, p, i; 913 char option[PETSC_MAX_PATH_LEN]; 914 PetscBool isNonlinear; 915 PetscErrorCode ierr; 916 917 PetscFunctionBegin; 918 919 ierr = PCGetDM(pc, &dm); CHKERRQ(ierr); 920 /* dofcounts section is cellcounts section * dofPerCell */ 921 ierr = PetscSectionGetStorageSize(cellCounts, &numCells);CHKERRQ(ierr); 922 ierr = PetscSectionGetStorageSize(patch->pointCounts, &numPoints);CHKERRQ(ierr); 923 numDofs = numCells * totalDofsPerCell; 924 ierr = PetscMalloc1(numDofs, &dofsArray);CHKERRQ(ierr); 925 ierr = PetscMalloc1(numPoints*Nf, &offsArray);CHKERRQ(ierr); 926 ierr = PetscMalloc1(numDofs, &asmArray);CHKERRQ(ierr); 927 ierr = PetscMalloc1(numCells, &newCellsArray);CHKERRQ(ierr); 928 ierr = PetscSectionGetChart(cellCounts, &vStart, &vEnd);CHKERRQ(ierr); 929 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCounts);CHKERRQ(ierr); 930 gtolCounts = patch->gtolCounts; 931 ierr = PetscSectionSetChart(gtolCounts, vStart, vEnd);CHKERRQ(ierr); 932 ierr = PetscObjectSetName((PetscObject) patch->gtolCounts, "Patch Global Index Section");CHKERRQ(ierr); 933 934 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) 935 { 936 ierr = PetscMalloc1(numPoints*Nf, &offsArrayWithArtificial);CHKERRQ(ierr); 937 ierr = PetscMalloc1(numDofs, &asmArrayWithArtificial);CHKERRQ(ierr); 938 ierr = PetscMalloc1(numDofs, &dofsArrayWithArtificial);CHKERRQ(ierr); 939 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithArtificial);CHKERRQ(ierr); 940 gtolCountsWithArtificial = patch->gtolCountsWithArtificial; 941 ierr = PetscSectionSetChart(gtolCountsWithArtificial, vStart, vEnd);CHKERRQ(ierr); 942 ierr = PetscObjectSetName((PetscObject) patch->gtolCountsWithArtificial, "Patch Global Index Section Including Artificial BCs");CHKERRQ(ierr); 943 } 944 945 isNonlinear = patch->isNonlinear; 946 if(isNonlinear) 947 { 948 ierr = PetscMalloc1(numPoints*Nf, &offsArrayWithAll);CHKERRQ(ierr); 949 ierr = PetscMalloc1(numDofs, &asmArrayWithAll);CHKERRQ(ierr); 950 ierr = PetscMalloc1(numDofs, &dofsArrayWithAll);CHKERRQ(ierr); 951 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->gtolCountsWithAll);CHKERRQ(ierr); 952 gtolCountsWithAll = patch->gtolCountsWithAll; 953 ierr = PetscSectionSetChart(gtolCountsWithAll, vStart, vEnd);CHKERRQ(ierr); 954 ierr = PetscObjectSetName((PetscObject) patch->gtolCountsWithAll, "Patch Global Index Section Including All BCs");CHKERRQ(ierr); 955 } 956 957 /* Outside the patch loop, get the dofs that are globally-enforced Dirichlet 958 conditions */ 959 ierr = PetscHSetICreate(&globalBcs);CHKERRQ(ierr); 960 ierr = ISGetIndices(patch->ghostBcNodes, &bcNodes); CHKERRQ(ierr); 961 ierr = ISGetSize(patch->ghostBcNodes, &numBcs); CHKERRQ(ierr); 962 for (i = 0; i < numBcs; ++i) { 963 ierr = PetscHSetIAdd(globalBcs, bcNodes[i]);CHKERRQ(ierr); /* these are already in concatenated numbering */ 964 } 965 ierr = ISRestoreIndices(patch->ghostBcNodes, &bcNodes); CHKERRQ(ierr); 966 ierr = ISDestroy(&patch->ghostBcNodes); CHKERRQ(ierr); /* memory optimisation */ 967 968 /* Hash tables for artificial BC construction */ 969 ierr = PetscHSetICreate(&ownedpts);CHKERRQ(ierr); 970 ierr = PetscHSetICreate(&seenpts);CHKERRQ(ierr); 971 ierr = PetscHSetICreate(&owneddofs);CHKERRQ(ierr); 972 ierr = PetscHSetICreate(&seendofs);CHKERRQ(ierr); 973 ierr = PetscHSetICreate(&artificialbcs);CHKERRQ(ierr); 974 975 ierr = ISGetIndices(cells, &cellsArray);CHKERRQ(ierr); 976 ierr = ISGetIndices(points, &pointsArray);CHKERRQ(ierr); 977 ierr = PetscHMapICreate(&ht);CHKERRQ(ierr); 978 ierr = PetscHMapICreate(&htWithArtificial);CHKERRQ(ierr); 979 ierr = PetscHMapICreate(&htWithAll);CHKERRQ(ierr); 980 for (v = vStart; v < vEnd; ++v) { 981 PetscInt localIndex = 0; 982 PetscInt localIndexWithArtificial = 0; 983 PetscInt localIndexWithAll = 0; 984 PetscInt dof, off, i, j, k, l; 985 986 ierr = PetscHMapIClear(ht);CHKERRQ(ierr); 987 ierr = PetscHMapIClear(htWithArtificial);CHKERRQ(ierr); 988 ierr = PetscHMapIClear(htWithAll);CHKERRQ(ierr); 989 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 990 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 991 if (dof <= 0) continue; 992 993 /* Calculate the global numbers of the artificial BC dofs here first */ 994 ierr = patch->patchconstructop((void*)patch, dm, v, ownedpts); CHKERRQ(ierr); 995 ierr = PCPatchCompleteCellPatch(pc, ownedpts, seenpts); CHKERRQ(ierr); 996 ierr = PCPatchGetPointDofs(pc, ownedpts, owneddofs, v, &patch->subspaces_to_exclude); CHKERRQ(ierr); 997 ierr = PCPatchGetPointDofs(pc, seenpts, seendofs, v, NULL); CHKERRQ(ierr); 998 ierr = PCPatchComputeSetDifference_Private(owneddofs, seendofs, artificialbcs); CHKERRQ(ierr); 999 if (patch->viewPatches) { 1000 PetscHSetI globalbcdofs; 1001 PetscHashIter hi; 1002 MPI_Comm comm = PetscObjectComm((PetscObject)pc); 1003 1004 ierr = PetscHSetICreate(&globalbcdofs);CHKERRQ(ierr); 1005 ierr = PetscSynchronizedPrintf(comm, "Patch %d: owned dofs:\n", v); CHKERRQ(ierr); 1006 PetscHashIterBegin(owneddofs, hi); 1007 while (!PetscHashIterAtEnd(owneddofs, hi)) { 1008 PetscInt globalDof; 1009 1010 PetscHashIterGetKey(owneddofs, hi, globalDof); 1011 PetscHashIterNext(owneddofs, hi); 1012 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 1013 } 1014 ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr); 1015 ierr = PetscSynchronizedPrintf(comm, "Patch %d: seen dofs:\n", v); CHKERRQ(ierr); 1016 PetscHashIterBegin(seendofs, hi); 1017 while (!PetscHashIterAtEnd(seendofs, hi)) { 1018 PetscInt globalDof; 1019 PetscBool flg; 1020 1021 PetscHashIterGetKey(seendofs, hi, globalDof); 1022 PetscHashIterNext(seendofs, hi); 1023 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 1024 1025 ierr = PetscHSetIHas(globalBcs, globalDof, &flg);CHKERRQ(ierr); 1026 if (flg) {ierr = PetscHSetIAdd(globalbcdofs, globalDof);CHKERRQ(ierr);} 1027 } 1028 ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr); 1029 ierr = PetscSynchronizedPrintf(comm, "Patch %d: global BCs:\n", v); CHKERRQ(ierr); 1030 ierr = PetscHSetIGetSize(globalbcdofs, &numBcs);CHKERRQ(ierr); 1031 if (numBcs > 0) { 1032 PetscHashIterBegin(globalbcdofs, hi); 1033 while (!PetscHashIterAtEnd(globalbcdofs, hi)) { 1034 PetscInt globalDof; 1035 PetscHashIterGetKey(globalbcdofs, hi, globalDof); 1036 PetscHashIterNext(globalbcdofs, hi); 1037 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof);CHKERRQ(ierr); 1038 } 1039 } 1040 ierr = PetscSynchronizedPrintf(comm, "\n");CHKERRQ(ierr); 1041 ierr = PetscSynchronizedPrintf(comm, "Patch %d: artificial BCs:\n", v);CHKERRQ(ierr); 1042 ierr = PetscHSetIGetSize(artificialbcs, &numBcs);CHKERRQ(ierr); 1043 if (numBcs > 0) { 1044 PetscHashIterBegin(artificialbcs, hi); 1045 while (!PetscHashIterAtEnd(artificialbcs, hi)) { 1046 PetscInt globalDof; 1047 PetscHashIterGetKey(artificialbcs, hi, globalDof); 1048 PetscHashIterNext(artificialbcs, hi); 1049 ierr = PetscSynchronizedPrintf(comm, "%d ", globalDof); CHKERRQ(ierr); 1050 } 1051 } 1052 ierr = PetscSynchronizedPrintf(comm, "\n\n"); CHKERRQ(ierr); 1053 ierr = PetscHSetIDestroy(&globalbcdofs);CHKERRQ(ierr); 1054 } 1055 for (k = 0; k < patch->nsubspaces; ++k) { 1056 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1057 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1058 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1059 PetscInt bs = patch->bs[k]; 1060 1061 for (i = off; i < off + dof; ++i) { 1062 /* Walk over the cells in this patch. */ 1063 const PetscInt c = cellsArray[i]; 1064 PetscInt cell = c; 1065 1066 /* TODO Change this to an IS */ 1067 if (cellNumbering) { 1068 ierr = PetscSectionGetDof(cellNumbering, c, &cell);CHKERRQ(ierr); 1069 if (cell <= 0) SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_OUTOFRANGE, "Cell %D doesn't appear in cell numbering map", c); 1070 ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr); 1071 } 1072 newCellsArray[i] = cell; 1073 for (j = 0; j < nodesPerCell; ++j) { 1074 /* For each global dof, map it into contiguous local storage. */ 1075 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + subspaceOffset; 1076 /* finally, loop over block size */ 1077 for (l = 0; l < bs; ++l) { 1078 PetscInt localDof; 1079 PetscBool isGlobalBcDof, isArtificialBcDof; 1080 1081 /* first, check if this is either a globally enforced or locally enforced BC dof */ 1082 ierr = PetscHSetIHas(globalBcs, globalDof + l, &isGlobalBcDof);CHKERRQ(ierr); 1083 ierr = PetscHSetIHas(artificialbcs, globalDof + l, &isArtificialBcDof);CHKERRQ(ierr); 1084 1085 /* if it's either, don't ever give it a local dof number */ 1086 if (isGlobalBcDof || isArtificialBcDof) { 1087 dofsArray[globalIndex] = -1; /* don't use this in assembly in this patch */ 1088 } else { 1089 ierr = PetscHMapIGet(ht, globalDof + l, &localDof);CHKERRQ(ierr); 1090 if (localDof == -1) { 1091 localDof = localIndex++; 1092 ierr = PetscHMapISet(ht, globalDof + l, localDof);CHKERRQ(ierr); 1093 } 1094 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 1095 /* And store. */ 1096 dofsArray[globalIndex] = localDof; 1097 } 1098 1099 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1100 if (isGlobalBcDof) { 1101 dofsArrayWithArtificial[globalIndex] = -1; /* don't use this in assembly in this patch */ 1102 } else { 1103 ierr = PetscHMapIGet(htWithArtificial, globalDof + l, &localDof);CHKERRQ(ierr); 1104 if (localDof == -1) { 1105 localDof = localIndexWithArtificial++; 1106 ierr = PetscHMapISet(htWithArtificial, globalDof + l, localDof);CHKERRQ(ierr); 1107 } 1108 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 1109 /* And store.*/ 1110 dofsArrayWithArtificial[globalIndex] = localDof; 1111 } 1112 } 1113 1114 if(isNonlinear) { 1115 /* Build the dofmap for the function space with _all_ dofs, 1116 including those in any kind of boundary condition */ 1117 ierr = PetscHMapIGet(htWithAll, globalDof + l, &localDof);CHKERRQ(ierr); 1118 if (localDof == -1) { 1119 localDof = localIndexWithAll++; 1120 ierr = PetscHMapISet(htWithAll, globalDof + l, localDof);CHKERRQ(ierr); 1121 } 1122 if ( globalIndex >= numDofs ) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Found more dofs %D than expected %D", globalIndex+1, numDofs); 1123 /* And store.*/ 1124 dofsArrayWithAll[globalIndex] = localDof; 1125 } 1126 globalIndex++; 1127 } 1128 } 1129 } 1130 } 1131 /*How many local dofs in this patch? */ 1132 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1133 ierr = PetscHMapIGetSize(htWithArtificial, &dof);CHKERRQ(ierr); 1134 ierr = PetscSectionSetDof(gtolCountsWithArtificial, v, dof);CHKERRQ(ierr); 1135 } 1136 if (isNonlinear) { 1137 ierr = PetscHMapIGetSize(htWithAll, &dof);CHKERRQ(ierr); 1138 ierr = PetscSectionSetDof(gtolCountsWithAll, v, dof);CHKERRQ(ierr); 1139 } 1140 ierr = PetscHMapIGetSize(ht, &dof);CHKERRQ(ierr); 1141 ierr = PetscSectionSetDof(gtolCounts, v, dof);CHKERRQ(ierr); 1142 } 1143 if (globalIndex != numDofs) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Expected number of dofs (%d) doesn't match found number (%d)", numDofs, globalIndex); 1144 ierr = PetscSectionSetUp(gtolCounts);CHKERRQ(ierr); 1145 ierr = PetscSectionGetStorageSize(gtolCounts, &numGlobalDofs);CHKERRQ(ierr); 1146 ierr = PetscMalloc1(numGlobalDofs, &globalDofsArray);CHKERRQ(ierr); 1147 1148 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1149 ierr = PetscSectionSetUp(gtolCountsWithArtificial);CHKERRQ(ierr); 1150 ierr = PetscSectionGetStorageSize(gtolCountsWithArtificial, &numGlobalDofsWithArtificial);CHKERRQ(ierr); 1151 ierr = PetscMalloc1(numGlobalDofsWithArtificial, &globalDofsArrayWithArtificial);CHKERRQ(ierr); 1152 } 1153 if (isNonlinear) { 1154 ierr = PetscSectionSetUp(gtolCountsWithAll);CHKERRQ(ierr); 1155 ierr = PetscSectionGetStorageSize(gtolCountsWithAll, &numGlobalDofsWithAll);CHKERRQ(ierr); 1156 ierr = PetscMalloc1(numGlobalDofsWithAll, &globalDofsArrayWithAll);CHKERRQ(ierr); 1157 } 1158 /* Now populate the global to local map. This could be merged into the above loop if we were willing to deal with reallocs. */ 1159 for (v = vStart; v < vEnd; ++v) { 1160 PetscHashIter hi; 1161 PetscInt dof, off, Np, ooff, i, j, k, l; 1162 1163 ierr = PetscHMapIClear(ht);CHKERRQ(ierr); 1164 ierr = PetscHMapIClear(htWithArtificial);CHKERRQ(ierr); 1165 ierr = PetscHMapIClear(htWithAll);CHKERRQ(ierr); 1166 ierr = PetscSectionGetDof(cellCounts, v, &dof);CHKERRQ(ierr); 1167 ierr = PetscSectionGetOffset(cellCounts, v, &off);CHKERRQ(ierr); 1168 ierr = PetscSectionGetDof(pointCounts, v, &Np);CHKERRQ(ierr); 1169 ierr = PetscSectionGetOffset(pointCounts, v, &ooff);CHKERRQ(ierr); 1170 if (dof <= 0) continue; 1171 1172 for (k = 0; k < patch->nsubspaces; ++k) { 1173 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1174 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1175 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1176 PetscInt bs = patch->bs[k]; 1177 PetscInt goff; 1178 1179 for (i = off; i < off + dof; ++i) { 1180 /* Reconstruct mapping of global-to-local on this patch. */ 1181 const PetscInt c = cellsArray[i]; 1182 PetscInt cell = c; 1183 1184 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 1185 for (j = 0; j < nodesPerCell; ++j) { 1186 for (l = 0; l < bs; ++l) { 1187 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1188 const PetscInt localDof = dofsArray[key]; 1189 if (localDof >= 0) {ierr = PetscHMapISet(ht, globalDof, localDof);CHKERRQ(ierr);} 1190 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1191 const PetscInt localDofWithArtificial = dofsArrayWithArtificial[key]; 1192 if (localDofWithArtificial >= 0) { 1193 ierr = PetscHMapISet(htWithArtificial, globalDof, localDofWithArtificial);CHKERRQ(ierr); 1194 } 1195 } 1196 if (isNonlinear) { 1197 const PetscInt localDofWithAll = dofsArrayWithAll[key]; 1198 if (localDofWithAll >= 0) { 1199 ierr = PetscHMapISet(htWithAll, globalDof, localDofWithAll);CHKERRQ(ierr); 1200 } 1201 } 1202 key++; 1203 } 1204 } 1205 } 1206 1207 /* Shove it in the output data structure. */ 1208 ierr = PetscSectionGetOffset(gtolCounts, v, &goff);CHKERRQ(ierr); 1209 PetscHashIterBegin(ht, hi); 1210 while (!PetscHashIterAtEnd(ht, hi)) { 1211 PetscInt globalDof, localDof; 1212 1213 PetscHashIterGetKey(ht, hi, globalDof); 1214 PetscHashIterGetVal(ht, hi, localDof); 1215 if (globalDof >= 0) globalDofsArray[goff + localDof] = globalDof; 1216 PetscHashIterNext(ht, hi); 1217 } 1218 1219 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1220 ierr = PetscSectionGetOffset(gtolCountsWithArtificial, v, &goff);CHKERRQ(ierr); 1221 PetscHashIterBegin(htWithArtificial, hi); 1222 while (!PetscHashIterAtEnd(htWithArtificial, hi)) { 1223 PetscInt globalDof, localDof; 1224 PetscHashIterGetKey(htWithArtificial, hi, globalDof); 1225 PetscHashIterGetVal(htWithArtificial, hi, localDof); 1226 if (globalDof >= 0) globalDofsArrayWithArtificial[goff + localDof] = globalDof; 1227 PetscHashIterNext(htWithArtificial, hi); 1228 } 1229 } 1230 if (isNonlinear) { 1231 ierr = PetscSectionGetOffset(gtolCountsWithAll, v, &goff);CHKERRQ(ierr); 1232 PetscHashIterBegin(htWithAll, hi); 1233 while (!PetscHashIterAtEnd(htWithAll, hi)) { 1234 PetscInt globalDof, localDof; 1235 PetscHashIterGetKey(htWithAll, hi, globalDof); 1236 PetscHashIterGetVal(htWithAll, hi, localDof); 1237 if (globalDof >= 0) globalDofsArrayWithAll[goff + localDof] = globalDof; 1238 PetscHashIterNext(htWithAll, hi); 1239 } 1240 } 1241 1242 for (p = 0; p < Np; ++p) { 1243 const PetscInt point = pointsArray[ooff + p]; 1244 PetscInt globalDof, localDof; 1245 1246 ierr = PCPatchGetGlobalDofs(pc, patch->dofSection, k, patch->combined, point, NULL, &globalDof);CHKERRQ(ierr); 1247 ierr = PetscHMapIGet(ht, globalDof, &localDof);CHKERRQ(ierr); 1248 offsArray[(ooff + p)*Nf + k] = localDof; 1249 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1250 ierr = PetscHMapIGet(htWithArtificial, globalDof, &localDof);CHKERRQ(ierr); 1251 offsArrayWithArtificial[(ooff + p)*Nf + k] = localDof; 1252 } 1253 if (isNonlinear) { 1254 ierr = PetscHMapIGet(htWithAll, globalDof, &localDof);CHKERRQ(ierr); 1255 offsArrayWithAll[(ooff + p)*Nf + k] = localDof; 1256 } 1257 } 1258 } 1259 1260 ierr = PetscHSetIDestroy(&globalBcs);CHKERRQ(ierr); 1261 ierr = PetscHSetIDestroy(&ownedpts);CHKERRQ(ierr); 1262 ierr = PetscHSetIDestroy(&seenpts);CHKERRQ(ierr); 1263 ierr = PetscHSetIDestroy(&owneddofs);CHKERRQ(ierr); 1264 ierr = PetscHSetIDestroy(&seendofs);CHKERRQ(ierr); 1265 ierr = PetscHSetIDestroy(&artificialbcs);CHKERRQ(ierr); 1266 1267 /* At this point, we have a hash table ht built that maps globalDof -> localDof. 1268 We need to create the dof table laid out cellwise first, then by subspace, 1269 as the assembler assembles cell-wise and we need to stuff the different 1270 contributions of the different function spaces to the right places. So we loop 1271 over cells, then over subspaces. */ 1272 if (patch->nsubspaces > 1) { /* for nsubspaces = 1, data we need is already in dofsArray */ 1273 for (i = off; i < off + dof; ++i) { 1274 const PetscInt c = cellsArray[i]; 1275 PetscInt cell = c; 1276 1277 if (cellNumbering) {ierr = PetscSectionGetOffset(cellNumbering, c, &cell);CHKERRQ(ierr);} 1278 for (k = 0; k < patch->nsubspaces; ++k) { 1279 const PetscInt *cellNodeMap = patch->cellNodeMap[k]; 1280 PetscInt nodesPerCell = patch->nodesPerCell[k]; 1281 PetscInt subspaceOffset = patch->subspaceOffsets[k]; 1282 PetscInt bs = patch->bs[k]; 1283 1284 for (j = 0; j < nodesPerCell; ++j) { 1285 for (l = 0; l < bs; ++l) { 1286 const PetscInt globalDof = cellNodeMap[cell*nodesPerCell + j]*bs + l + subspaceOffset; 1287 PetscInt localDof; 1288 1289 ierr = PetscHMapIGet(ht, globalDof, &localDof);CHKERRQ(ierr); 1290 /* If it's not in the hash table, i.e. is a BC dof, 1291 then the PetscHSetIMap above gives -1, which matches 1292 exactly the convention for PETSc's matrix assembly to 1293 ignore the dof. So we don't need to do anything here */ 1294 asmArray[asmKey] = localDof; 1295 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1296 ierr = PetscHMapIGet(htWithArtificial, globalDof, &localDof);CHKERRQ(ierr); 1297 asmArrayWithArtificial[asmKey] = localDof; 1298 } 1299 if (isNonlinear) { 1300 ierr = PetscHMapIGet(htWithAll, globalDof, &localDof);CHKERRQ(ierr); 1301 asmArrayWithAll[asmKey] = localDof; 1302 } 1303 asmKey++; 1304 } 1305 } 1306 } 1307 } 1308 } 1309 } 1310 if (1 == patch->nsubspaces) { 1311 ierr = PetscMemcpy(asmArray, dofsArray, numDofs * sizeof(PetscInt));CHKERRQ(ierr); 1312 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1313 ierr = PetscMemcpy(asmArrayWithArtificial, dofsArrayWithArtificial, numDofs * sizeof(PetscInt));CHKERRQ(ierr); 1314 } 1315 if (isNonlinear) { 1316 ierr = PetscMemcpy(asmArrayWithAll, dofsArrayWithAll, numDofs * sizeof(PetscInt));CHKERRQ(ierr); 1317 } 1318 } 1319 1320 ierr = PetscHMapIDestroy(&ht);CHKERRQ(ierr); 1321 ierr = PetscHMapIDestroy(&htWithArtificial);CHKERRQ(ierr); 1322 ierr = PetscHMapIDestroy(&htWithAll);CHKERRQ(ierr); 1323 ierr = ISRestoreIndices(cells, &cellsArray);CHKERRQ(ierr); 1324 ierr = ISRestoreIndices(points, &pointsArray);CHKERRQ(ierr); 1325 ierr = PetscFree(dofsArray);CHKERRQ(ierr); 1326 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1327 ierr = PetscFree(dofsArrayWithArtificial);CHKERRQ(ierr); 1328 } 1329 if (isNonlinear) { 1330 ierr = PetscFree(dofsArrayWithAll);CHKERRQ(ierr); 1331 } 1332 /* Create placeholder section for map from points to patch dofs */ 1333 ierr = PetscSectionCreate(PETSC_COMM_SELF, &patch->patchSection);CHKERRQ(ierr); 1334 ierr = PetscSectionSetNumFields(patch->patchSection, patch->nsubspaces);CHKERRQ(ierr); 1335 if (patch->combined) { 1336 PetscInt numFields; 1337 ierr = PetscSectionGetNumFields(patch->dofSection[0], &numFields);CHKERRQ(ierr); 1338 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); 1339 ierr = PetscSectionGetChart(patch->dofSection[0], &pStart, &pEnd);CHKERRQ(ierr); 1340 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1341 for (p = pStart; p < pEnd; ++p) { 1342 PetscInt dof, fdof, f; 1343 1344 ierr = PetscSectionGetDof(patch->dofSection[0], p, &dof);CHKERRQ(ierr); 1345 ierr = PetscSectionSetDof(patch->patchSection, p, dof);CHKERRQ(ierr); 1346 for (f = 0; f < patch->nsubspaces; ++f) { 1347 ierr = PetscSectionGetFieldDof(patch->dofSection[0], p, f, &fdof);CHKERRQ(ierr); 1348 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1349 } 1350 } 1351 } else { 1352 PetscInt pStartf, pEndf, f; 1353 pStart = PETSC_MAX_INT; 1354 pEnd = PETSC_MIN_INT; 1355 for (f = 0; f < patch->nsubspaces; ++f) { 1356 ierr = PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);CHKERRQ(ierr); 1357 pStart = PetscMin(pStart, pStartf); 1358 pEnd = PetscMax(pEnd, pEndf); 1359 } 1360 ierr = PetscSectionSetChart(patch->patchSection, pStart, pEnd);CHKERRQ(ierr); 1361 for (f = 0; f < patch->nsubspaces; ++f) { 1362 ierr = PetscSectionGetChart(patch->dofSection[f], &pStartf, &pEndf);CHKERRQ(ierr); 1363 for (p = pStartf; p < pEndf; ++p) { 1364 PetscInt fdof; 1365 ierr = PetscSectionGetDof(patch->dofSection[f], p, &fdof);CHKERRQ(ierr); 1366 ierr = PetscSectionAddDof(patch->patchSection, p, fdof);CHKERRQ(ierr); 1367 ierr = PetscSectionSetFieldDof(patch->patchSection, p, f, fdof);CHKERRQ(ierr); 1368 } 1369 } 1370 } 1371 ierr = PetscSectionSetUp(patch->patchSection);CHKERRQ(ierr); 1372 ierr = PetscSectionSetUseFieldOffsets(patch->patchSection, PETSC_TRUE);CHKERRQ(ierr); 1373 /* Replace cell indices with firedrake-numbered ones. */ 1374 ierr = ISGeneralSetIndices(cells, numCells, (const PetscInt *) newCellsArray, PETSC_OWN_POINTER);CHKERRQ(ierr); 1375 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofs, globalDofsArray, PETSC_OWN_POINTER, &patch->gtol);CHKERRQ(ierr); 1376 ierr = PetscObjectSetName((PetscObject) patch->gtol, "Global Indices");CHKERRQ(ierr); 1377 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_g2l_view", patch->classname);CHKERRQ(ierr); 1378 ierr = PetscSectionViewFromOptions(patch->gtolCounts, (PetscObject) pc, option);CHKERRQ(ierr); 1379 ierr = ISViewFromOptions(patch->gtol, (PetscObject) pc, option);CHKERRQ(ierr); 1380 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArray, PETSC_OWN_POINTER, &patch->dofs);CHKERRQ(ierr); 1381 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArray, PETSC_OWN_POINTER, &patch->offs);CHKERRQ(ierr); 1382 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1383 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithArtificial, globalDofsArrayWithArtificial, PETSC_OWN_POINTER, &patch->gtolWithArtificial);CHKERRQ(ierr); 1384 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithArtificial, PETSC_OWN_POINTER, &patch->dofsWithArtificial);CHKERRQ(ierr); 1385 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithArtificial, PETSC_OWN_POINTER, &patch->offsWithArtificial);CHKERRQ(ierr); 1386 } 1387 if (isNonlinear) { 1388 ierr = ISCreateGeneral(PETSC_COMM_SELF, numGlobalDofsWithAll, globalDofsArrayWithAll, PETSC_OWN_POINTER, &patch->gtolWithAll);CHKERRQ(ierr); 1389 ierr = ISCreateGeneral(PETSC_COMM_SELF, numDofs, asmArrayWithAll, PETSC_OWN_POINTER, &patch->dofsWithAll);CHKERRQ(ierr); 1390 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPoints*Nf, offsArrayWithAll, PETSC_OWN_POINTER, &patch->offsWithAll);CHKERRQ(ierr); 1391 } 1392 PetscFunctionReturn(0); 1393 } 1394 1395 static PetscErrorCode PCPatchZeroFillMatrix_Private(Mat mat, const PetscInt ncell, const PetscInt ndof, const PetscInt *dof) 1396 { 1397 PetscScalar *values = NULL; 1398 PetscInt rows, c, i; 1399 PetscErrorCode ierr; 1400 1401 PetscFunctionBegin; 1402 ierr = PetscCalloc1(ndof*ndof, &values);CHKERRQ(ierr); 1403 for (c = 0; c < ncell; ++c) { 1404 const PetscInt *idx = &dof[ndof*c]; 1405 ierr = MatSetValues(mat, ndof, idx, ndof, idx, values, INSERT_VALUES);CHKERRQ(ierr); 1406 } 1407 ierr = MatGetLocalSize(mat, &rows, NULL);CHKERRQ(ierr); 1408 for (i = 0; i < rows; ++i) { 1409 ierr = MatSetValues(mat, 1, &i, 1, &i, values, INSERT_VALUES);CHKERRQ(ierr); 1410 } 1411 ierr = MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1412 ierr = MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1413 ierr = PetscFree(values);CHKERRQ(ierr); 1414 PetscFunctionReturn(0); 1415 } 1416 1417 static PetscErrorCode PCPatchCreateMatrix_Private(PC pc, PetscInt point, Mat *mat, PetscBool withArtificial) 1418 { 1419 PC_PATCH *patch = (PC_PATCH *) pc->data; 1420 Vec x, y; 1421 PetscBool flg; 1422 PetscInt csize, rsize; 1423 const char *prefix = NULL; 1424 PetscErrorCode ierr; 1425 1426 PetscFunctionBegin; 1427 if(withArtificial) { 1428 /* would be nice if we could create a rectangular matrix of size numDofsWithArtificial x numDofs here */ 1429 x = patch->patchRHSWithArtificial[point]; 1430 y = patch->patchRHSWithArtificial[point]; 1431 } else { 1432 x = patch->patchRHS[point]; 1433 y = patch->patchUpdate[point]; 1434 } 1435 1436 ierr = VecGetSize(x, &csize);CHKERRQ(ierr); 1437 ierr = VecGetSize(y, &rsize);CHKERRQ(ierr); 1438 ierr = MatCreate(PETSC_COMM_SELF, mat);CHKERRQ(ierr); 1439 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1440 ierr = MatSetOptionsPrefix(*mat, prefix);CHKERRQ(ierr); 1441 ierr = MatAppendOptionsPrefix(*mat, "pc_patch_sub_");CHKERRQ(ierr); 1442 if (patch->sub_mat_type) {ierr = MatSetType(*mat, patch->sub_mat_type);CHKERRQ(ierr);} 1443 else if (!patch->sub_mat_type) {ierr = MatSetType(*mat, MATDENSE);CHKERRQ(ierr);} 1444 ierr = MatSetSizes(*mat, rsize, csize, rsize, csize);CHKERRQ(ierr); 1445 ierr = PetscObjectTypeCompare((PetscObject) *mat, MATDENSE, &flg);CHKERRQ(ierr); 1446 if (!flg) {ierr = PetscObjectTypeCompare((PetscObject)*mat, MATSEQDENSE, &flg);CHKERRQ(ierr);} 1447 /* Sparse patch matrices */ 1448 if (!flg) { 1449 PetscBT bt; 1450 PetscInt *dnnz = NULL; 1451 const PetscInt *dofsArray = NULL; 1452 PetscInt pStart, pEnd, ncell, offset, c, i, j; 1453 1454 if(withArtificial) { 1455 ierr = ISGetIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 1456 } else { 1457 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1458 } 1459 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1460 point += pStart; 1461 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1462 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1463 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1464 ierr = PetscCalloc1(rsize, &dnnz);CHKERRQ(ierr); 1465 ierr = PetscLogEventBegin(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1466 /* XXX: This uses N^2 bits to store the sparsity pattern on a 1467 * patch. This is probably OK if the patches are not too big, 1468 * but could use quite a bit of memory for planes in 3D. 1469 * Should we switch based on the value of rsize to a 1470 * hash-table (slower, but more memory efficient) approach? */ 1471 ierr = PetscBTCreate(rsize*rsize, &bt);CHKERRQ(ierr); 1472 for (c = 0; c < ncell; ++c) { 1473 const PetscInt *idx = dofsArray + (offset + c)*patch->totalDofsPerCell; 1474 for (i = 0; i < patch->totalDofsPerCell; ++i) { 1475 const PetscInt row = idx[i]; 1476 if (row < 0) continue; 1477 for (j = 0; j < patch->totalDofsPerCell; ++j) { 1478 const PetscInt col = idx[j]; 1479 const PetscInt key = row*rsize + col; 1480 if (col < 0) continue; 1481 if (!PetscBTLookupSet(bt, key)) ++dnnz[row]; 1482 } 1483 } 1484 } 1485 ierr = PetscBTDestroy(&bt);CHKERRQ(ierr); 1486 ierr = MatXAIJSetPreallocation(*mat, 1, dnnz, NULL, NULL, NULL);CHKERRQ(ierr); 1487 ierr = PetscFree(dnnz);CHKERRQ(ierr); 1488 ierr = PCPatchZeroFillMatrix_Private(*mat, ncell, patch->totalDofsPerCell, &dofsArray[offset*patch->totalDofsPerCell]);CHKERRQ(ierr); 1489 ierr = PetscLogEventEnd(PC_Patch_Prealloc, pc, 0, 0, 0);CHKERRQ(ierr); 1490 if(withArtificial) { 1491 ierr = ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 1492 } else { 1493 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1494 } 1495 } 1496 ierr = MatSetUp(*mat);CHKERRQ(ierr); 1497 PetscFunctionReturn(0); 1498 } 1499 1500 static PetscErrorCode PCPatchComputeFunction_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Vec F, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx) 1501 { 1502 PC_PATCH *patch = (PC_PATCH *) pc->data; 1503 DM dm; 1504 PetscSection s; 1505 const PetscInt *parray, *oarray; 1506 PetscInt Nf = patch->nsubspaces, Np, poff, p, f; 1507 PetscErrorCode ierr; 1508 1509 PetscFunctionBegin; 1510 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1511 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1512 /* Set offset into patch */ 1513 ierr = PetscSectionGetDof(patch->pointCounts, patchNum, &Np);CHKERRQ(ierr); 1514 ierr = PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);CHKERRQ(ierr); 1515 ierr = ISGetIndices(patch->points, &parray);CHKERRQ(ierr); 1516 ierr = ISGetIndices(patch->offs, &oarray);CHKERRQ(ierr); 1517 for (f = 0; f < Nf; ++f) { 1518 for (p = 0; p < Np; ++p) { 1519 const PetscInt point = parray[poff+p]; 1520 PetscInt dof; 1521 1522 ierr = PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);CHKERRQ(ierr); 1523 ierr = PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr); 1524 if (patch->nsubspaces == 1) {ierr = PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr);} 1525 else {ierr = PetscSectionSetOffset(patch->patchSection, point, -1);CHKERRQ(ierr);} 1526 } 1527 } 1528 ierr = ISRestoreIndices(patch->points, &parray);CHKERRQ(ierr); 1529 ierr = ISRestoreIndices(patch->offs, &oarray);CHKERRQ(ierr); 1530 if (patch->viewSection) {ierr = ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);CHKERRQ(ierr);} 1531 ierr = DMPlexComputeResidual_Patch_Internal(pc->dm, patch->patchSection, cellIS, 0.0, x, NULL, F, ctx);CHKERRQ(ierr); 1532 PetscFunctionReturn(0); 1533 } 1534 1535 PetscErrorCode PCPatchComputeFunction_Internal(PC pc, Vec x, Vec F, PetscInt point) 1536 { 1537 PC_PATCH *patch = (PC_PATCH *) pc->data; 1538 const PetscInt *dofsArray; 1539 const PetscInt *dofsArrayWithAll; 1540 const PetscInt *cellsArray; 1541 PetscInt ncell, offset, pStart, pEnd; 1542 PetscErrorCode ierr; 1543 1544 PetscFunctionBegin; 1545 ierr = PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1546 if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n"); 1547 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1548 ierr = ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 1549 ierr = ISGetIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1550 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1551 1552 point += pStart; 1553 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1554 1555 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1556 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1557 if (ncell <= 0) { 1558 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1559 PetscFunctionReturn(0); 1560 } 1561 PetscStackPush("PCPatch user callback"); 1562 /* Cannot reuse the same IS because the geometry info is being cached in it */ 1563 ierr = ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);CHKERRQ(ierr); 1564 ierr = patch->usercomputef(pc, point, x, F, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, 1565 dofsArrayWithAll + offset*patch->totalDofsPerCell, 1566 patch->usercomputefctx);CHKERRQ(ierr); 1567 PetscStackPop; 1568 ierr = ISDestroy(&patch->cellIS);CHKERRQ(ierr); 1569 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1570 ierr = ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 1571 ierr = ISRestoreIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1572 if (patch->viewMatrix) { 1573 char name[PETSC_MAX_PATH_LEN]; 1574 1575 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch vector for Point %D", point);CHKERRQ(ierr); 1576 ierr = PetscObjectSetName((PetscObject) F, name);CHKERRQ(ierr); 1577 ierr = ObjectView((PetscObject) F, patch->viewerMatrix, patch->formatMatrix);CHKERRQ(ierr); 1578 } 1579 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1580 PetscFunctionReturn(0); 1581 } 1582 1583 static PetscErrorCode PCPatchComputeOperator_DMPlex_Private(PC pc, PetscInt patchNum, Vec x, Mat J, IS cellIS, PetscInt n, const PetscInt *l2p, const PetscInt *l2pWithAll, void *ctx) 1584 { 1585 PC_PATCH *patch = (PC_PATCH *) pc->data; 1586 DM dm; 1587 PetscSection s; 1588 const PetscInt *parray, *oarray; 1589 PetscInt Nf = patch->nsubspaces, Np, poff, p, f; 1590 PetscErrorCode ierr; 1591 1592 PetscFunctionBegin; 1593 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1594 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1595 /* Set offset into patch */ 1596 ierr = PetscSectionGetDof(patch->pointCounts, patchNum, &Np);CHKERRQ(ierr); 1597 ierr = PetscSectionGetOffset(patch->pointCounts, patchNum, &poff);CHKERRQ(ierr); 1598 ierr = ISGetIndices(patch->points, &parray);CHKERRQ(ierr); 1599 ierr = ISGetIndices(patch->offs, &oarray);CHKERRQ(ierr); 1600 for (f = 0; f < Nf; ++f) { 1601 for (p = 0; p < Np; ++p) { 1602 const PetscInt point = parray[poff+p]; 1603 PetscInt dof; 1604 1605 ierr = PetscSectionGetFieldDof(patch->patchSection, point, f, &dof);CHKERRQ(ierr); 1606 ierr = PetscSectionSetFieldOffset(patch->patchSection, point, f, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr); 1607 if (patch->nsubspaces == 1) {ierr = PetscSectionSetOffset(patch->patchSection, point, oarray[(poff+p)*Nf+f]);CHKERRQ(ierr);} 1608 else {ierr = PetscSectionSetOffset(patch->patchSection, point, -1);CHKERRQ(ierr);} 1609 } 1610 } 1611 ierr = ISRestoreIndices(patch->points, &parray);CHKERRQ(ierr); 1612 ierr = ISRestoreIndices(patch->offs, &oarray);CHKERRQ(ierr); 1613 if (patch->viewSection) {ierr = ObjectView((PetscObject) patch->patchSection, patch->viewerSection, patch->formatSection);CHKERRQ(ierr);} 1614 /* TODO Shut off MatViewFromOptions() in MatAssemblyEnd() here */ 1615 ierr = DMPlexComputeJacobian_Patch_Internal(pc->dm, patch->patchSection, patch->patchSection, cellIS, 0.0, 0.0, x, NULL, J, J, ctx);CHKERRQ(ierr); 1616 PetscFunctionReturn(0); 1617 } 1618 1619 PetscErrorCode PCPatchComputeOperator_Internal(PC pc, Vec x, Mat mat, PetscInt point, PetscBool withArtificial) 1620 { 1621 PC_PATCH *patch = (PC_PATCH *) pc->data; 1622 const PetscInt *dofsArray; 1623 const PetscInt *dofsArrayWithArtificial = NULL; 1624 const PetscInt *dofsArrayWithAll = NULL; 1625 const PetscInt *cellsArray; 1626 PetscInt ncell, offset, pStart, pEnd; 1627 PetscBool isNonlinear; 1628 PetscErrorCode ierr; 1629 1630 PetscFunctionBegin; 1631 ierr = PetscLogEventBegin(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1632 isNonlinear = patch->isNonlinear; 1633 if (!patch->usercomputeop) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call PCPatchSetComputeOperator() to set user callback\n"); 1634 if(withArtificial) { 1635 ierr = ISGetIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 1636 } else { 1637 ierr = ISGetIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1638 } 1639 if (isNonlinear) { 1640 ierr = ISGetIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 1641 } 1642 ierr = ISGetIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1643 ierr = PetscSectionGetChart(patch->cellCounts, &pStart, &pEnd);CHKERRQ(ierr); 1644 1645 point += pStart; 1646 if (point >= pEnd) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Operator point %D not in [%D, %D)\n", point, pStart, pEnd);CHKERRQ(ierr); 1647 1648 ierr = PetscSectionGetDof(patch->cellCounts, point, &ncell);CHKERRQ(ierr); 1649 ierr = PetscSectionGetOffset(patch->cellCounts, point, &offset);CHKERRQ(ierr); 1650 if (ncell <= 0) { 1651 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1652 PetscFunctionReturn(0); 1653 } 1654 PetscStackPush("PCPatch user callback"); 1655 /* Cannot reuse the same IS because the geometry info is being cached in it */ 1656 ierr = ISCreateGeneral(PETSC_COMM_SELF, ncell, cellsArray + offset, PETSC_USE_POINTER, &patch->cellIS);CHKERRQ(ierr); 1657 ierr = patch->usercomputeop(pc, point, x, mat, patch->cellIS, ncell*patch->totalDofsPerCell, dofsArray + offset*patch->totalDofsPerCell, dofsArrayWithAll ? dofsArrayWithAll + offset*patch->totalDofsPerCell : NULL, patch->usercomputeopctx);CHKERRQ(ierr); 1658 PetscStackPop; 1659 ierr = ISDestroy(&patch->cellIS);CHKERRQ(ierr); 1660 if(withArtificial) { 1661 ierr = ISRestoreIndices(patch->dofsWithArtificial, &dofsArray);CHKERRQ(ierr); 1662 } else { 1663 ierr = ISRestoreIndices(patch->dofs, &dofsArray);CHKERRQ(ierr); 1664 } 1665 if (isNonlinear) { 1666 ierr = ISRestoreIndices(patch->dofsWithAll, &dofsArrayWithAll);CHKERRQ(ierr); 1667 } 1668 ierr = ISRestoreIndices(patch->cells, &cellsArray);CHKERRQ(ierr); 1669 if (patch->viewMatrix) { 1670 char name[PETSC_MAX_PATH_LEN]; 1671 1672 ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "Patch matrix for Point %D", point);CHKERRQ(ierr); 1673 ierr = PetscObjectSetName((PetscObject) mat, name);CHKERRQ(ierr); 1674 ierr = ObjectView((PetscObject) mat, patch->viewerMatrix, patch->formatMatrix);CHKERRQ(ierr); 1675 } 1676 ierr = PetscLogEventEnd(PC_Patch_ComputeOp, pc, 0, 0, 0);CHKERRQ(ierr); 1677 PetscFunctionReturn(0); 1678 } 1679 1680 PetscErrorCode PCPatch_ScatterLocal_Private(PC pc, PetscInt p, Vec x, Vec y, InsertMode mode, ScatterMode scat, PatchScatterType scattertype) 1681 { 1682 PC_PATCH *patch = (PC_PATCH *) pc->data; 1683 const PetscScalar *xArray = NULL; 1684 PetscScalar *yArray = NULL; 1685 const PetscInt *gtolArray = NULL; 1686 PetscInt dof, offset, lidx; 1687 PetscErrorCode ierr; 1688 1689 PetscFunctionBeginHot; 1690 ierr = PetscLogEventBegin(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 1691 ierr = VecGetArrayRead(x, &xArray);CHKERRQ(ierr); 1692 ierr = VecGetArray(y, &yArray);CHKERRQ(ierr); 1693 if (scattertype == SCATTER_WITHARTIFICIAL) { 1694 ierr = PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);CHKERRQ(ierr); 1695 ierr = PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offset);CHKERRQ(ierr); 1696 ierr = ISGetIndices(patch->gtolWithArtificial, >olArray);CHKERRQ(ierr); 1697 } else if (scattertype == SCATTER_WITHALL) { 1698 ierr = PetscSectionGetDof(patch->gtolCountsWithAll, p, &dof);CHKERRQ(ierr); 1699 ierr = PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offset);CHKERRQ(ierr); 1700 ierr = ISGetIndices(patch->gtolWithAll, >olArray);CHKERRQ(ierr); 1701 } else { 1702 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 1703 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 1704 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1705 } 1706 if (mode == INSERT_VALUES && scat != SCATTER_FORWARD) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't insert if not scattering forward\n"); 1707 if (mode == ADD_VALUES && scat != SCATTER_REVERSE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can't add if not scattering reverse\n"); 1708 for (lidx = 0; lidx < dof; ++lidx) { 1709 const PetscInt gidx = gtolArray[offset+lidx]; 1710 1711 if (mode == INSERT_VALUES) yArray[lidx] = xArray[gidx]; /* Forward */ 1712 else yArray[gidx] += xArray[lidx]; /* Reverse */ 1713 } 1714 if (scattertype == SCATTER_WITHARTIFICIAL) { 1715 ierr = ISRestoreIndices(patch->gtolWithArtificial, >olArray);CHKERRQ(ierr); 1716 } else if (scattertype == SCATTER_WITHALL) { 1717 ierr = ISRestoreIndices(patch->gtolWithAll, >olArray);CHKERRQ(ierr); 1718 } else { 1719 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1720 } 1721 ierr = VecRestoreArrayRead(x, &xArray);CHKERRQ(ierr); 1722 ierr = VecRestoreArray(y, &yArray);CHKERRQ(ierr); 1723 ierr = PetscLogEventEnd(PC_Patch_Scatter, pc, 0, 0, 0);CHKERRQ(ierr); 1724 PetscFunctionReturn(0); 1725 } 1726 1727 static PetscErrorCode PCSetUp_PATCH_Linear(PC pc) 1728 { 1729 PC_PATCH *patch = (PC_PATCH *) pc->data; 1730 const char *prefix; 1731 PetscInt i; 1732 PetscErrorCode ierr; 1733 1734 PetscFunctionBegin; 1735 if (!pc->setupcalled) { 1736 ierr = PetscMalloc1(patch->npatch, &patch->solver);CHKERRQ(ierr); 1737 ierr = PCGetOptionsPrefix(pc, &prefix);CHKERRQ(ierr); 1738 for (i = 0; i < patch->npatch; ++i) { 1739 KSP ksp; 1740 PC subpc; 1741 1742 ierr = KSPCreate(PETSC_COMM_SELF, &ksp);CHKERRQ(ierr); 1743 ierr = KSPSetOptionsPrefix(ksp, prefix);CHKERRQ(ierr); 1744 ierr = KSPAppendOptionsPrefix(ksp, "sub_");CHKERRQ(ierr); 1745 ierr = PetscObjectIncrementTabLevel((PetscObject) ksp, (PetscObject) pc, 1);CHKERRQ(ierr); 1746 ierr = KSPGetPC(ksp, &subpc);CHKERRQ(ierr); 1747 ierr = PetscObjectIncrementTabLevel((PetscObject) subpc, (PetscObject) pc, 1);CHKERRQ(ierr); 1748 ierr = PetscLogObjectParent((PetscObject) pc, (PetscObject) ksp);CHKERRQ(ierr); 1749 patch->solver[i] = (PetscObject) ksp; 1750 } 1751 } 1752 if (patch->save_operators) { 1753 for (i = 0; i < patch->npatch; ++i) { 1754 ierr = MatZeroEntries(patch->mat[i]);CHKERRQ(ierr); 1755 ierr = PCPatchComputeOperator_Internal(pc, NULL, patch->mat[i], i, PETSC_FALSE);CHKERRQ(ierr); 1756 ierr = KSPSetOperators((KSP) patch->solver[i], patch->mat[i], patch->mat[i]);CHKERRQ(ierr); 1757 } 1758 } 1759 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1760 for (i = 0; i < patch->npatch; ++i) { 1761 /* Instead of padding patch->patchUpdate with zeros to get */ 1762 /* patch->patchUpdateWithArtificial and then multiplying with the matrix, */ 1763 /* just get rid of the columns that correspond to the dofs with */ 1764 /* artificial bcs. That's of course fairly inefficient, hopefully we */ 1765 /* can just assemble the rectangular matrix in the first place. */ 1766 Mat matSquare; 1767 IS rowis; 1768 PetscInt dof; 1769 1770 ierr = MatGetSize(patch->mat[i], &dof, NULL);CHKERRQ(ierr); 1771 if (dof == 0) { 1772 patch->matWithArtificial[i] = NULL; 1773 continue; 1774 } 1775 1776 ierr = PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);CHKERRQ(ierr); 1777 ierr = MatZeroEntries(matSquare);CHKERRQ(ierr); 1778 ierr = PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);CHKERRQ(ierr); 1779 1780 ierr = MatGetSize(matSquare, &dof, NULL);CHKERRQ(ierr); 1781 ierr = ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis); CHKERRQ(ierr); 1782 if(pc->setupcalled) { 1783 ierr = MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_REUSE_MATRIX, &patch->matWithArtificial[i]); CHKERRQ(ierr); 1784 } else { 1785 ierr = MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &patch->matWithArtificial[i]); CHKERRQ(ierr); 1786 } 1787 ierr = ISDestroy(&rowis); CHKERRQ(ierr); 1788 ierr = MatDestroy(&matSquare);CHKERRQ(ierr); 1789 } 1790 } 1791 PetscFunctionReturn(0); 1792 } 1793 1794 static PetscErrorCode PCSetUp_PATCH(PC pc) 1795 { 1796 PC_PATCH *patch = (PC_PATCH *) pc->data; 1797 PetscInt i; 1798 PetscBool isNonlinear; 1799 PetscErrorCode ierr; 1800 1801 PetscFunctionBegin; 1802 if (!pc->setupcalled) { 1803 PetscInt pStart, pEnd, p; 1804 PetscInt localSize; 1805 1806 ierr = PetscLogEventBegin(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 1807 1808 isNonlinear = patch->isNonlinear; 1809 if (!patch->nsubspaces) { 1810 DM dm; 1811 PetscDS prob; 1812 PetscSection s; 1813 PetscInt cStart, cEnd, c, Nf, f, numGlobalBcs = 0, *globalBcs, *Nb, totNb = 0, **cellDofs; 1814 1815 ierr = PCGetDM(pc, &dm);CHKERRQ(ierr); 1816 if (!dm) SETERRQ(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Must set DM for PCPATCH or call PCPatchSetDiscretisationInfo()"); 1817 ierr = DMGetDefaultSection(dm, &s);CHKERRQ(ierr); 1818 ierr = PetscSectionGetNumFields(s, &Nf);CHKERRQ(ierr); 1819 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr); 1820 for (p = pStart; p < pEnd; ++p) { 1821 PetscInt cdof; 1822 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 1823 numGlobalBcs += cdof; 1824 } 1825 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1826 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 1827 ierr = PetscMalloc3(Nf, &Nb, Nf, &cellDofs, numGlobalBcs, &globalBcs);CHKERRQ(ierr); 1828 for (f = 0; f < Nf; ++f) { 1829 PetscFE fe; 1830 PetscDualSpace sp; 1831 PetscInt cdoff = 0; 1832 1833 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1834 /* ierr = PetscFEGetNumComponents(fe, &Nc[f]);CHKERRQ(ierr); */ 1835 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 1836 ierr = PetscDualSpaceGetDimension(sp, &Nb[f]);CHKERRQ(ierr); 1837 totNb += Nb[f]; 1838 1839 ierr = PetscMalloc1((cEnd-cStart)*Nb[f], &cellDofs[f]);CHKERRQ(ierr); 1840 for (c = cStart; c < cEnd; ++c) { 1841 PetscInt *closure = NULL; 1842 PetscInt clSize = 0, cl; 1843 1844 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1845 for (cl = 0; cl < clSize*2; cl += 2) { 1846 const PetscInt p = closure[cl]; 1847 PetscInt fdof, d, foff; 1848 1849 ierr = PetscSectionGetFieldDof(s, p, f, &fdof);CHKERRQ(ierr); 1850 ierr = PetscSectionGetFieldOffset(s, p, f, &foff);CHKERRQ(ierr); 1851 for (d = 0; d < fdof; ++d, ++cdoff) cellDofs[f][cdoff] = foff + d; 1852 } 1853 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 1854 } 1855 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]); 1856 } 1857 numGlobalBcs = 0; 1858 for (p = pStart; p < pEnd; ++p) { 1859 const PetscInt *ind; 1860 PetscInt off, cdof, d; 1861 1862 ierr = PetscSectionGetOffset(s, p, &off);CHKERRQ(ierr); 1863 ierr = PetscSectionGetConstraintDof(s, p, &cdof);CHKERRQ(ierr); 1864 ierr = PetscSectionGetConstraintIndices(s, p, &ind);CHKERRQ(ierr); 1865 for (d = 0; d < cdof; ++d) globalBcs[numGlobalBcs++] = off + ind[d]; 1866 } 1867 1868 ierr = PCPatchSetDiscretisationInfoCombined(pc, dm, Nb, (const PetscInt **) cellDofs, numGlobalBcs, globalBcs, numGlobalBcs, globalBcs);CHKERRQ(ierr); 1869 for (f = 0; f < Nf; ++f) { 1870 ierr = PetscFree(cellDofs[f]);CHKERRQ(ierr); 1871 } 1872 ierr = PetscFree3(Nb, cellDofs, globalBcs);CHKERRQ(ierr); 1873 ierr = PCPatchSetComputeFunction(pc, PCPatchComputeFunction_DMPlex_Private, NULL);CHKERRQ(ierr); 1874 ierr = PCPatchSetComputeOperator(pc, PCPatchComputeOperator_DMPlex_Private, NULL);CHKERRQ(ierr); 1875 } 1876 1877 localSize = patch->subspaceOffsets[patch->nsubspaces]; 1878 ierr = VecCreateSeq(PETSC_COMM_SELF, localSize, &patch->localRHS);CHKERRQ(ierr); 1879 ierr = VecSetUp(patch->localRHS);CHKERRQ(ierr); 1880 ierr = VecDuplicate(patch->localRHS, &patch->localUpdate);CHKERRQ(ierr); 1881 ierr = PCPatchCreateCellPatches(pc);CHKERRQ(ierr); 1882 ierr = PCPatchCreateCellPatchDiscretisationInfo(pc);CHKERRQ(ierr); 1883 1884 /* OK, now build the work vectors */ 1885 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, &pEnd);CHKERRQ(ierr); 1886 ierr = PetscMalloc1(patch->npatch, &patch->patchRHS);CHKERRQ(ierr); 1887 ierr = PetscMalloc1(patch->npatch, &patch->patchUpdate);CHKERRQ(ierr); 1888 1889 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1890 ierr = PetscMalloc1(patch->npatch, &patch->patchRHSWithArtificial);CHKERRQ(ierr); 1891 ierr = PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithArtificial);CHKERRQ(ierr); 1892 } 1893 if (isNonlinear) { 1894 ierr = PetscMalloc1(patch->npatch, &patch->dofMappingWithoutToWithAll);CHKERRQ(ierr); 1895 } 1896 for (p = pStart; p < pEnd; ++p) { 1897 PetscInt dof; 1898 1899 ierr = PetscSectionGetDof(patch->gtolCounts, p, &dof);CHKERRQ(ierr); 1900 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchRHS[p-pStart]);CHKERRQ(ierr); 1901 ierr = VecSetUp(patch->patchRHS[p-pStart]);CHKERRQ(ierr); 1902 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchUpdate[p-pStart]);CHKERRQ(ierr); 1903 ierr = VecSetUp(patch->patchUpdate[p-pStart]);CHKERRQ(ierr); 1904 if (patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 1905 const PetscInt *gtolArray, *gtolArrayWithArtificial = NULL; 1906 PetscInt numPatchDofs, offset; 1907 PetscInt numPatchDofsWithArtificial, offsetWithArtificial; 1908 PetscInt dofWithoutArtificialCounter = 0; 1909 PetscInt *patchWithoutArtificialToWithArtificialArray; 1910 1911 ierr = PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &dof);CHKERRQ(ierr); 1912 ierr = VecCreateSeq(PETSC_COMM_SELF, dof, &patch->patchRHSWithArtificial[p-pStart]);CHKERRQ(ierr); 1913 ierr = VecSetUp(patch->patchRHSWithArtificial[p-pStart]);CHKERRQ(ierr); 1914 1915 /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */ 1916 /* the index in the patch with all dofs */ 1917 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1918 1919 ierr = PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);CHKERRQ(ierr); 1920 1921 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 1922 ierr = ISGetIndices(patch->gtolWithArtificial, >olArrayWithArtificial);CHKERRQ(ierr); 1923 ierr = PetscSectionGetDof(patch->gtolCountsWithArtificial, p, &numPatchDofsWithArtificial);CHKERRQ(ierr); 1924 ierr = PetscSectionGetOffset(patch->gtolCountsWithArtificial, p, &offsetWithArtificial);CHKERRQ(ierr); 1925 1926 ierr = PetscMalloc1(numPatchDofs, &patchWithoutArtificialToWithArtificialArray);CHKERRQ(ierr); 1927 1928 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutArtificialToWithArtificialArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithArtificial[p-pStart]);CHKERRQ(ierr); 1929 if (numPatchDofs == 0) continue; 1930 for (i=0; i<numPatchDofsWithArtificial; i++) { 1931 if (gtolArrayWithArtificial[i+offsetWithArtificial] == gtolArray[offset+dofWithoutArtificialCounter]) { 1932 patchWithoutArtificialToWithArtificialArray[dofWithoutArtificialCounter] = i; 1933 dofWithoutArtificialCounter++; 1934 if (dofWithoutArtificialCounter == numPatchDofs) 1935 break; 1936 } 1937 } 1938 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1939 ierr = ISRestoreIndices(patch->gtolWithArtificial, >olArrayWithArtificial);CHKERRQ(ierr); 1940 } 1941 if (isNonlinear) { 1942 const PetscInt *gtolArray, *gtolArrayWithAll = NULL; 1943 PetscInt numPatchDofs, offset; 1944 PetscInt numPatchDofsWithAll, offsetWithAll; 1945 PetscInt dofWithoutAllCounter = 0; 1946 PetscInt *patchWithoutAllToWithAllArray; 1947 1948 /* Now build the mapping that for a dof in a patch WITHOUT dofs that have artificial bcs gives the */ 1949 /* the index in the patch with all dofs */ 1950 ierr = ISGetIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1951 1952 ierr = PetscSectionGetDof(patch->gtolCounts, p, &numPatchDofs);CHKERRQ(ierr); 1953 1954 ierr = PetscSectionGetOffset(patch->gtolCounts, p, &offset);CHKERRQ(ierr); 1955 ierr = ISGetIndices(patch->gtolWithAll, >olArrayWithAll);CHKERRQ(ierr); 1956 ierr = PetscSectionGetDof(patch->gtolCountsWithAll, p, &numPatchDofsWithAll);CHKERRQ(ierr); 1957 ierr = PetscSectionGetOffset(patch->gtolCountsWithAll, p, &offsetWithAll);CHKERRQ(ierr); 1958 1959 ierr = PetscMalloc1(numPatchDofs, &patchWithoutAllToWithAllArray);CHKERRQ(ierr); 1960 1961 ierr = ISCreateGeneral(PETSC_COMM_SELF, numPatchDofs, patchWithoutAllToWithAllArray, PETSC_OWN_POINTER, &patch->dofMappingWithoutToWithAll[p-pStart]);CHKERRQ(ierr); 1962 if (numPatchDofs == 0) continue; 1963 for (i=0; i<numPatchDofsWithAll; i++) { 1964 if (gtolArrayWithAll[i+offsetWithAll] == gtolArray[offset+dofWithoutAllCounter]) { 1965 patchWithoutAllToWithAllArray[dofWithoutAllCounter] = i; 1966 dofWithoutAllCounter++; 1967 if (dofWithoutAllCounter == numPatchDofs) 1968 break; 1969 } 1970 } 1971 ierr = ISRestoreIndices(patch->gtol, >olArray);CHKERRQ(ierr); 1972 ierr = ISRestoreIndices(patch->gtolWithAll, >olArrayWithAll);CHKERRQ(ierr); 1973 } 1974 } 1975 if (patch->save_operators) { 1976 ierr = PetscMalloc1(patch->npatch, &patch->mat);CHKERRQ(ierr); 1977 for (i = 0; i < patch->npatch; ++i) { 1978 ierr = PCPatchCreateMatrix_Private(pc, i, &patch->mat[i], PETSC_FALSE);CHKERRQ(ierr); 1979 } 1980 } 1981 ierr = PetscLogEventEnd(PC_Patch_CreatePatches, pc, 0, 0, 0);CHKERRQ(ierr); 1982 1983 /* If desired, calculate weights for dof multiplicity */ 1984 if (patch->partition_of_unity) { 1985 PetscScalar *input = NULL; 1986 PetscScalar *output = NULL; 1987 Vec global; 1988 1989 ierr = VecDuplicate(patch->localRHS, &patch->dof_weights);CHKERRQ(ierr); 1990 if(patch->local_composition_type == PC_COMPOSITE_ADDITIVE) { 1991 for (i = 0; i < patch->npatch; ++i) { 1992 PetscInt dof; 1993 1994 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &dof);CHKERRQ(ierr); 1995 if (dof <= 0) continue; 1996 ierr = VecSet(patch->patchRHS[i], 1.0);CHKERRQ(ierr); 1997 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchRHS[i], patch->dof_weights, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);CHKERRQ(ierr); 1998 } 1999 } else { 2000 /* multiplicative is actually only locally multiplicative and globally additive. need the pou where the mesh decomposition overlaps */ 2001 ierr = VecSet(patch->dof_weights, 1.0);CHKERRQ(ierr); 2002 } 2003 2004 VecDuplicate(patch->dof_weights, &global); 2005 VecSet(global, 0.); 2006 2007 ierr = VecGetArray(patch->dof_weights, &input);CHKERRQ(ierr); 2008 ierr = VecGetArray(global, &output);CHKERRQ(ierr); 2009 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);CHKERRQ(ierr); 2010 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, input, output, MPI_SUM);CHKERRQ(ierr); 2011 ierr = VecRestoreArray(patch->dof_weights, &input);CHKERRQ(ierr); 2012 ierr = VecRestoreArray(global, &output);CHKERRQ(ierr); 2013 2014 ierr = VecReciprocal(global);CHKERRQ(ierr); 2015 2016 ierr = VecGetArray(patch->dof_weights, &output);CHKERRQ(ierr); 2017 ierr = VecGetArray(global, &input);CHKERRQ(ierr); 2018 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, input, output);CHKERRQ(ierr); 2019 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, input, output);CHKERRQ(ierr); 2020 ierr = VecRestoreArray(patch->dof_weights, &output);CHKERRQ(ierr); 2021 ierr = VecRestoreArray(global, &input);CHKERRQ(ierr); 2022 ierr = VecDestroy(&global);CHKERRQ(ierr); 2023 } 2024 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE && patch->save_operators) { 2025 ierr = PetscMalloc1(patch->npatch, &patch->matWithArtificial);CHKERRQ(ierr); 2026 } 2027 } 2028 ierr = (*patch->setupsolver)(pc);CHKERRQ(ierr); 2029 PetscFunctionReturn(0); 2030 } 2031 2032 static PetscErrorCode PCApply_PATCH_Linear(PC pc, PetscInt i, Vec x, Vec y) 2033 { 2034 PC_PATCH *patch = (PC_PATCH *) pc->data; 2035 KSP ksp = (KSP) patch->solver[i]; 2036 PetscErrorCode ierr; 2037 2038 PetscFunctionBegin; 2039 if (!patch->save_operators) { 2040 Mat mat; 2041 2042 ierr = PCPatchCreateMatrix_Private(pc, i, &mat, PETSC_FALSE);CHKERRQ(ierr); 2043 /* Populate operator here. */ 2044 ierr = PCPatchComputeOperator_Internal(pc, NULL, mat, i, PETSC_FALSE);CHKERRQ(ierr); 2045 ierr = KSPSetOperators(ksp, mat, mat);CHKERRQ(ierr); 2046 /* Drop reference so the KSPSetOperators below will blow it away. */ 2047 ierr = MatDestroy(&mat);CHKERRQ(ierr); 2048 } 2049 ierr = PetscLogEventBegin(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 2050 if (!ksp->setfromoptionscalled) { 2051 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 2052 } 2053 ierr = KSPSolve(ksp, x, y);CHKERRQ(ierr); 2054 ierr = PetscLogEventEnd(PC_Patch_Solve, pc, 0, 0, 0);CHKERRQ(ierr); 2055 if (!patch->save_operators) { 2056 PC pc; 2057 ierr = KSPSetOperators(ksp, NULL, NULL);CHKERRQ(ierr); 2058 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 2059 /* Destroy PC context too, otherwise the factored matrix hangs around. */ 2060 ierr = PCReset(pc);CHKERRQ(ierr); 2061 } 2062 PetscFunctionReturn(0); 2063 } 2064 2065 static PetscErrorCode PCUpdateMultiplicative_PATCH_Linear(PC pc, PetscInt i, PetscInt pStart) 2066 { 2067 PC_PATCH *patch = (PC_PATCH *) pc->data; 2068 Mat multMat; 2069 PetscErrorCode ierr; 2070 2071 PetscFunctionBegin; 2072 2073 if (patch->save_operators) { 2074 multMat = patch->matWithArtificial[i]; 2075 } else { 2076 /*Very inefficient, hopefully we can just assemble the rectangular matrix in the first place.*/ 2077 Mat matSquare; 2078 PetscInt dof; 2079 IS rowis; 2080 ierr = PCPatchCreateMatrix_Private(pc, i, &matSquare, PETSC_TRUE);CHKERRQ(ierr); 2081 ierr = MatZeroEntries(matSquare);CHKERRQ(ierr); 2082 ierr = PCPatchComputeOperator_Internal(pc, NULL, matSquare, i, PETSC_TRUE);CHKERRQ(ierr); 2083 ierr = MatGetSize(matSquare, &dof, NULL);CHKERRQ(ierr); 2084 ierr = ISCreateStride(PETSC_COMM_SELF, dof, 0, 1, &rowis); CHKERRQ(ierr); 2085 ierr = MatCreateSubMatrix(matSquare, rowis, patch->dofMappingWithoutToWithArtificial[i], MAT_INITIAL_MATRIX, &multMat); CHKERRQ(ierr); 2086 ierr = MatDestroy(&matSquare);CHKERRQ(ierr); 2087 ierr = ISDestroy(&rowis); CHKERRQ(ierr); 2088 } 2089 ierr = MatMult(multMat, patch->patchUpdate[i], patch->patchRHSWithArtificial[i]); CHKERRQ(ierr); 2090 ierr = VecScale(patch->patchRHSWithArtificial[i], -1.0); CHKERRQ(ierr); 2091 ierr = PCPatch_ScatterLocal_Private(pc, i + pStart, patch->patchRHSWithArtificial[i], patch->localRHS, ADD_VALUES, SCATTER_REVERSE, SCATTER_WITHARTIFICIAL); CHKERRQ(ierr); 2092 if (!patch->save_operators) { 2093 ierr = MatDestroy(&multMat); CHKERRQ(ierr); 2094 } 2095 PetscFunctionReturn(0); 2096 } 2097 2098 static PetscErrorCode PCApply_PATCH(PC pc, Vec x, Vec y) 2099 { 2100 PC_PATCH *patch = (PC_PATCH *) pc->data; 2101 const PetscScalar *globalRHS = NULL; 2102 PetscScalar *localRHS = NULL; 2103 PetscScalar *globalUpdate = NULL; 2104 const PetscInt *bcNodes = NULL; 2105 PetscInt nsweep = patch->symmetrise_sweep ? 2 : 1; 2106 PetscInt start[2] = {0, 0}; 2107 PetscInt end[2] = {-1, -1}; 2108 const PetscInt inc[2] = {1, -1}; 2109 const PetscScalar *localUpdate; 2110 const PetscInt *iterationSet; 2111 PetscInt pStart, numBcs, n, sweep, bc, j; 2112 PetscErrorCode ierr; 2113 2114 PetscFunctionBegin; 2115 ierr = PetscLogEventBegin(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 2116 ierr = PetscOptionsPushGetViewerOff(PETSC_TRUE);CHKERRQ(ierr); 2117 /* start, end, inc have 2 entries to manage a second backward sweep if we symmetrize */ 2118 end[0] = patch->npatch; 2119 start[1] = patch->npatch-1; 2120 if (patch->user_patches) { 2121 ierr = ISGetLocalSize(patch->iterationSet, &end[0]);CHKERRQ(ierr); 2122 start[1] = end[0] - 1; 2123 ierr = ISGetIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr); 2124 } 2125 /* Scatter from global space into overlapped local spaces */ 2126 ierr = VecGetArrayRead(x, &globalRHS);CHKERRQ(ierr); 2127 ierr = VecGetArray(patch->localRHS, &localRHS);CHKERRQ(ierr); 2128 ierr = PetscSFBcastBegin(patch->defaultSF, MPIU_SCALAR, globalRHS, localRHS);CHKERRQ(ierr); 2129 ierr = PetscSFBcastEnd(patch->defaultSF, MPIU_SCALAR, globalRHS, localRHS);CHKERRQ(ierr); 2130 ierr = VecRestoreArrayRead(x, &globalRHS);CHKERRQ(ierr); 2131 ierr = VecRestoreArray(patch->localRHS, &localRHS);CHKERRQ(ierr); 2132 2133 ierr = VecSet(patch->localUpdate, 0.0);CHKERRQ(ierr); 2134 ierr = PetscSectionGetChart(patch->gtolCounts, &pStart, NULL);CHKERRQ(ierr); 2135 for (sweep = 0; sweep < nsweep; sweep++) { 2136 for (j = start[sweep]; j*inc[sweep] < end[sweep]*inc[sweep]; j += inc[sweep]) { 2137 PetscInt i = patch->user_patches ? iterationSet[j] : j; 2138 PetscInt start, len; 2139 2140 ierr = PetscSectionGetDof(patch->gtolCounts, i+pStart, &len);CHKERRQ(ierr); 2141 ierr = PetscSectionGetOffset(patch->gtolCounts, i+pStart, &start);CHKERRQ(ierr); 2142 /* TODO: Squash out these guys in the setup as well. */ 2143 if (len <= 0) continue; 2144 /* TODO: Do we need different scatters for X and Y? */ 2145 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->localRHS, patch->patchRHS[i], INSERT_VALUES, SCATTER_FORWARD, SCATTER_INTERIOR);CHKERRQ(ierr); 2146 ierr = (*patch->applysolver)(pc, i, patch->patchRHS[i], patch->patchUpdate[i]);CHKERRQ(ierr); 2147 ierr = PCPatch_ScatterLocal_Private(pc, i+pStart, patch->patchUpdate[i], patch->localUpdate, ADD_VALUES, SCATTER_REVERSE, SCATTER_INTERIOR);CHKERRQ(ierr); 2148 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 2149 ierr = (*patch->updatemultiplicative)(pc, i, pStart);CHKERRQ(ierr); 2150 } 2151 } 2152 } 2153 if (patch->user_patches) {ierr = ISRestoreIndices(patch->iterationSet, &iterationSet);CHKERRQ(ierr);} 2154 /* XXX: should we do this on the global vector? */ 2155 if (patch->partition_of_unity) { 2156 ierr = VecPointwiseMult(patch->localUpdate, patch->localUpdate, patch->dof_weights);CHKERRQ(ierr); 2157 } 2158 /* Now patch->localUpdate contains the solution of the patch solves, so we need to combine them all. */ 2159 ierr = VecSet(y, 0.0);CHKERRQ(ierr); 2160 ierr = VecGetArray(y, &globalUpdate);CHKERRQ(ierr); 2161 ierr = VecGetArrayRead(patch->localUpdate, &localUpdate);CHKERRQ(ierr); 2162 ierr = PetscSFReduceBegin(patch->defaultSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);CHKERRQ(ierr); 2163 ierr = PetscSFReduceEnd(patch->defaultSF, MPIU_SCALAR, localUpdate, globalUpdate, MPI_SUM);CHKERRQ(ierr); 2164 ierr = VecRestoreArrayRead(patch->localUpdate, &localUpdate);CHKERRQ(ierr); 2165 2166 /* Now we need to send the global BC values through */ 2167 ierr = VecGetArrayRead(x, &globalRHS);CHKERRQ(ierr); 2168 ierr = ISGetSize(patch->globalBcNodes, &numBcs);CHKERRQ(ierr); 2169 ierr = ISGetIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 2170 ierr = VecGetLocalSize(x, &n);CHKERRQ(ierr); 2171 for (bc = 0; bc < numBcs; ++bc) { 2172 const PetscInt idx = bcNodes[bc]; 2173 if (idx < n) globalUpdate[idx] = globalRHS[idx]; 2174 } 2175 2176 ierr = ISRestoreIndices(patch->globalBcNodes, &bcNodes);CHKERRQ(ierr); 2177 ierr = VecRestoreArrayRead(x, &globalRHS);CHKERRQ(ierr); 2178 ierr = VecRestoreArray(y, &globalUpdate);CHKERRQ(ierr); 2179 2180 ierr = PetscOptionsPopGetViewerOff();CHKERRQ(ierr); 2181 ierr = PetscLogEventEnd(PC_Patch_Apply, pc, 0, 0, 0);CHKERRQ(ierr); 2182 PetscFunctionReturn(0); 2183 } 2184 2185 static PetscErrorCode PCReset_PATCH_Linear(PC pc) 2186 { 2187 PC_PATCH *patch = (PC_PATCH *) pc->data; 2188 PetscInt i; 2189 PetscErrorCode ierr; 2190 2191 PetscFunctionBegin; 2192 if (patch->solver) { 2193 for (i = 0; i < patch->npatch; ++i) {ierr = KSPReset((KSP) patch->solver[i]);CHKERRQ(ierr);} 2194 } 2195 PetscFunctionReturn(0); 2196 } 2197 2198 static PetscErrorCode PCReset_PATCH(PC pc) 2199 { 2200 PC_PATCH *patch = (PC_PATCH *) pc->data; 2201 PetscInt i; 2202 PetscErrorCode ierr; 2203 2204 PetscFunctionBegin; 2205 /* TODO: Get rid of all these ifs */ 2206 ierr = PetscSFDestroy(&patch->defaultSF);CHKERRQ(ierr); 2207 ierr = PetscSectionDestroy(&patch->cellCounts);CHKERRQ(ierr); 2208 ierr = PetscSectionDestroy(&patch->pointCounts);CHKERRQ(ierr); 2209 ierr = PetscSectionDestroy(&patch->cellNumbering);CHKERRQ(ierr); 2210 ierr = PetscSectionDestroy(&patch->gtolCounts);CHKERRQ(ierr); 2211 ierr = ISDestroy(&patch->gtol);CHKERRQ(ierr); 2212 ierr = ISDestroy(&patch->cells);CHKERRQ(ierr); 2213 ierr = ISDestroy(&patch->points);CHKERRQ(ierr); 2214 ierr = ISDestroy(&patch->dofs);CHKERRQ(ierr); 2215 ierr = ISDestroy(&patch->offs);CHKERRQ(ierr); 2216 ierr = PetscSectionDestroy(&patch->patchSection);CHKERRQ(ierr); 2217 ierr = ISDestroy(&patch->ghostBcNodes);CHKERRQ(ierr); 2218 ierr = ISDestroy(&patch->globalBcNodes);CHKERRQ(ierr); 2219 ierr = PetscSectionDestroy(&patch->gtolCountsWithArtificial);CHKERRQ(ierr); 2220 ierr = ISDestroy(&patch->gtolWithArtificial);CHKERRQ(ierr); 2221 ierr = ISDestroy(&patch->dofsWithArtificial);CHKERRQ(ierr); 2222 ierr = ISDestroy(&patch->offsWithArtificial);CHKERRQ(ierr); 2223 ierr = PetscSectionDestroy(&patch->gtolCountsWithAll);CHKERRQ(ierr); 2224 ierr = ISDestroy(&patch->gtolWithAll);CHKERRQ(ierr); 2225 ierr = ISDestroy(&patch->dofsWithAll);CHKERRQ(ierr); 2226 ierr = ISDestroy(&patch->offsWithAll);CHKERRQ(ierr); 2227 2228 2229 if (patch->dofSection) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscSectionDestroy(&patch->dofSection[i]);CHKERRQ(ierr);} 2230 ierr = PetscFree(patch->dofSection);CHKERRQ(ierr); 2231 ierr = PetscFree(patch->bs);CHKERRQ(ierr); 2232 ierr = PetscFree(patch->nodesPerCell);CHKERRQ(ierr); 2233 if (patch->cellNodeMap) for (i = 0; i < patch->nsubspaces; i++) {ierr = PetscFree(patch->cellNodeMap[i]);CHKERRQ(ierr);} 2234 ierr = PetscFree(patch->cellNodeMap);CHKERRQ(ierr); 2235 ierr = PetscFree(patch->subspaceOffsets);CHKERRQ(ierr); 2236 2237 ierr = (*patch->resetsolver)(pc);CHKERRQ(ierr); 2238 2239 if (patch->subspaces_to_exclude) { 2240 PetscHSetIDestroy(&patch->subspaces_to_exclude); 2241 } 2242 2243 ierr = VecDestroy(&patch->localRHS);CHKERRQ(ierr); 2244 ierr = VecDestroy(&patch->localUpdate);CHKERRQ(ierr); 2245 if (patch->patchRHS) { 2246 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchRHS[i]);CHKERRQ(ierr);} 2247 ierr = PetscFree(patch->patchRHS);CHKERRQ(ierr); 2248 } 2249 if (patch->patchUpdate) { 2250 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchUpdate[i]);CHKERRQ(ierr);} 2251 ierr = PetscFree(patch->patchUpdate);CHKERRQ(ierr); 2252 } 2253 ierr = VecDestroy(&patch->dof_weights);CHKERRQ(ierr); 2254 if (patch->patch_dof_weights) { 2255 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patch_dof_weights[i]);CHKERRQ(ierr);} 2256 ierr = PetscFree(patch->patch_dof_weights);CHKERRQ(ierr); 2257 } 2258 if (patch->mat) { 2259 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->mat[i]);CHKERRQ(ierr);} 2260 ierr = PetscFree(patch->mat);CHKERRQ(ierr); 2261 } 2262 if (patch->matWithArtificial) { 2263 for (i = 0; i < patch->npatch; ++i) {ierr = MatDestroy(&patch->matWithArtificial[i]);CHKERRQ(ierr);} 2264 ierr = PetscFree(patch->matWithArtificial);CHKERRQ(ierr); 2265 } 2266 if (patch->patchRHSWithArtificial) { 2267 for (i = 0; i < patch->npatch; ++i) {ierr = VecDestroy(&patch->patchRHSWithArtificial[i]);CHKERRQ(ierr);} 2268 ierr = PetscFree(patch->patchRHSWithArtificial);CHKERRQ(ierr); 2269 } 2270 if(patch->dofMappingWithoutToWithArtificial) { 2271 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->dofMappingWithoutToWithArtificial[i]);CHKERRQ(ierr);} 2272 ierr = PetscFree(patch->dofMappingWithoutToWithArtificial);CHKERRQ(ierr); 2273 2274 } 2275 if(patch->dofMappingWithoutToWithAll) { 2276 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->dofMappingWithoutToWithAll[i]);CHKERRQ(ierr);} 2277 ierr = PetscFree(patch->dofMappingWithoutToWithAll);CHKERRQ(ierr); 2278 2279 } 2280 ierr = PetscFree(patch->sub_mat_type);CHKERRQ(ierr); 2281 if (patch->userIS) { 2282 for (i = 0; i < patch->npatch; ++i) {ierr = ISDestroy(&patch->userIS[i]);CHKERRQ(ierr);} 2283 ierr = PetscFree(patch->userIS);CHKERRQ(ierr); 2284 } 2285 patch->bs = 0; 2286 patch->cellNodeMap = NULL; 2287 patch->nsubspaces = 0; 2288 ierr = ISDestroy(&patch->iterationSet);CHKERRQ(ierr); 2289 2290 ierr = PetscViewerDestroy(&patch->viewerSection);CHKERRQ(ierr); 2291 PetscFunctionReturn(0); 2292 } 2293 2294 static PetscErrorCode PCDestroy_PATCH_Linear(PC pc) 2295 { 2296 PC_PATCH *patch = (PC_PATCH *) pc->data; 2297 PetscInt i; 2298 PetscErrorCode ierr; 2299 2300 PetscFunctionBegin; 2301 if (patch->solver) { 2302 for (i = 0; i < patch->npatch; ++i) {ierr = KSPDestroy((KSP *) &patch->solver[i]);CHKERRQ(ierr);} 2303 ierr = PetscFree(patch->solver);CHKERRQ(ierr); 2304 } 2305 PetscFunctionReturn(0); 2306 } 2307 2308 static PetscErrorCode PCDestroy_PATCH(PC pc) 2309 { 2310 PC_PATCH *patch = (PC_PATCH *) pc->data; 2311 PetscErrorCode ierr; 2312 2313 PetscFunctionBegin; 2314 ierr = PCReset_PATCH(pc);CHKERRQ(ierr); 2315 ierr = (*patch->destroysolver)(pc);CHKERRQ(ierr); 2316 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2317 PetscFunctionReturn(0); 2318 } 2319 2320 static PetscErrorCode PCSetFromOptions_PATCH(PetscOptionItems *PetscOptionsObject, PC pc) 2321 { 2322 PC_PATCH *patch = (PC_PATCH *) pc->data; 2323 PCPatchConstructType patchConstructionType = PC_PATCH_STAR; 2324 char sub_mat_type[PETSC_MAX_PATH_LEN]; 2325 char option[PETSC_MAX_PATH_LEN]; 2326 const char *prefix; 2327 PetscBool flg, dimflg, codimflg; 2328 MPI_Comm comm; 2329 PetscInt *ifields, nfields, k; 2330 PetscErrorCode ierr; 2331 PCCompositeType loctype = PC_COMPOSITE_ADDITIVE; 2332 2333 PetscFunctionBegin; 2334 ierr = PetscObjectGetComm((PetscObject) pc, &comm);CHKERRQ(ierr); 2335 ierr = PetscObjectGetOptionsPrefix((PetscObject) pc, &prefix);CHKERRQ(ierr); 2336 ierr = PetscOptionsHead(PetscOptionsObject, "Patch solver options");CHKERRQ(ierr); 2337 2338 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_save_operators", patch->classname); 2339 ierr = PetscOptionsBool(option, "Store all patch operators for lifetime of object?", "PCPatchSetSaveOperators", patch->save_operators, &patch->save_operators, &flg);CHKERRQ(ierr); 2340 2341 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_partition_of_unity", patch->classname); 2342 ierr = PetscOptionsBool(option, "Weight contributions by dof multiplicity?", "PCPatchSetPartitionOfUnity", patch->partition_of_unity, &patch->partition_of_unity, &flg);CHKERRQ(ierr); 2343 2344 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_local_type", patch->classname); 2345 ierr = PetscOptionsEnum(option,"Type of local solver composition (additive or multiplicative)","PCPatchSetLocalComposition",PCCompositeTypes,(PetscEnum)loctype,(PetscEnum*)&loctype,&flg);CHKERRQ(ierr); 2346 if(flg) { ierr = PCPatchSetLocalComposition(pc, loctype);CHKERRQ(ierr);} 2347 2348 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_dim", patch->classname); 2349 ierr = PetscOptionsInt(option, "What dimension of mesh point to construct patches by? (0 = vertices)", "PCPATCH", patch->dim, &patch->dim, &dimflg);CHKERRQ(ierr); 2350 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_codim", patch->classname); 2351 ierr = PetscOptionsInt(option, "What co-dimension of mesh point to construct patches by? (0 = cells)", "PCPATCH", patch->codim, &patch->codim, &codimflg);CHKERRQ(ierr); 2352 if (dimflg && codimflg) {SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Can only set one of dimension or co-dimension");CHKERRQ(ierr);} 2353 2354 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_construct_type", patch->classname); 2355 ierr = PetscOptionsEnum(option, "How should the patches be constructed?", "PCPatchSetConstructType", PCPatchConstructTypes, (PetscEnum) patchConstructionType, (PetscEnum *) &patchConstructionType, &flg);CHKERRQ(ierr); 2356 if (flg) {ierr = PCPatchSetConstructType(pc, patchConstructionType, NULL, NULL);CHKERRQ(ierr);} 2357 2358 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_vanka_dim", patch->classname); 2359 ierr = PetscOptionsInt(option, "Topological dimension of entities for Vanka to ignore", "PCPATCH", patch->vankadim, &patch->vankadim, &flg);CHKERRQ(ierr); 2360 2361 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_ignore_dim", patch->classname); 2362 ierr = PetscOptionsInt(option, "Topological dimension of entities for completion to ignore", "PCPATCH", patch->ignoredim, &patch->ignoredim, &flg);CHKERRQ(ierr); 2363 2364 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_sub_mat_type", patch->classname); 2365 ierr = PetscOptionsFList(option, "Matrix type for patch solves", "PCPatchSetSubMatType", MatList, NULL, sub_mat_type, PETSC_MAX_PATH_LEN, &flg);CHKERRQ(ierr); 2366 if (flg) {ierr = PCPatchSetSubMatType(pc, sub_mat_type);CHKERRQ(ierr);} 2367 2368 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_symmetrise_sweep", patch->classname); 2369 ierr = PetscOptionsBool(option, "Go start->end, end->start?", "PCPATCH", patch->symmetrise_sweep, &patch->symmetrise_sweep, &flg);CHKERRQ(ierr); 2370 2371 /* If the user has set the number of subspaces, use that for the buffer size, 2372 otherwise use a large number */ 2373 if (patch->nsubspaces <= 0) { 2374 nfields = 128; 2375 } else { 2376 nfields = patch->nsubspaces; 2377 } 2378 ierr = PetscMalloc1(nfields, &ifields);CHKERRQ(ierr); 2379 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_exclude_subspaces", patch->classname); 2380 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->options,((PetscObject)pc)->prefix,option,ifields,&nfields,&flg);CHKERRQ(ierr); 2381 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"); 2382 if (flg) { 2383 PetscHSetIClear(patch->subspaces_to_exclude); 2384 for (k = 0; k < nfields; k++) { 2385 PetscHSetIAdd(patch->subspaces_to_exclude, ifields[k]); 2386 } 2387 } 2388 ierr = PetscFree(ifields);CHKERRQ(ierr); 2389 2390 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_patches_view", patch->classname); 2391 ierr = PetscOptionsBool(option, "Print out information during patch construction", "PCPATCH", patch->viewPatches, &patch->viewPatches, &flg);CHKERRQ(ierr); 2392 2393 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_cells_view", patch->classname); 2394 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, option, &patch->viewerCells, &patch->formatCells, &patch->viewCells);CHKERRQ(ierr); 2395 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_points_view", patch->classname); 2396 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, option, &patch->viewerPoints, &patch->formatPoints, &patch->viewPoints);CHKERRQ(ierr); 2397 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_section_view", patch->classname); 2398 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, option, &patch->viewerSection, &patch->formatSection, &patch->viewSection);CHKERRQ(ierr); 2399 ierr = PetscSNPrintf(option, PETSC_MAX_PATH_LEN, "-%s_patch_mat_view", patch->classname); 2400 ierr = PetscOptionsGetViewer(comm,((PetscObject) pc)->options,prefix, option, &patch->viewerMatrix, &patch->formatMatrix, &patch->viewMatrix);CHKERRQ(ierr); 2401 ierr = PetscOptionsTail();CHKERRQ(ierr); 2402 patch->optionsSet = PETSC_TRUE; 2403 PetscFunctionReturn(0); 2404 } 2405 2406 static PetscErrorCode PCSetUpOnBlocks_PATCH(PC pc) 2407 { 2408 PC_PATCH *patch = (PC_PATCH*) pc->data; 2409 KSPConvergedReason reason; 2410 PetscInt i; 2411 PetscErrorCode ierr; 2412 2413 PetscFunctionBegin; 2414 if (!patch->save_operators) { 2415 /* Can't do this here because the sub KSPs don't have an operator attached yet. */ 2416 PetscFunctionReturn(0); 2417 } 2418 for (i = 0; i < patch->npatch; ++i) { 2419 if (!((KSP) patch->solver[i])->setfromoptionscalled) { 2420 ierr = KSPSetFromOptions((KSP) patch->solver[i]);CHKERRQ(ierr); 2421 } 2422 ierr = KSPSetUp((KSP) patch->solver[i]);CHKERRQ(ierr); 2423 ierr = KSPGetConvergedReason((KSP) patch->solver[i], &reason);CHKERRQ(ierr); 2424 if (reason == KSP_DIVERGED_PCSETUP_FAILED) pc->failedreason = PC_SUBPC_ERROR; 2425 } 2426 PetscFunctionReturn(0); 2427 } 2428 2429 static PetscErrorCode PCView_PATCH(PC pc, PetscViewer viewer) 2430 { 2431 PC_PATCH *patch = (PC_PATCH *) pc->data; 2432 PetscViewer sviewer; 2433 PetscBool isascii; 2434 PetscMPIInt rank; 2435 PetscErrorCode ierr; 2436 2437 PetscFunctionBegin; 2438 /* TODO Redo tabbing with set tbas in new style */ 2439 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 2440 if (!isascii) PetscFunctionReturn(0); 2441 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) pc), &rank);CHKERRQ(ierr); 2442 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 2443 ierr = PetscViewerASCIIPrintf(viewer, "Subspace Correction preconditioner with %d patches\n", patch->npatch);CHKERRQ(ierr); 2444 if(patch->local_composition_type == PC_COMPOSITE_MULTIPLICATIVE) { 2445 ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: multiplicative\n");CHKERRQ(ierr); 2446 } else { 2447 ierr = PetscViewerASCIIPrintf(viewer, "Schwarz type: additive\n");CHKERRQ(ierr); 2448 } 2449 if (patch->partition_of_unity) {ierr = PetscViewerASCIIPrintf(viewer, "Weighting by partition of unity\n");CHKERRQ(ierr);} 2450 else {ierr = PetscViewerASCIIPrintf(viewer, "Not weighting by partition of unity\n");CHKERRQ(ierr);} 2451 if (patch->symmetrise_sweep) {ierr = PetscViewerASCIIPrintf(viewer, "Symmetrising sweep (start->end, then end->start)\n");CHKERRQ(ierr);} 2452 else {ierr = PetscViewerASCIIPrintf(viewer, "Not symmetrising sweep\n");CHKERRQ(ierr);} 2453 if (!patch->save_operators) {ierr = PetscViewerASCIIPrintf(viewer, "Not saving patch operators (rebuilt every PCApply)\n");CHKERRQ(ierr);} 2454 else {ierr = PetscViewerASCIIPrintf(viewer, "Saving patch operators (rebuilt every PCSetUp)\n");CHKERRQ(ierr);} 2455 if (patch->patchconstructop == PCPatchConstruct_Star) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: star\n");CHKERRQ(ierr);} 2456 else if (patch->patchconstructop == PCPatchConstruct_Vanka) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: Vanka\n");CHKERRQ(ierr);} 2457 else if (patch->patchconstructop == PCPatchConstruct_User) {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: user-specified\n");CHKERRQ(ierr);} 2458 else {ierr = PetscViewerASCIIPrintf(viewer, "Patch construction operator: unknown\n");CHKERRQ(ierr);} 2459 ierr = PetscViewerASCIIPrintf(viewer, "Solver on patches (all same):\n");CHKERRQ(ierr); 2460 if (patch->solver) { 2461 ierr = PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 2462 if (!rank) { 2463 ierr = PetscViewerASCIIPushTab(sviewer);CHKERRQ(ierr); 2464 ierr = PetscObjectView(patch->solver[0], sviewer);CHKERRQ(ierr); 2465 ierr = PetscViewerASCIIPopTab(sviewer);CHKERRQ(ierr); 2466 } 2467 ierr = PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);CHKERRQ(ierr); 2468 } else { 2469 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 2470 ierr = PetscViewerASCIIPrintf(viewer, "Solver not yet set.\n");CHKERRQ(ierr); 2471 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 2472 } 2473 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 2474 PetscFunctionReturn(0); 2475 } 2476 2477 /*MC 2478 PCPATCH - A PC object that encapsulates flexible definition of blocks for overlapping and non-overlapping 2479 small block additive preconditioners. Block definition is based on topology from 2480 a DM and equation numbering from a PetscSection. 2481 2482 Options Database Keys: 2483 + -pc_patch_cells_view - Views the process local cell numbers for each patch 2484 . -pc_patch_points_view - Views the process local mesh point numbers for each patch 2485 . -pc_patch_g2l_view - Views the map between global dofs and patch local dofs for each patch 2486 . -pc_patch_patches_view - Views the global dofs associated with each patch and its boundary 2487 - -pc_patch_sub_mat_view - Views the matrix associated with each patch 2488 2489 Level: intermediate 2490 2491 .seealso: PCType, PCCreate(), PCSetType() 2492 M*/ 2493 PETSC_EXTERN PetscErrorCode PCCreate_Patch(PC pc) 2494 { 2495 PC_PATCH *patch; 2496 PetscErrorCode ierr; 2497 2498 PetscFunctionBegin; 2499 ierr = PetscNewLog(pc, &patch);CHKERRQ(ierr); 2500 2501 if (patch->subspaces_to_exclude) { 2502 PetscHSetIDestroy(&patch->subspaces_to_exclude); 2503 } 2504 PetscHSetICreate(&patch->subspaces_to_exclude); 2505 2506 patch->classname = "pc"; 2507 patch->isNonlinear = PETSC_FALSE; 2508 2509 /* Set some defaults */ 2510 patch->combined = PETSC_FALSE; 2511 patch->save_operators = PETSC_TRUE; 2512 patch->local_composition_type = PC_COMPOSITE_ADDITIVE; 2513 patch->partition_of_unity = PETSC_FALSE; 2514 patch->codim = -1; 2515 patch->dim = -1; 2516 patch->vankadim = -1; 2517 patch->ignoredim = -1; 2518 patch->patchconstructop = PCPatchConstruct_Star; 2519 patch->symmetrise_sweep = PETSC_FALSE; 2520 patch->npatch = 0; 2521 patch->userIS = NULL; 2522 patch->optionsSet = PETSC_FALSE; 2523 patch->iterationSet = NULL; 2524 patch->user_patches = PETSC_FALSE; 2525 ierr = PetscStrallocpy(MATDENSE, (char **) &patch->sub_mat_type);CHKERRQ(ierr); 2526 patch->viewPatches = PETSC_FALSE; 2527 patch->viewCells = PETSC_FALSE; 2528 patch->viewPoints = PETSC_FALSE; 2529 patch->viewSection = PETSC_FALSE; 2530 patch->viewMatrix = PETSC_FALSE; 2531 patch->setupsolver = PCSetUp_PATCH_Linear; 2532 patch->applysolver = PCApply_PATCH_Linear; 2533 patch->resetsolver = PCReset_PATCH_Linear; 2534 patch->destroysolver = PCDestroy_PATCH_Linear; 2535 patch->updatemultiplicative = PCUpdateMultiplicative_PATCH_Linear; 2536 2537 pc->data = (void *) patch; 2538 pc->ops->apply = PCApply_PATCH; 2539 pc->ops->applytranspose = 0; /* PCApplyTranspose_PATCH; */ 2540 pc->ops->setup = PCSetUp_PATCH; 2541 pc->ops->reset = PCReset_PATCH; 2542 pc->ops->destroy = PCDestroy_PATCH; 2543 pc->ops->setfromoptions = PCSetFromOptions_PATCH; 2544 pc->ops->setuponblocks = PCSetUpOnBlocks_PATCH; 2545 pc->ops->view = PCView_PATCH; 2546 pc->ops->applyrichardson = 0; 2547 2548 PetscFunctionReturn(0); 2549 } 2550