1 /* 2 This file defines an "generalized" additive Schwarz preconditioner for any Mat implementation. 3 In this version each MPI rank may intersect multiple subdomains and any subdomain may 4 intersect multiple MPI ranks. Intersections of subdomains with MPI ranks are called *local 5 subdomains*. 6 7 N - total number of distinct global subdomains (set explicitly in PCGASMSetTotalSubdomains() or implicitly PCGASMSetSubdomains() and then calculated in PCSetUp_GASM()) 8 n - actual number of local subdomains on this rank (set in PCGASMSetSubdomains() or calculated in PCGASMSetTotalSubdomains()) 9 nmax - maximum number of local subdomains per rank (calculated in PCSetUp_GASM()) 10 */ 11 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 12 #include <petscdm.h> 13 14 typedef struct { 15 PetscInt N, n, nmax; 16 PetscInt overlap; /* overlap requested by user */ 17 PCGASMType type; /* use reduced interpolation, restriction or both */ 18 PetscBool type_set; /* if user set this value (so won't change it for symmetric problems) */ 19 PetscBool same_subdomain_solvers; /* flag indicating whether all local solvers are same */ 20 PetscBool sort_indices; /* flag to sort subdomain indices */ 21 PetscBool user_subdomains; /* whether the user set explicit subdomain index sets -- keep them on PCReset() */ 22 PetscBool dm_subdomains; /* whether DM is allowed to define subdomains */ 23 PetscBool hierarchicalpartitioning; 24 IS *ois; /* index sets that define the outer (conceptually, overlapping) subdomains */ 25 IS *iis; /* index sets that define the inner (conceptually, nonoverlapping) subdomains */ 26 KSP *ksp; /* linear solvers for each subdomain */ 27 Mat *pmat; /* subdomain block matrices */ 28 Vec gx, gy; /* Merged work vectors */ 29 Vec *x, *y; /* Split work vectors; storage aliases pieces of storage of the above merged vectors. */ 30 VecScatter gorestriction; /* merged restriction to disjoint union of outer subdomains */ 31 VecScatter girestriction; /* merged restriction to disjoint union of inner subdomains */ 32 VecScatter pctoouter; 33 IS permutationIS; 34 Mat permutationP; 35 Mat pcmat; 36 Vec pcx, pcy; 37 } PC_GASM; 38 39 static PetscErrorCode PCGASMComputeGlobalSubdomainNumbering_Private(PC pc, PetscInt **numbering, PetscInt **permutation) 40 { 41 PC_GASM *osm = (PC_GASM *)pc->data; 42 PetscInt i; 43 44 PetscFunctionBegin; 45 /* Determine the number of globally-distinct subdomains and compute a global numbering for them. */ 46 PetscCall(PetscMalloc2(osm->n, numbering, osm->n, permutation)); 47 PetscCall(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->iis, NULL, *numbering)); 48 for (i = 0; i < osm->n; ++i) (*permutation)[i] = i; 49 PetscCall(PetscSortIntWithPermutation(osm->n, *numbering, *permutation)); 50 PetscFunctionReturn(PETSC_SUCCESS); 51 } 52 53 static PetscErrorCode PCGASMSubdomainView_Private(PC pc, PetscInt i, PetscViewer viewer) 54 { 55 PC_GASM *osm = (PC_GASM *)pc->data; 56 PetscInt j, nidx; 57 const PetscInt *idx; 58 PetscViewer sviewer; 59 char *cidx; 60 61 PetscFunctionBegin; 62 PetscCheck(i >= -1 && i < osm->n, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Invalid subdomain %" PetscInt_FMT ": must nonnegative and less than %" PetscInt_FMT, i, osm->n); 63 64 /* Inner subdomains. */ 65 /* 66 No more than 15 characters per index plus a space. 67 PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx, 68 in case nidx == 0. That will take care of the space for the trailing '\0' as well. 69 For nidx == 0, the whole string 16 '\0'. 70 */ 71 PetscCall(PetscViewerASCIIPrintf(viewer, "Inner subdomain:\n")); 72 PetscCall(PetscViewerFlush(viewer)); 73 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 74 if (i > -1) { 75 PetscCall(ISGetLocalSize(osm->iis[i], &nidx)); 76 PetscCall(PetscMalloc1(16 * (nidx + 1) + 1, &cidx)); 77 PetscCall(PetscViewerStringOpen(PETSC_COMM_SELF, cidx, 16 * (nidx + 1) + 1, &sviewer)); 78 PetscCall(ISGetIndices(osm->iis[i], &idx)); 79 for (j = 0; j < nidx; ++j) PetscCall(PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j])); 80 PetscCall(ISRestoreIndices(osm->iis[i], &idx)); 81 PetscCall(PetscViewerDestroy(&sviewer)); 82 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx)); 83 PetscCall(PetscFree(cidx)); 84 } 85 PetscCall(PetscViewerFlush(viewer)); 86 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 87 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 88 PetscCall(PetscViewerFlush(viewer)); 89 90 /* Outer subdomains. */ 91 /* 92 No more than 15 characters per index plus a space. 93 PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx, 94 in case nidx == 0. That will take care of the space for the trailing '\0' as well. 95 For nidx == 0, the whole string 16 '\0'. 96 */ 97 PetscCall(PetscViewerASCIIPrintf(viewer, "Outer subdomain:\n")); 98 PetscCall(PetscViewerFlush(viewer)); 99 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 100 if (i > -1) { 101 PetscCall(ISGetLocalSize(osm->ois[i], &nidx)); 102 PetscCall(PetscMalloc1(16 * (nidx + 1) + 1, &cidx)); 103 PetscCall(PetscViewerStringOpen(PETSC_COMM_SELF, cidx, 16 * (nidx + 1) + 1, &sviewer)); 104 PetscCall(ISGetIndices(osm->ois[i], &idx)); 105 for (j = 0; j < nidx; ++j) PetscCall(PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j])); 106 PetscCall(PetscViewerDestroy(&sviewer)); 107 PetscCall(ISRestoreIndices(osm->ois[i], &idx)); 108 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx)); 109 PetscCall(PetscFree(cidx)); 110 } 111 PetscCall(PetscViewerFlush(viewer)); 112 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 113 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 114 PetscCall(PetscViewerFlush(viewer)); 115 PetscFunctionReturn(PETSC_SUCCESS); 116 } 117 118 static PetscErrorCode PCGASMPrintSubdomains(PC pc) 119 { 120 PC_GASM *osm = (PC_GASM *)pc->data; 121 const char *prefix; 122 char fname[PETSC_MAX_PATH_LEN + 1]; 123 PetscInt l, d, count; 124 PetscBool found; 125 PetscViewer viewer; 126 PetscInt *numbering, *permutation; /* global numbering of locally-supported subdomains and the permutation from the local ordering */ 127 128 PetscFunctionBegin; 129 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 130 PetscCall(PetscOptionsHasName(NULL, prefix, "-pc_gasm_print_subdomains", &found)); 131 if (!found) PetscFunctionReturn(PETSC_SUCCESS); 132 PetscCall(PetscOptionsGetString(NULL, prefix, "-pc_gasm_print_subdomains", fname, sizeof(fname), &found)); 133 if (!found) PetscCall(PetscStrncpy(fname, "stdout", sizeof(fname))); 134 PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc), fname, &viewer)); 135 /* 136 Make sure the viewer has a name. Otherwise this may cause a deadlock or other weird errors when creating a subcomm viewer: 137 the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed collectively on the comm. 138 */ 139 PetscCall(PetscObjectName((PetscObject)viewer)); 140 l = 0; 141 PetscCall(PCGASMComputeGlobalSubdomainNumbering_Private(pc, &numbering, &permutation)); 142 for (count = 0; count < osm->N; ++count) { 143 /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */ 144 if (l < osm->n) { 145 d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */ 146 if (numbering[d] == count) l++; 147 else d = -1; 148 } else d = -1; 149 PetscCall(PCGASMSubdomainView_Private(pc, d, viewer)); 150 } 151 PetscCall(PetscFree2(numbering, permutation)); 152 PetscCall(PetscViewerDestroy(&viewer)); 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155 156 static PetscErrorCode PCView_GASM(PC pc, PetscViewer viewer) 157 { 158 PC_GASM *osm = (PC_GASM *)pc->data; 159 const char *prefix; 160 PetscMPIInt rank, size; 161 PetscInt bsz; 162 PetscBool iascii, view_subdomains = PETSC_FALSE; 163 PetscViewer sviewer; 164 PetscInt count, l; 165 char overlap[256] = "user-defined overlap"; 166 char gsubdomains[256] = "unknown total number of subdomains"; 167 char msubdomains[256] = "unknown max number of local subdomains"; 168 PetscInt *numbering, *permutation; /* global numbering of locally-supported subdomains and the permutation from the local ordering */ 169 170 PetscFunctionBegin; 171 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 172 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 173 174 if (osm->overlap >= 0) PetscCall(PetscSNPrintf(overlap, sizeof(overlap), "requested amount of overlap = %" PetscInt_FMT, osm->overlap)); 175 if (osm->N != PETSC_DETERMINE) PetscCall(PetscSNPrintf(gsubdomains, sizeof(gsubdomains), "total number of subdomains = %" PetscInt_FMT, osm->N)); 176 if (osm->nmax != PETSC_DETERMINE) PetscCall(PetscSNPrintf(msubdomains, sizeof(msubdomains), "max number of local subdomains = %" PetscInt_FMT, osm->nmax)); 177 178 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 179 PetscCall(PetscOptionsGetBool(NULL, prefix, "-pc_gasm_view_subdomains", &view_subdomains, NULL)); 180 181 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 182 if (iascii) { 183 /* 184 Make sure the viewer has a name. Otherwise this may cause a deadlock when creating a subcomm viewer: 185 the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed 186 collectively on the comm. 187 */ 188 PetscCall(PetscObjectName((PetscObject)viewer)); 189 PetscCall(PetscViewerASCIIPrintf(viewer, " Restriction/interpolation type: %s\n", PCGASMTypes[osm->type])); 190 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", overlap)); 191 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", gsubdomains)); 192 PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", msubdomains)); 193 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 194 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d|%d] number of locally-supported subdomains = %" PetscInt_FMT "\n", rank, size, osm->n)); 195 PetscCall(PetscViewerFlush(viewer)); 196 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 197 /* Cannot take advantage of osm->same_subdomain_solvers without a global numbering of subdomains. */ 198 PetscCall(PetscViewerASCIIPrintf(viewer, " Subdomain solver info is as follows:\n")); 199 PetscCall(PetscViewerASCIIPushTab(viewer)); 200 PetscCall(PetscViewerASCIIPrintf(viewer, " - - - - - - - - - - - - - - - - - -\n")); 201 /* Make sure that everybody waits for the banner to be printed. */ 202 PetscCallMPI(MPI_Barrier(PetscObjectComm((PetscObject)viewer))); 203 /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */ 204 PetscCall(PCGASMComputeGlobalSubdomainNumbering_Private(pc, &numbering, &permutation)); 205 l = 0; 206 for (count = 0; count < osm->N; ++count) { 207 PetscMPIInt srank, ssize; 208 if (l < osm->n) { 209 PetscInt d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */ 210 if (numbering[d] == count) { 211 PetscCallMPI(MPI_Comm_size(((PetscObject)osm->ois[d])->comm, &ssize)); 212 PetscCallMPI(MPI_Comm_rank(((PetscObject)osm->ois[d])->comm, &srank)); 213 PetscCall(PetscViewerGetSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer)); 214 PetscCall(ISGetLocalSize(osm->ois[d], &bsz)); 215 PetscCall(PetscViewerASCIISynchronizedPrintf(sviewer, " [%d|%d] (subcomm [%d|%d]) local subdomain number %" PetscInt_FMT ", local size = %" PetscInt_FMT "\n", rank, size, srank, ssize, d, bsz)); 216 PetscCall(PetscViewerFlush(sviewer)); 217 if (view_subdomains) PetscCall(PCGASMSubdomainView_Private(pc, d, sviewer)); 218 if (!pc->setupcalled) { 219 PetscCall(PetscViewerASCIIPrintf(sviewer, " Solver not set up yet: PCSetUp() not yet called\n")); 220 } else { 221 PetscCall(KSPView(osm->ksp[d], sviewer)); 222 } 223 PetscCall(PetscViewerASCIIPrintf(sviewer, " - - - - - - - - - - - - - - - - - -\n")); 224 PetscCall(PetscViewerFlush(sviewer)); 225 PetscCall(PetscViewerRestoreSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer)); 226 ++l; 227 } 228 } 229 PetscCallMPI(MPI_Barrier(PetscObjectComm((PetscObject)pc))); 230 } 231 PetscCall(PetscFree2(numbering, permutation)); 232 PetscCall(PetscViewerASCIIPopTab(viewer)); 233 PetscCall(PetscViewerFlush(viewer)); 234 /* this line is needed to match the extra PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */ 235 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 236 } 237 PetscFunctionReturn(PETSC_SUCCESS); 238 } 239 240 PETSC_INTERN PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[]); 241 242 static PetscErrorCode PCGASMSetHierarchicalPartitioning(PC pc) 243 { 244 PC_GASM *osm = (PC_GASM *)pc->data; 245 MatPartitioning part; 246 MPI_Comm comm; 247 PetscMPIInt size; 248 PetscInt nlocalsubdomains, fromrows_localsize; 249 IS partitioning, fromrows, isn; 250 Vec outervec; 251 252 PetscFunctionBegin; 253 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 254 PetscCallMPI(MPI_Comm_size(comm, &size)); 255 /* we do not need a hierarchical partitioning when 256 * the total number of subdomains is consistent with 257 * the number of MPI tasks. 258 * For the following cases, we do not need to use HP 259 * */ 260 if (osm->N == PETSC_DETERMINE || osm->N >= size || osm->N == 1) PetscFunctionReturn(PETSC_SUCCESS); 261 PetscCheck(size % osm->N == 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "have to specify the total number of subdomains %" PetscInt_FMT " to be a factor of the number of ranks %d ", osm->N, size); 262 nlocalsubdomains = size / osm->N; 263 osm->n = 1; 264 PetscCall(MatPartitioningCreate(comm, &part)); 265 PetscCall(MatPartitioningSetAdjacency(part, pc->pmat)); 266 PetscCall(MatPartitioningSetType(part, MATPARTITIONINGHIERARCH)); 267 PetscCall(MatPartitioningHierarchicalSetNcoarseparts(part, osm->N)); 268 PetscCall(MatPartitioningHierarchicalSetNfineparts(part, nlocalsubdomains)); 269 PetscCall(MatPartitioningSetFromOptions(part)); 270 /* get new rank owner number of each vertex */ 271 PetscCall(MatPartitioningApply(part, &partitioning)); 272 PetscCall(ISBuildTwoSided(partitioning, NULL, &fromrows)); 273 PetscCall(ISPartitioningToNumbering(partitioning, &isn)); 274 PetscCall(ISDestroy(&isn)); 275 PetscCall(ISGetLocalSize(fromrows, &fromrows_localsize)); 276 PetscCall(MatPartitioningDestroy(&part)); 277 PetscCall(MatCreateVecs(pc->pmat, &outervec, NULL)); 278 PetscCall(VecCreateMPI(comm, fromrows_localsize, PETSC_DETERMINE, &(osm->pcx))); 279 PetscCall(VecDuplicate(osm->pcx, &(osm->pcy))); 280 PetscCall(VecScatterCreate(osm->pcx, NULL, outervec, fromrows, &(osm->pctoouter))); 281 PetscCall(MatCreateSubMatrix(pc->pmat, fromrows, fromrows, MAT_INITIAL_MATRIX, &(osm->permutationP))); 282 PetscCall(PetscObjectReference((PetscObject)fromrows)); 283 osm->permutationIS = fromrows; 284 osm->pcmat = pc->pmat; 285 PetscCall(PetscObjectReference((PetscObject)osm->permutationP)); 286 pc->pmat = osm->permutationP; 287 PetscCall(VecDestroy(&outervec)); 288 PetscCall(ISDestroy(&fromrows)); 289 PetscCall(ISDestroy(&partitioning)); 290 osm->n = PETSC_DETERMINE; 291 PetscFunctionReturn(PETSC_SUCCESS); 292 } 293 294 static PetscErrorCode PCSetUp_GASM(PC pc) 295 { 296 PC_GASM *osm = (PC_GASM *)pc->data; 297 PetscInt i, nInnerIndices, nTotalInnerIndices; 298 PetscMPIInt rank, size; 299 MatReuse scall = MAT_REUSE_MATRIX; 300 KSP ksp; 301 PC subpc; 302 const char *prefix, *pprefix; 303 Vec x, y; 304 PetscInt oni; /* Number of indices in the i-th local outer subdomain. */ 305 const PetscInt *oidxi; /* Indices from the i-th subdomain local outer subdomain. */ 306 PetscInt on; /* Number of indices in the disjoint union of local outer subdomains. */ 307 PetscInt *oidx; /* Indices in the disjoint union of local outer subdomains. */ 308 IS gois; /* Disjoint union the global indices of outer subdomains. */ 309 IS goid; /* Identity IS of the size of the disjoint union of outer subdomains. */ 310 PetscScalar *gxarray, *gyarray; 311 PetscInt gostart; /* Start of locally-owned indices in the vectors -- osm->gx,osm->gy -- over the disjoint union of outer subdomains. */ 312 PetscInt num_subdomains = 0; 313 DM *subdomain_dm = NULL; 314 char **subdomain_names = NULL; 315 PetscInt *numbering; 316 317 PetscFunctionBegin; 318 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 319 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 320 if (!pc->setupcalled) { 321 /* use a hierarchical partitioning */ 322 if (osm->hierarchicalpartitioning) PetscCall(PCGASMSetHierarchicalPartitioning(pc)); 323 if (osm->n == PETSC_DETERMINE) { 324 if (osm->N != PETSC_DETERMINE) { 325 /* No local subdomains given, but the desired number of total subdomains is known, so construct them accordingly. */ 326 PetscCall(PCGASMCreateSubdomains(pc->pmat, osm->N, &osm->n, &osm->iis)); 327 } else if (osm->dm_subdomains && pc->dm) { 328 /* try pc->dm next, if allowed */ 329 PetscInt d; 330 IS *inner_subdomain_is, *outer_subdomain_is; 331 PetscCall(DMCreateDomainDecomposition(pc->dm, &num_subdomains, &subdomain_names, &inner_subdomain_is, &outer_subdomain_is, &subdomain_dm)); 332 if (num_subdomains) PetscCall(PCGASMSetSubdomains(pc, num_subdomains, inner_subdomain_is, outer_subdomain_is)); 333 for (d = 0; d < num_subdomains; ++d) { 334 if (inner_subdomain_is) PetscCall(ISDestroy(&inner_subdomain_is[d])); 335 if (outer_subdomain_is) PetscCall(ISDestroy(&outer_subdomain_is[d])); 336 } 337 PetscCall(PetscFree(inner_subdomain_is)); 338 PetscCall(PetscFree(outer_subdomain_is)); 339 } else { 340 /* still no subdomains; use one per rank */ 341 osm->nmax = osm->n = 1; 342 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 343 osm->N = size; 344 PetscCall(PCGASMCreateLocalSubdomains(pc->pmat, osm->n, &osm->iis)); 345 } 346 } 347 if (!osm->iis) { 348 /* 349 osm->n was set in PCGASMSetSubdomains(), but the actual subdomains have not been supplied. 350 We create the requisite number of local inner subdomains and then expand them into 351 out subdomains, if necessary. 352 */ 353 PetscCall(PCGASMCreateLocalSubdomains(pc->pmat, osm->n, &osm->iis)); 354 } 355 if (!osm->ois) { 356 /* 357 Initially make outer subdomains the same as inner subdomains. If nonzero additional overlap 358 has been requested, copy the inner subdomains over so they can be modified. 359 */ 360 PetscCall(PetscMalloc1(osm->n, &osm->ois)); 361 for (i = 0; i < osm->n; ++i) { 362 if (osm->overlap > 0 && osm->N > 1) { /* With positive overlap, osm->iis[i] will be modified */ 363 PetscCall(ISDuplicate(osm->iis[i], (osm->ois) + i)); 364 PetscCall(ISCopy(osm->iis[i], osm->ois[i])); 365 } else { 366 PetscCall(PetscObjectReference((PetscObject)((osm->iis)[i]))); 367 osm->ois[i] = osm->iis[i]; 368 } 369 } 370 if (osm->overlap > 0 && osm->N > 1) { 371 /* Extend the "overlapping" regions by a number of steps */ 372 PetscCall(MatIncreaseOverlapSplit(pc->pmat, osm->n, osm->ois, osm->overlap)); 373 } 374 } 375 376 /* Now the subdomains are defined. Determine their global and max local numbers, if necessary. */ 377 if (osm->nmax == PETSC_DETERMINE) { 378 PetscMPIInt inwork, outwork; 379 /* determine global number of subdomains and the max number of local subdomains */ 380 inwork = osm->n; 381 PetscCall(MPIU_Allreduce(&inwork, &outwork, 1, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)pc))); 382 osm->nmax = outwork; 383 } 384 if (osm->N == PETSC_DETERMINE) { 385 /* Determine the number of globally-distinct subdomains and compute a global numbering for them. */ 386 PetscCall(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->ois, &osm->N, NULL)); 387 } 388 389 if (osm->sort_indices) { 390 for (i = 0; i < osm->n; i++) { 391 PetscCall(ISSort(osm->ois[i])); 392 PetscCall(ISSort(osm->iis[i])); 393 } 394 } 395 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 396 PetscCall(PCGASMPrintSubdomains(pc)); 397 398 /* 399 Merge the ISs, create merged vectors and restrictions. 400 */ 401 /* Merge outer subdomain ISs and construct a restriction onto the disjoint union of local outer subdomains. */ 402 on = 0; 403 for (i = 0; i < osm->n; i++) { 404 PetscCall(ISGetLocalSize(osm->ois[i], &oni)); 405 on += oni; 406 } 407 PetscCall(PetscMalloc1(on, &oidx)); 408 on = 0; 409 /* Merge local indices together */ 410 for (i = 0; i < osm->n; i++) { 411 PetscCall(ISGetLocalSize(osm->ois[i], &oni)); 412 PetscCall(ISGetIndices(osm->ois[i], &oidxi)); 413 PetscCall(PetscArraycpy(oidx + on, oidxi, oni)); 414 PetscCall(ISRestoreIndices(osm->ois[i], &oidxi)); 415 on += oni; 416 } 417 PetscCall(ISCreateGeneral(((PetscObject)(pc))->comm, on, oidx, PETSC_OWN_POINTER, &gois)); 418 nTotalInnerIndices = 0; 419 for (i = 0; i < osm->n; i++) { 420 PetscCall(ISGetLocalSize(osm->iis[i], &nInnerIndices)); 421 nTotalInnerIndices += nInnerIndices; 422 } 423 PetscCall(VecCreateMPI(((PetscObject)(pc))->comm, nTotalInnerIndices, PETSC_DETERMINE, &x)); 424 PetscCall(VecDuplicate(x, &y)); 425 426 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pc), on, PETSC_DECIDE, &osm->gx)); 427 PetscCall(VecDuplicate(osm->gx, &osm->gy)); 428 PetscCall(VecGetOwnershipRange(osm->gx, &gostart, NULL)); 429 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), on, gostart, 1, &goid)); 430 /* gois might indices not on local */ 431 PetscCall(VecScatterCreate(x, gois, osm->gx, goid, &(osm->gorestriction))); 432 PetscCall(PetscMalloc1(osm->n, &numbering)); 433 PetscCall(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->ois, NULL, numbering)); 434 PetscCall(VecDestroy(&x)); 435 PetscCall(ISDestroy(&gois)); 436 437 /* Merge inner subdomain ISs and construct a restriction onto the disjoint union of local inner subdomains. */ 438 { 439 PetscInt ini; /* Number of indices the i-th a local inner subdomain. */ 440 PetscInt in; /* Number of indices in the disjoint union of local inner subdomains. */ 441 PetscInt *iidx; /* Global indices in the merged local inner subdomain. */ 442 PetscInt *ioidx; /* Global indices of the disjoint union of inner subdomains within the disjoint union of outer subdomains. */ 443 IS giis; /* IS for the disjoint union of inner subdomains. */ 444 IS giois; /* IS for the disjoint union of inner subdomains within the disjoint union of outer subdomains. */ 445 PetscScalar *array; 446 const PetscInt *indices; 447 PetscInt k; 448 on = 0; 449 for (i = 0; i < osm->n; i++) { 450 PetscCall(ISGetLocalSize(osm->ois[i], &oni)); 451 on += oni; 452 } 453 PetscCall(PetscMalloc1(on, &iidx)); 454 PetscCall(PetscMalloc1(on, &ioidx)); 455 PetscCall(VecGetArray(y, &array)); 456 /* set communicator id to determine where overlap is */ 457 in = 0; 458 for (i = 0; i < osm->n; i++) { 459 PetscCall(ISGetLocalSize(osm->iis[i], &ini)); 460 for (k = 0; k < ini; ++k) array[in + k] = numbering[i]; 461 in += ini; 462 } 463 PetscCall(VecRestoreArray(y, &array)); 464 PetscCall(VecScatterBegin(osm->gorestriction, y, osm->gy, INSERT_VALUES, SCATTER_FORWARD)); 465 PetscCall(VecScatterEnd(osm->gorestriction, y, osm->gy, INSERT_VALUES, SCATTER_FORWARD)); 466 PetscCall(VecGetOwnershipRange(osm->gy, &gostart, NULL)); 467 PetscCall(VecGetArray(osm->gy, &array)); 468 on = 0; 469 in = 0; 470 for (i = 0; i < osm->n; i++) { 471 PetscCall(ISGetLocalSize(osm->ois[i], &oni)); 472 PetscCall(ISGetIndices(osm->ois[i], &indices)); 473 for (k = 0; k < oni; k++) { 474 /* skip overlapping indices to get inner domain */ 475 if (PetscRealPart(array[on + k]) != numbering[i]) continue; 476 iidx[in] = indices[k]; 477 ioidx[in++] = gostart + on + k; 478 } 479 PetscCall(ISRestoreIndices(osm->ois[i], &indices)); 480 on += oni; 481 } 482 PetscCall(VecRestoreArray(osm->gy, &array)); 483 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), in, iidx, PETSC_OWN_POINTER, &giis)); 484 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), in, ioidx, PETSC_OWN_POINTER, &giois)); 485 PetscCall(VecScatterCreate(y, giis, osm->gy, giois, &osm->girestriction)); 486 PetscCall(VecDestroy(&y)); 487 PetscCall(ISDestroy(&giis)); 488 PetscCall(ISDestroy(&giois)); 489 } 490 PetscCall(ISDestroy(&goid)); 491 PetscCall(PetscFree(numbering)); 492 493 /* Create the subdomain work vectors. */ 494 PetscCall(PetscMalloc1(osm->n, &osm->x)); 495 PetscCall(PetscMalloc1(osm->n, &osm->y)); 496 PetscCall(VecGetArray(osm->gx, &gxarray)); 497 PetscCall(VecGetArray(osm->gy, &gyarray)); 498 for (i = 0, on = 0; i < osm->n; ++i, on += oni) { 499 PetscInt oNi; 500 PetscCall(ISGetLocalSize(osm->ois[i], &oni)); 501 /* on a sub communicator */ 502 PetscCall(ISGetSize(osm->ois[i], &oNi)); 503 PetscCall(VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm, 1, oni, oNi, gxarray + on, &osm->x[i])); 504 PetscCall(VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm, 1, oni, oNi, gyarray + on, &osm->y[i])); 505 } 506 PetscCall(VecRestoreArray(osm->gx, &gxarray)); 507 PetscCall(VecRestoreArray(osm->gy, &gyarray)); 508 /* Create the subdomain solvers */ 509 PetscCall(PetscMalloc1(osm->n, &osm->ksp)); 510 for (i = 0; i < osm->n; i++) { 511 char subprefix[PETSC_MAX_PATH_LEN + 1]; 512 PetscCall(KSPCreate(((PetscObject)(osm->ois[i]))->comm, &ksp)); 513 PetscCall(KSPSetNestLevel(ksp, pc->kspnestlevel)); 514 PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure)); 515 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1)); 516 PetscCall(KSPSetType(ksp, KSPPREONLY)); 517 PetscCall(KSPGetPC(ksp, &subpc)); /* Why do we need this here? */ 518 if (subdomain_dm) { 519 PetscCall(KSPSetDM(ksp, subdomain_dm[i])); 520 PetscCall(DMDestroy(subdomain_dm + i)); 521 } 522 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 523 PetscCall(KSPSetOptionsPrefix(ksp, prefix)); 524 if (subdomain_names && subdomain_names[i]) { 525 PetscCall(PetscSNPrintf(subprefix, PETSC_MAX_PATH_LEN, "sub_%s_", subdomain_names[i])); 526 PetscCall(KSPAppendOptionsPrefix(ksp, subprefix)); 527 PetscCall(PetscFree(subdomain_names[i])); 528 } 529 PetscCall(KSPAppendOptionsPrefix(ksp, "sub_")); 530 osm->ksp[i] = ksp; 531 } 532 PetscCall(PetscFree(subdomain_dm)); 533 PetscCall(PetscFree(subdomain_names)); 534 scall = MAT_INITIAL_MATRIX; 535 } else { /* if (pc->setupcalled) */ 536 /* 537 Destroy the submatrices from the previous iteration 538 */ 539 if (pc->flag == DIFFERENT_NONZERO_PATTERN) { 540 PetscCall(MatDestroyMatrices(osm->n, &osm->pmat)); 541 scall = MAT_INITIAL_MATRIX; 542 } 543 if (osm->permutationIS) { 544 PetscCall(MatCreateSubMatrix(pc->pmat, osm->permutationIS, osm->permutationIS, scall, &osm->permutationP)); 545 PetscCall(PetscObjectReference((PetscObject)osm->permutationP)); 546 osm->pcmat = pc->pmat; 547 pc->pmat = osm->permutationP; 548 } 549 } 550 551 /* 552 Extract the submatrices. 553 */ 554 if (size > 1) { 555 PetscCall(MatCreateSubMatricesMPI(pc->pmat, osm->n, osm->ois, osm->ois, scall, &osm->pmat)); 556 } else { 557 PetscCall(MatCreateSubMatrices(pc->pmat, osm->n, osm->ois, osm->ois, scall, &osm->pmat)); 558 } 559 if (scall == MAT_INITIAL_MATRIX) { 560 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc->pmat, &pprefix)); 561 for (i = 0; i < osm->n; i++) PetscCall(PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i], pprefix)); 562 } 563 564 /* Return control to the user so that the submatrices can be modified (e.g., to apply 565 different boundary conditions for the submatrices than for the global problem) */ 566 PetscCall(PCModifySubMatrices(pc, osm->n, osm->ois, osm->ois, osm->pmat, pc->modifysubmatricesP)); 567 568 /* 569 Loop over submatrices putting them into local ksps 570 */ 571 for (i = 0; i < osm->n; i++) { 572 PetscCall(KSPSetOperators(osm->ksp[i], osm->pmat[i], osm->pmat[i])); 573 PetscCall(KSPGetOptionsPrefix(osm->ksp[i], &prefix)); 574 PetscCall(MatSetOptionsPrefix(osm->pmat[i], prefix)); 575 if (!pc->setupcalled) PetscCall(KSPSetFromOptions(osm->ksp[i])); 576 } 577 if (osm->pcmat) { 578 PetscCall(MatDestroy(&pc->pmat)); 579 pc->pmat = osm->pcmat; 580 osm->pcmat = NULL; 581 } 582 PetscFunctionReturn(PETSC_SUCCESS); 583 } 584 585 static PetscErrorCode PCSetUpOnBlocks_GASM(PC pc) 586 { 587 PC_GASM *osm = (PC_GASM *)pc->data; 588 PetscInt i; 589 590 PetscFunctionBegin; 591 for (i = 0; i < osm->n; i++) PetscCall(KSPSetUp(osm->ksp[i])); 592 PetscFunctionReturn(PETSC_SUCCESS); 593 } 594 595 static PetscErrorCode PCApply_GASM(PC pc, Vec xin, Vec yout) 596 { 597 PC_GASM *osm = (PC_GASM *)pc->data; 598 PetscInt i; 599 Vec x, y; 600 ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE; 601 602 PetscFunctionBegin; 603 if (osm->pctoouter) { 604 PetscCall(VecScatterBegin(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 605 PetscCall(VecScatterEnd(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 606 x = osm->pcx; 607 y = osm->pcy; 608 } else { 609 x = xin; 610 y = yout; 611 } 612 /* 613 support for limiting the restriction or interpolation only to the inner 614 subdomain values (leaving the other values 0). 615 */ 616 if (!(osm->type & PC_GASM_RESTRICT)) { 617 /* have to zero the work RHS since scatter may leave some slots empty */ 618 PetscCall(VecZeroEntries(osm->gx)); 619 PetscCall(VecScatterBegin(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 620 } else { 621 PetscCall(VecScatterBegin(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 622 } 623 PetscCall(VecZeroEntries(osm->gy)); 624 if (!(osm->type & PC_GASM_RESTRICT)) { 625 PetscCall(VecScatterEnd(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 626 } else { 627 PetscCall(VecScatterEnd(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 628 } 629 /* do the subdomain solves */ 630 for (i = 0; i < osm->n; ++i) { 631 PetscCall(KSPSolve(osm->ksp[i], osm->x[i], osm->y[i])); 632 PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i])); 633 } 634 /* do we need to zero y? */ 635 PetscCall(VecZeroEntries(y)); 636 if (!(osm->type & PC_GASM_INTERPOLATE)) { 637 PetscCall(VecScatterBegin(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 638 PetscCall(VecScatterEnd(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 639 } else { 640 PetscCall(VecScatterBegin(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 641 PetscCall(VecScatterEnd(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 642 } 643 if (osm->pctoouter) { 644 PetscCall(VecScatterBegin(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 645 PetscCall(VecScatterEnd(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 646 } 647 PetscFunctionReturn(PETSC_SUCCESS); 648 } 649 650 static PetscErrorCode PCMatApply_GASM(PC pc, Mat Xin, Mat Yout) 651 { 652 PC_GASM *osm = (PC_GASM *)pc->data; 653 Mat X, Y, O = NULL, Z, W; 654 Vec x, y; 655 PetscInt i, m, M, N; 656 ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE; 657 658 PetscFunctionBegin; 659 PetscCheck(osm->n == 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not yet implemented"); 660 PetscCall(MatGetSize(Xin, NULL, &N)); 661 if (osm->pctoouter) { 662 PetscCall(VecGetLocalSize(osm->pcx, &m)); 663 PetscCall(VecGetSize(osm->pcx, &M)); 664 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &O)); 665 for (i = 0; i < N; ++i) { 666 PetscCall(MatDenseGetColumnVecRead(Xin, i, &x)); 667 PetscCall(MatDenseGetColumnVecWrite(O, i, &y)); 668 PetscCall(VecScatterBegin(osm->pctoouter, x, y, INSERT_VALUES, SCATTER_REVERSE)); 669 PetscCall(VecScatterEnd(osm->pctoouter, x, y, INSERT_VALUES, SCATTER_REVERSE)); 670 PetscCall(MatDenseRestoreColumnVecWrite(O, i, &y)); 671 PetscCall(MatDenseRestoreColumnVecRead(Xin, i, &x)); 672 } 673 X = Y = O; 674 } else { 675 X = Xin; 676 Y = Yout; 677 } 678 /* 679 support for limiting the restriction or interpolation only to the inner 680 subdomain values (leaving the other values 0). 681 */ 682 PetscCall(VecGetLocalSize(osm->x[0], &m)); 683 PetscCall(VecGetSize(osm->x[0], &M)); 684 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &Z)); 685 for (i = 0; i < N; ++i) { 686 PetscCall(MatDenseGetColumnVecRead(X, i, &x)); 687 PetscCall(MatDenseGetColumnVecWrite(Z, i, &y)); 688 if (!(osm->type & PC_GASM_RESTRICT)) { 689 /* have to zero the work RHS since scatter may leave some slots empty */ 690 PetscCall(VecZeroEntries(y)); 691 PetscCall(VecScatterBegin(osm->girestriction, x, y, INSERT_VALUES, forward)); 692 PetscCall(VecScatterEnd(osm->girestriction, x, y, INSERT_VALUES, forward)); 693 } else { 694 PetscCall(VecScatterBegin(osm->gorestriction, x, y, INSERT_VALUES, forward)); 695 PetscCall(VecScatterEnd(osm->gorestriction, x, y, INSERT_VALUES, forward)); 696 } 697 PetscCall(MatDenseRestoreColumnVecWrite(Z, i, &y)); 698 PetscCall(MatDenseRestoreColumnVecRead(X, i, &x)); 699 } 700 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &W)); 701 PetscCall(MatSetOption(Z, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 702 PetscCall(MatAssemblyBegin(Z, MAT_FINAL_ASSEMBLY)); 703 PetscCall(MatAssemblyEnd(Z, MAT_FINAL_ASSEMBLY)); 704 /* do the subdomain solve */ 705 PetscCall(KSPMatSolve(osm->ksp[0], Z, W)); 706 PetscCall(KSPCheckSolve(osm->ksp[0], pc, NULL)); 707 PetscCall(MatDestroy(&Z)); 708 /* do we need to zero y? */ 709 PetscCall(MatZeroEntries(Y)); 710 for (i = 0; i < N; ++i) { 711 PetscCall(MatDenseGetColumnVecWrite(Y, i, &y)); 712 PetscCall(MatDenseGetColumnVecRead(W, i, &x)); 713 if (!(osm->type & PC_GASM_INTERPOLATE)) { 714 PetscCall(VecScatterBegin(osm->girestriction, x, y, ADD_VALUES, reverse)); 715 PetscCall(VecScatterEnd(osm->girestriction, x, y, ADD_VALUES, reverse)); 716 } else { 717 PetscCall(VecScatterBegin(osm->gorestriction, x, y, ADD_VALUES, reverse)); 718 PetscCall(VecScatterEnd(osm->gorestriction, x, y, ADD_VALUES, reverse)); 719 } 720 PetscCall(MatDenseRestoreColumnVecRead(W, i, &x)); 721 if (osm->pctoouter) { 722 PetscCall(MatDenseGetColumnVecWrite(Yout, i, &x)); 723 PetscCall(VecScatterBegin(osm->pctoouter, y, x, INSERT_VALUES, SCATTER_FORWARD)); 724 PetscCall(VecScatterEnd(osm->pctoouter, y, x, INSERT_VALUES, SCATTER_FORWARD)); 725 PetscCall(MatDenseRestoreColumnVecRead(Yout, i, &x)); 726 } 727 PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &y)); 728 } 729 PetscCall(MatDestroy(&W)); 730 PetscCall(MatDestroy(&O)); 731 PetscFunctionReturn(PETSC_SUCCESS); 732 } 733 734 static PetscErrorCode PCApplyTranspose_GASM(PC pc, Vec xin, Vec yout) 735 { 736 PC_GASM *osm = (PC_GASM *)pc->data; 737 PetscInt i; 738 Vec x, y; 739 ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE; 740 741 PetscFunctionBegin; 742 if (osm->pctoouter) { 743 PetscCall(VecScatterBegin(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 744 PetscCall(VecScatterEnd(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE)); 745 x = osm->pcx; 746 y = osm->pcy; 747 } else { 748 x = xin; 749 y = yout; 750 } 751 /* 752 Support for limiting the restriction or interpolation to only local 753 subdomain values (leaving the other values 0). 754 755 Note: these are reversed from the PCApply_GASM() because we are applying the 756 transpose of the three terms 757 */ 758 if (!(osm->type & PC_GASM_INTERPOLATE)) { 759 /* have to zero the work RHS since scatter may leave some slots empty */ 760 PetscCall(VecZeroEntries(osm->gx)); 761 PetscCall(VecScatterBegin(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 762 } else { 763 PetscCall(VecScatterBegin(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 764 } 765 PetscCall(VecZeroEntries(osm->gy)); 766 if (!(osm->type & PC_GASM_INTERPOLATE)) { 767 PetscCall(VecScatterEnd(osm->girestriction, x, osm->gx, INSERT_VALUES, forward)); 768 } else { 769 PetscCall(VecScatterEnd(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward)); 770 } 771 /* do the local solves */ 772 for (i = 0; i < osm->n; ++i) { /* Note that the solves are local, so we can go to osm->n, rather than osm->nmax. */ 773 PetscCall(KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i])); 774 PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i])); 775 } 776 PetscCall(VecZeroEntries(y)); 777 if (!(osm->type & PC_GASM_RESTRICT)) { 778 PetscCall(VecScatterBegin(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 779 PetscCall(VecScatterEnd(osm->girestriction, osm->gy, y, ADD_VALUES, reverse)); 780 } else { 781 PetscCall(VecScatterBegin(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 782 PetscCall(VecScatterEnd(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse)); 783 } 784 if (osm->pctoouter) { 785 PetscCall(VecScatterBegin(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 786 PetscCall(VecScatterEnd(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD)); 787 } 788 PetscFunctionReturn(PETSC_SUCCESS); 789 } 790 791 static PetscErrorCode PCReset_GASM(PC pc) 792 { 793 PC_GASM *osm = (PC_GASM *)pc->data; 794 PetscInt i; 795 796 PetscFunctionBegin; 797 if (osm->ksp) { 798 for (i = 0; i < osm->n; i++) PetscCall(KSPReset(osm->ksp[i])); 799 } 800 if (osm->pmat) { 801 if (osm->n > 0) { 802 PetscMPIInt size; 803 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 804 if (size > 1) { 805 /* osm->pmat is created by MatCreateSubMatricesMPI(), cannot use MatDestroySubMatrices() */ 806 PetscCall(MatDestroyMatrices(osm->n, &osm->pmat)); 807 } else { 808 PetscCall(MatDestroySubMatrices(osm->n, &osm->pmat)); 809 } 810 } 811 } 812 if (osm->x) { 813 for (i = 0; i < osm->n; i++) { 814 PetscCall(VecDestroy(&osm->x[i])); 815 PetscCall(VecDestroy(&osm->y[i])); 816 } 817 } 818 PetscCall(VecDestroy(&osm->gx)); 819 PetscCall(VecDestroy(&osm->gy)); 820 821 PetscCall(VecScatterDestroy(&osm->gorestriction)); 822 PetscCall(VecScatterDestroy(&osm->girestriction)); 823 if (!osm->user_subdomains) { 824 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->ois, &osm->iis)); 825 osm->N = PETSC_DETERMINE; 826 osm->nmax = PETSC_DETERMINE; 827 } 828 if (osm->pctoouter) PetscCall(VecScatterDestroy(&(osm->pctoouter))); 829 if (osm->permutationIS) PetscCall(ISDestroy(&(osm->permutationIS))); 830 if (osm->pcx) PetscCall(VecDestroy(&(osm->pcx))); 831 if (osm->pcy) PetscCall(VecDestroy(&(osm->pcy))); 832 if (osm->permutationP) PetscCall(MatDestroy(&(osm->permutationP))); 833 if (osm->pcmat) PetscCall(MatDestroy(&osm->pcmat)); 834 PetscFunctionReturn(PETSC_SUCCESS); 835 } 836 837 static PetscErrorCode PCDestroy_GASM(PC pc) 838 { 839 PC_GASM *osm = (PC_GASM *)pc->data; 840 PetscInt i; 841 842 PetscFunctionBegin; 843 PetscCall(PCReset_GASM(pc)); 844 /* PCReset will not destroy subdomains, if user_subdomains is true. */ 845 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->ois, &osm->iis)); 846 if (osm->ksp) { 847 for (i = 0; i < osm->n; i++) PetscCall(KSPDestroy(&osm->ksp[i])); 848 PetscCall(PetscFree(osm->ksp)); 849 } 850 PetscCall(PetscFree(osm->x)); 851 PetscCall(PetscFree(osm->y)); 852 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", NULL)); 853 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", NULL)); 854 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", NULL)); 855 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", NULL)); 856 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", NULL)); 857 PetscCall(PetscFree(pc->data)); 858 PetscFunctionReturn(PETSC_SUCCESS); 859 } 860 861 static PetscErrorCode PCSetFromOptions_GASM(PC pc, PetscOptionItems *PetscOptionsObject) 862 { 863 PC_GASM *osm = (PC_GASM *)pc->data; 864 PetscInt blocks, ovl; 865 PetscBool flg; 866 PCGASMType gasmtype; 867 868 PetscFunctionBegin; 869 PetscOptionsHeadBegin(PetscOptionsObject, "Generalized additive Schwarz options"); 870 PetscCall(PetscOptionsBool("-pc_gasm_use_dm_subdomains", "If subdomains aren't set, use DMCreateDomainDecomposition() to define subdomains.", "PCGASMSetUseDMSubdomains", osm->dm_subdomains, &osm->dm_subdomains, &flg)); 871 PetscCall(PetscOptionsInt("-pc_gasm_total_subdomains", "Total number of subdomains across communicator", "PCGASMSetTotalSubdomains", osm->N, &blocks, &flg)); 872 if (flg) PetscCall(PCGASMSetTotalSubdomains(pc, blocks)); 873 PetscCall(PetscOptionsInt("-pc_gasm_overlap", "Number of overlapping degrees of freedom", "PCGASMSetOverlap", osm->overlap, &ovl, &flg)); 874 if (flg) { 875 PetscCall(PCGASMSetOverlap(pc, ovl)); 876 osm->dm_subdomains = PETSC_FALSE; 877 } 878 flg = PETSC_FALSE; 879 PetscCall(PetscOptionsEnum("-pc_gasm_type", "Type of restriction/extension", "PCGASMSetType", PCGASMTypes, (PetscEnum)osm->type, (PetscEnum *)&gasmtype, &flg)); 880 if (flg) PetscCall(PCGASMSetType(pc, gasmtype)); 881 PetscCall(PetscOptionsBool("-pc_gasm_use_hierachical_partitioning", "use hierarchical partitioning", NULL, osm->hierarchicalpartitioning, &osm->hierarchicalpartitioning, &flg)); 882 PetscOptionsHeadEnd(); 883 PetscFunctionReturn(PETSC_SUCCESS); 884 } 885 886 /*@ 887 PCGASMSetTotalSubdomains - sets the total number of subdomains to use across the communicator for `PCGASM` 888 889 Logically Collective 890 891 Input Parameters: 892 + pc - the preconditioner 893 - N - total number of subdomains 894 895 Level: beginner 896 897 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()` 898 `PCGASMCreateSubdomains2D()` 899 @*/ 900 PetscErrorCode PCGASMSetTotalSubdomains(PC pc, PetscInt N) 901 { 902 PC_GASM *osm = (PC_GASM *)pc->data; 903 PetscMPIInt size, rank; 904 905 PetscFunctionBegin; 906 PetscCheck(N >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Total number of subdomains must be 1 or more, got N = %" PetscInt_FMT, N); 907 PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetTotalSubdomains() should be called before calling PCSetUp()."); 908 909 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois)); 910 osm->ois = osm->iis = NULL; 911 912 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); 913 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank)); 914 osm->N = N; 915 osm->n = PETSC_DETERMINE; 916 osm->nmax = PETSC_DETERMINE; 917 osm->dm_subdomains = PETSC_FALSE; 918 PetscFunctionReturn(PETSC_SUCCESS); 919 } 920 921 static PetscErrorCode PCGASMSetSubdomains_GASM(PC pc, PetscInt n, IS iis[], IS ois[]) 922 { 923 PC_GASM *osm = (PC_GASM *)pc->data; 924 PetscInt i; 925 926 PetscFunctionBegin; 927 PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Each MPI rank must have 1 or more subdomains, got n = %" PetscInt_FMT, n); 928 PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetSubdomains() should be called before calling PCSetUp()."); 929 930 PetscCall(PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois)); 931 osm->iis = osm->ois = NULL; 932 osm->n = n; 933 osm->N = PETSC_DETERMINE; 934 osm->nmax = PETSC_DETERMINE; 935 if (ois) { 936 PetscCall(PetscMalloc1(n, &osm->ois)); 937 for (i = 0; i < n; i++) { 938 PetscCall(PetscObjectReference((PetscObject)ois[i])); 939 osm->ois[i] = ois[i]; 940 } 941 /* 942 Since the user set the outer subdomains, even if nontrivial overlap was requested via PCGASMSetOverlap(), 943 it will be ignored. To avoid confusion later on (e.g., when viewing the PC), the overlap size is set to -1. 944 */ 945 osm->overlap = -1; 946 /* inner subdomains must be provided */ 947 PetscCheck(iis, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "inner indices have to be provided "); 948 } /* end if */ 949 if (iis) { 950 PetscCall(PetscMalloc1(n, &osm->iis)); 951 for (i = 0; i < n; i++) { 952 PetscCall(PetscObjectReference((PetscObject)iis[i])); 953 osm->iis[i] = iis[i]; 954 } 955 if (!ois) { 956 osm->ois = NULL; 957 /* if user does not provide outer indices, we will create the corresponding outer indices using osm->overlap =1 in PCSetUp_GASM */ 958 } 959 } 960 if (PetscDefined(USE_DEBUG)) { 961 PetscInt j, rstart, rend, *covered, lsize; 962 const PetscInt *indices; 963 /* check if the inner indices cover and only cover the local portion of the preconditioning matrix */ 964 PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend)); 965 PetscCall(PetscCalloc1(rend - rstart, &covered)); 966 /* check if the current MPI rank owns indices from others */ 967 for (i = 0; i < n; i++) { 968 PetscCall(ISGetIndices(osm->iis[i], &indices)); 969 PetscCall(ISGetLocalSize(osm->iis[i], &lsize)); 970 for (j = 0; j < lsize; j++) { 971 PetscCheck(indices[j] >= rstart && indices[j] < rend, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "inner subdomains can not own an index %" PetscInt_FMT " from other ranks", indices[j]); 972 PetscCheck(covered[indices[j] - rstart] != 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "inner subdomains can not have an overlapping index %" PetscInt_FMT " ", indices[j]); 973 covered[indices[j] - rstart] = 1; 974 } 975 PetscCall(ISRestoreIndices(osm->iis[i], &indices)); 976 } 977 /* check if we miss any indices */ 978 for (i = rstart; i < rend; i++) PetscCheck(covered[i - rstart], PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "local entity %" PetscInt_FMT " was not covered by inner subdomains", i); 979 PetscCall(PetscFree(covered)); 980 } 981 if (iis) osm->user_subdomains = PETSC_TRUE; 982 PetscFunctionReturn(PETSC_SUCCESS); 983 } 984 985 static PetscErrorCode PCGASMSetOverlap_GASM(PC pc, PetscInt ovl) 986 { 987 PC_GASM *osm = (PC_GASM *)pc->data; 988 989 PetscFunctionBegin; 990 PetscCheck(ovl >= 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap value requested"); 991 PetscCheck(!pc->setupcalled || ovl == osm->overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetOverlap() should be called before PCSetUp()."); 992 if (!pc->setupcalled) osm->overlap = ovl; 993 PetscFunctionReturn(PETSC_SUCCESS); 994 } 995 996 static PetscErrorCode PCGASMSetType_GASM(PC pc, PCGASMType type) 997 { 998 PC_GASM *osm = (PC_GASM *)pc->data; 999 1000 PetscFunctionBegin; 1001 osm->type = type; 1002 osm->type_set = PETSC_TRUE; 1003 PetscFunctionReturn(PETSC_SUCCESS); 1004 } 1005 1006 static PetscErrorCode PCGASMSetSortIndices_GASM(PC pc, PetscBool doSort) 1007 { 1008 PC_GASM *osm = (PC_GASM *)pc->data; 1009 1010 PetscFunctionBegin; 1011 osm->sort_indices = doSort; 1012 PetscFunctionReturn(PETSC_SUCCESS); 1013 } 1014 1015 /* 1016 FIXME: This routine might need to be modified now that multiple ranks per subdomain are allowed. 1017 In particular, it would upset the global subdomain number calculation. 1018 */ 1019 static PetscErrorCode PCGASMGetSubKSP_GASM(PC pc, PetscInt *n, PetscInt *first, KSP **ksp) 1020 { 1021 PC_GASM *osm = (PC_GASM *)pc->data; 1022 1023 PetscFunctionBegin; 1024 PetscCheck(osm->n >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Need to call PCSetUp() on PC (or KSPSetUp() on the outer KSP object) before calling here"); 1025 1026 if (n) *n = osm->n; 1027 if (first) { 1028 PetscCallMPI(MPI_Scan(&osm->n, first, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc))); 1029 *first -= osm->n; 1030 } 1031 if (ksp) { 1032 /* Assume that local solves are now different; not necessarily 1033 true, though! This flag is used only for PCView_GASM() */ 1034 *ksp = osm->ksp; 1035 osm->same_subdomain_solvers = PETSC_FALSE; 1036 } 1037 PetscFunctionReturn(PETSC_SUCCESS); 1038 } /* PCGASMGetSubKSP_GASM() */ 1039 1040 /*@C 1041 PCGASMSetSubdomains - Sets the subdomains for this MPI rank 1042 for the additive Schwarz preconditioner with multiple MPI ranks per subdomain, `PCGASM` 1043 1044 Collective 1045 1046 Input Parameters: 1047 + pc - the preconditioner object 1048 . n - the number of subdomains for this MPI rank 1049 . iis - the index sets that define the inner subdomains (or `NULL` for PETSc to determine subdomains) 1050 - ois - the index sets that define the outer subdomains (or `NULL` to use the same as `iis`, or to construct by expanding `iis` by 1051 the requested overlap) 1052 1053 Level: advanced 1054 1055 Notes: 1056 The `IS` indices use the parallel, global numbering of the vector entries. 1057 1058 Inner subdomains are those where the correction is applied. 1059 1060 Outer subdomains are those where the residual necessary to obtain the 1061 corrections is obtained (see `PCGASMType` for the use of inner/outer subdomains). 1062 1063 Both inner and outer subdomains can extend over several MPI ranks. 1064 This rank's portion of a subdomain is known as a local subdomain. 1065 1066 Inner subdomains can not overlap with each other, do not have any entities from remote ranks, 1067 and have to cover the entire local subdomain owned by the current rank. The index sets on each 1068 rank should be ordered such that the ith local subdomain is connected to the ith remote subdomain 1069 on another MPI rank. 1070 1071 By default the `PGASM` preconditioner uses 1 (local) subdomain per MPI rank. 1072 1073 The `iis` and `ois` arrays may be freed after this call using `PCGASMDestroySubdomains()` 1074 1075 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, `PCGASMDestroySubdomains()`, 1076 `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()` 1077 @*/ 1078 PetscErrorCode PCGASMSetSubdomains(PC pc, PetscInt n, IS iis[], IS ois[]) 1079 { 1080 PC_GASM *osm = (PC_GASM *)pc->data; 1081 1082 PetscFunctionBegin; 1083 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1084 PetscTryMethod(pc, "PCGASMSetSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, n, iis, ois)); 1085 osm->dm_subdomains = PETSC_FALSE; 1086 PetscFunctionReturn(PETSC_SUCCESS); 1087 } 1088 1089 /*@ 1090 PCGASMSetOverlap - Sets the overlap between a pair of subdomains for the 1091 additive Schwarz preconditioner `PCGASM`. Either all or no MPI ranks in the 1092 pc communicator must call this routine. 1093 1094 Logically Collective 1095 1096 Input Parameters: 1097 + pc - the preconditioner context 1098 - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 0) 1099 1100 Options Database Key: 1101 . -pc_gasm_overlap <overlap> - Sets overlap 1102 1103 Level: intermediate 1104 1105 Notes: 1106 By default the `PCGASM` preconditioner uses 1 subdomain per rank. To use 1107 multiple subdomain per perocessor or "straddling" subdomains that intersect 1108 multiple ranks use `PCGASMSetSubdomains()` (or option `-pc_gasm_total_subdomains` <n>). 1109 1110 The overlap defaults to 0, so if one desires that no additional 1111 overlap be computed beyond what may have been set with a call to 1112 `PCGASMSetSubdomains()`, then `ovl` must be set to be 0. In particular, if one does 1113 not explicitly set the subdomains in application code, then all overlap would be computed 1114 internally by PETSc, and using an overlap of 0 would result in an `PCGASM` 1115 variant that is equivalent to the block Jacobi preconditioner. 1116 1117 One can define initial index sets with any overlap via 1118 `PCGASMSetSubdomains()`; the routine `PCGASMSetOverlap()` merely allows 1119 PETSc to extend that overlap further, if desired. 1120 1121 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1122 `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()` 1123 @*/ 1124 PetscErrorCode PCGASMSetOverlap(PC pc, PetscInt ovl) 1125 { 1126 PC_GASM *osm = (PC_GASM *)pc->data; 1127 1128 PetscFunctionBegin; 1129 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1130 PetscValidLogicalCollectiveInt(pc, ovl, 2); 1131 PetscTryMethod(pc, "PCGASMSetOverlap_C", (PC, PetscInt), (pc, ovl)); 1132 osm->dm_subdomains = PETSC_FALSE; 1133 PetscFunctionReturn(PETSC_SUCCESS); 1134 } 1135 1136 /*@ 1137 PCGASMSetType - Sets the type of restriction and interpolation used 1138 for local problems in the `PCGASM` additive Schwarz method. 1139 1140 Logically Collective 1141 1142 Input Parameters: 1143 + pc - the preconditioner context 1144 - type - variant of `PCGASM`, one of 1145 .vb 1146 `PC_GASM_BASIC` - full interpolation and restriction 1147 `PC_GASM_RESTRICT` - full restriction, local MPI rank interpolation 1148 `PC_GASM_INTERPOLATE` - full interpolation, local MPI rank restriction 1149 `PC_GASM_NONE` - local MPI rank restriction and interpolation 1150 .ve 1151 1152 Options Database Key: 1153 . -pc_gasm_type [basic,restrict,interpolate,none] - Sets `PCGASM` type 1154 1155 Level: intermediate 1156 1157 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1158 `PCGASMCreateSubdomains2D()`, `PCASM`, `PCASMSetType()` 1159 @*/ 1160 PetscErrorCode PCGASMSetType(PC pc, PCGASMType type) 1161 { 1162 PetscFunctionBegin; 1163 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1164 PetscValidLogicalCollectiveEnum(pc, type, 2); 1165 PetscTryMethod(pc, "PCGASMSetType_C", (PC, PCGASMType), (pc, type)); 1166 PetscFunctionReturn(PETSC_SUCCESS); 1167 } 1168 1169 /*@ 1170 PCGASMSetSortIndices - Determines whether subdomain indices are sorted. 1171 1172 Logically Collective 1173 1174 Input Parameters: 1175 + pc - the preconditioner context 1176 - doSort - sort the subdomain indices 1177 1178 Level: intermediate 1179 1180 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, 1181 `PCGASMCreateSubdomains2D()` 1182 @*/ 1183 PetscErrorCode PCGASMSetSortIndices(PC pc, PetscBool doSort) 1184 { 1185 PetscFunctionBegin; 1186 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1187 PetscValidLogicalCollectiveBool(pc, doSort, 2); 1188 PetscTryMethod(pc, "PCGASMSetSortIndices_C", (PC, PetscBool), (pc, doSort)); 1189 PetscFunctionReturn(PETSC_SUCCESS); 1190 } 1191 1192 /*@C 1193 PCGASMGetSubKSP - Gets the local `KSP` contexts for all subdomains on 1194 this MPI rank. 1195 1196 Collective iff first_local is requested 1197 1198 Input Parameter: 1199 . pc - the preconditioner context 1200 1201 Output Parameters: 1202 + n_local - the number of blocks on this MPI rank or `NULL` 1203 . first_local - the global number of the first block on this rank or `NULL`, 1204 all ranks must request or all must pass `NULL` 1205 - ksp - the array of `KSP` contexts 1206 1207 Level: advanced 1208 1209 Note: 1210 After `PCGASMGetSubKSP()` the array of `KSP`es is not to be freed 1211 1212 Currently for some matrix implementations only 1 block per MPI process 1213 is supported. 1214 1215 You must call `KSPSetUp()` before calling `PCGASMGetSubKSP()`. 1216 1217 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()`, 1218 `PCGASMCreateSubdomains2D()`, 1219 @*/ 1220 PetscErrorCode PCGASMGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[]) 1221 { 1222 PetscFunctionBegin; 1223 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1224 PetscUseMethod(pc, "PCGASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp)); 1225 PetscFunctionReturn(PETSC_SUCCESS); 1226 } 1227 1228 /*MC 1229 PCGASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with 1230 its own `KSP` object on a subset of MPI ranks 1231 1232 Options Database Keys: 1233 + -pc_gasm_total_subdomains <n> - Sets total number of local subdomains to be distributed among MPI ranks 1234 . -pc_gasm_view_subdomains - activates the printing of subdomain indices in `PCView()`, -ksp_view or -snes_view 1235 . -pc_gasm_print_subdomains - activates the printing of subdomain indices in `PCSetUp()` 1236 . -pc_gasm_overlap <ovl> - Sets overlap by which to (automatically) extend local subdomains 1237 - -pc_gasm_type [basic,restrict,interpolate,none] - Sets `PCGASMType` 1238 1239 Level: beginner 1240 1241 Notes: 1242 To set options on the solvers for each block append `-sub_` to all the `KSP`, and `PC` 1243 options database keys. For example, `-sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly` 1244 1245 To set the options on the solvers separate for each block call `PCGASMGetSubKSP()` 1246 and set the options directly on the resulting `KSP` object (you can access its `PC` 1247 with `KSPGetPC()`) 1248 1249 References: 1250 + * - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions 1251 Courant Institute, New York University Technical report 1252 - * - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations, 1253 Cambridge University Press. 1254 1255 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCASM`, `PCGASMType`, `PCGASMSetType()`, 1256 `PCBJACOBI`, `PCGASMGetSubKSP()`, `PCGASMSetSubdomains()`, 1257 `PCSetModifySubMatrices()`, `PCGASMSetOverlap()`, `PCGASMSetType()` 1258 M*/ 1259 1260 PETSC_EXTERN PetscErrorCode PCCreate_GASM(PC pc) 1261 { 1262 PC_GASM *osm; 1263 1264 PetscFunctionBegin; 1265 PetscCall(PetscNew(&osm)); 1266 1267 osm->N = PETSC_DETERMINE; 1268 osm->n = PETSC_DECIDE; 1269 osm->nmax = PETSC_DETERMINE; 1270 osm->overlap = 0; 1271 osm->ksp = NULL; 1272 osm->gorestriction = NULL; 1273 osm->girestriction = NULL; 1274 osm->pctoouter = NULL; 1275 osm->gx = NULL; 1276 osm->gy = NULL; 1277 osm->x = NULL; 1278 osm->y = NULL; 1279 osm->pcx = NULL; 1280 osm->pcy = NULL; 1281 osm->permutationIS = NULL; 1282 osm->permutationP = NULL; 1283 osm->pcmat = NULL; 1284 osm->ois = NULL; 1285 osm->iis = NULL; 1286 osm->pmat = NULL; 1287 osm->type = PC_GASM_RESTRICT; 1288 osm->same_subdomain_solvers = PETSC_TRUE; 1289 osm->sort_indices = PETSC_TRUE; 1290 osm->dm_subdomains = PETSC_FALSE; 1291 osm->hierarchicalpartitioning = PETSC_FALSE; 1292 1293 pc->data = (void *)osm; 1294 pc->ops->apply = PCApply_GASM; 1295 pc->ops->matapply = PCMatApply_GASM; 1296 pc->ops->applytranspose = PCApplyTranspose_GASM; 1297 pc->ops->setup = PCSetUp_GASM; 1298 pc->ops->reset = PCReset_GASM; 1299 pc->ops->destroy = PCDestroy_GASM; 1300 pc->ops->setfromoptions = PCSetFromOptions_GASM; 1301 pc->ops->setuponblocks = PCSetUpOnBlocks_GASM; 1302 pc->ops->view = PCView_GASM; 1303 pc->ops->applyrichardson = NULL; 1304 1305 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", PCGASMSetSubdomains_GASM)); 1306 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", PCGASMSetOverlap_GASM)); 1307 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", PCGASMSetType_GASM)); 1308 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", PCGASMSetSortIndices_GASM)); 1309 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", PCGASMGetSubKSP_GASM)); 1310 PetscFunctionReturn(PETSC_SUCCESS); 1311 } 1312 1313 PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[]) 1314 { 1315 MatPartitioning mpart; 1316 const char *prefix; 1317 PetscInt i, j, rstart, rend, bs; 1318 PetscBool hasop, isbaij = PETSC_FALSE, foundpart = PETSC_FALSE; 1319 Mat Ad = NULL, adj; 1320 IS ispart, isnumb, *is; 1321 1322 PetscFunctionBegin; 1323 PetscCheck(nloc >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "number of local subdomains must > 0, got nloc = %" PetscInt_FMT, nloc); 1324 1325 /* Get prefix, row distribution, and block size */ 1326 PetscCall(MatGetOptionsPrefix(A, &prefix)); 1327 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 1328 PetscCall(MatGetBlockSize(A, &bs)); 1329 PetscCheck(rstart / bs * bs == rstart && rend / bs * bs == rend, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "bad row distribution [%" PetscInt_FMT ",%" PetscInt_FMT ") for matrix block size %" PetscInt_FMT, rstart, rend, bs); 1330 1331 /* Get diagonal block from matrix if possible */ 1332 PetscCall(MatHasOperation(A, MATOP_GET_DIAGONAL_BLOCK, &hasop)); 1333 if (hasop) PetscCall(MatGetDiagonalBlock(A, &Ad)); 1334 if (Ad) { 1335 PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQBAIJ, &isbaij)); 1336 if (!isbaij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQSBAIJ, &isbaij)); 1337 } 1338 if (Ad && nloc > 1) { 1339 PetscBool match, done; 1340 /* Try to setup a good matrix partitioning if available */ 1341 PetscCall(MatPartitioningCreate(PETSC_COMM_SELF, &mpart)); 1342 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix)); 1343 PetscCall(MatPartitioningSetFromOptions(mpart)); 1344 PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGCURRENT, &match)); 1345 if (!match) PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGSQUARE, &match)); 1346 if (!match) { /* assume a "good" partitioner is available */ 1347 PetscInt na; 1348 const PetscInt *ia, *ja; 1349 PetscCall(MatGetRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done)); 1350 if (done) { 1351 /* Build adjacency matrix by hand. Unfortunately a call to 1352 MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will 1353 remove the block-aij structure and we cannot expect 1354 MatPartitioning to split vertices as we need */ 1355 PetscInt i, j, len, nnz, cnt, *iia = NULL, *jja = NULL; 1356 const PetscInt *row; 1357 nnz = 0; 1358 for (i = 0; i < na; i++) { /* count number of nonzeros */ 1359 len = ia[i + 1] - ia[i]; 1360 row = ja + ia[i]; 1361 for (j = 0; j < len; j++) { 1362 if (row[j] == i) { /* don't count diagonal */ 1363 len--; 1364 break; 1365 } 1366 } 1367 nnz += len; 1368 } 1369 PetscCall(PetscMalloc1(na + 1, &iia)); 1370 PetscCall(PetscMalloc1(nnz, &jja)); 1371 nnz = 0; 1372 iia[0] = 0; 1373 for (i = 0; i < na; i++) { /* fill adjacency */ 1374 cnt = 0; 1375 len = ia[i + 1] - ia[i]; 1376 row = ja + ia[i]; 1377 for (j = 0; j < len; j++) { 1378 if (row[j] != i) jja[nnz + cnt++] = row[j]; /* if not diagonal */ 1379 } 1380 nnz += cnt; 1381 iia[i + 1] = nnz; 1382 } 1383 /* Partitioning of the adjacency matrix */ 1384 PetscCall(MatCreateMPIAdj(PETSC_COMM_SELF, na, na, iia, jja, NULL, &adj)); 1385 PetscCall(MatPartitioningSetAdjacency(mpart, adj)); 1386 PetscCall(MatPartitioningSetNParts(mpart, nloc)); 1387 PetscCall(MatPartitioningApply(mpart, &ispart)); 1388 PetscCall(ISPartitioningToNumbering(ispart, &isnumb)); 1389 PetscCall(MatDestroy(&adj)); 1390 foundpart = PETSC_TRUE; 1391 } 1392 PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done)); 1393 } 1394 PetscCall(MatPartitioningDestroy(&mpart)); 1395 } 1396 PetscCall(PetscMalloc1(nloc, &is)); 1397 if (!foundpart) { 1398 /* Partitioning by contiguous chunks of rows */ 1399 1400 PetscInt mbs = (rend - rstart) / bs; 1401 PetscInt start = rstart; 1402 for (i = 0; i < nloc; i++) { 1403 PetscInt count = (mbs / nloc + ((mbs % nloc) > i)) * bs; 1404 PetscCall(ISCreateStride(PETSC_COMM_SELF, count, start, 1, &is[i])); 1405 start += count; 1406 } 1407 1408 } else { 1409 /* Partitioning by adjacency of diagonal block */ 1410 1411 const PetscInt *numbering; 1412 PetscInt *count, nidx, *indices, *newidx, start = 0; 1413 /* Get node count in each partition */ 1414 PetscCall(PetscMalloc1(nloc, &count)); 1415 PetscCall(ISPartitioningCount(ispart, nloc, count)); 1416 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1417 for (i = 0; i < nloc; i++) count[i] *= bs; 1418 } 1419 /* Build indices from node numbering */ 1420 PetscCall(ISGetLocalSize(isnumb, &nidx)); 1421 PetscCall(PetscMalloc1(nidx, &indices)); 1422 for (i = 0; i < nidx; i++) indices[i] = i; /* needs to be initialized */ 1423 PetscCall(ISGetIndices(isnumb, &numbering)); 1424 PetscCall(PetscSortIntWithPermutation(nidx, numbering, indices)); 1425 PetscCall(ISRestoreIndices(isnumb, &numbering)); 1426 if (isbaij && bs > 1) { /* adjust for the block-aij case */ 1427 PetscCall(PetscMalloc1(nidx * bs, &newidx)); 1428 for (i = 0; i < nidx; i++) { 1429 for (j = 0; j < bs; j++) newidx[i * bs + j] = indices[i] * bs + j; 1430 } 1431 PetscCall(PetscFree(indices)); 1432 nidx *= bs; 1433 indices = newidx; 1434 } 1435 /* Shift to get global indices */ 1436 for (i = 0; i < nidx; i++) indices[i] += rstart; 1437 1438 /* Build the index sets for each block */ 1439 for (i = 0; i < nloc; i++) { 1440 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, count[i], &indices[start], PETSC_COPY_VALUES, &is[i])); 1441 PetscCall(ISSort(is[i])); 1442 start += count[i]; 1443 } 1444 1445 PetscCall(PetscFree(count)); 1446 PetscCall(PetscFree(indices)); 1447 PetscCall(ISDestroy(&isnumb)); 1448 PetscCall(ISDestroy(&ispart)); 1449 } 1450 *iis = is; 1451 PetscFunctionReturn(PETSC_SUCCESS); 1452 } 1453 1454 PETSC_INTERN PetscErrorCode PCGASMCreateStraddlingSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[]) 1455 { 1456 PetscFunctionBegin; 1457 PetscCall(MatSubdomainsCreateCoalesce(A, N, n, iis)); 1458 PetscFunctionReturn(PETSC_SUCCESS); 1459 } 1460 1461 /*@C 1462 PCGASMCreateSubdomains - Creates `n` index sets defining `n` nonoverlapping subdomains on this MPI process for the `PCGASM` additive 1463 Schwarz preconditioner for a any problem based on its matrix. 1464 1465 Collective 1466 1467 Input Parameters: 1468 + A - The global matrix operator 1469 - N - the number of global subdomains requested 1470 1471 Output Parameters: 1472 + n - the number of subdomains created on this MPI rank 1473 - iis - the array of index sets defining the local inner subdomains (on which the correction is applied) 1474 1475 Level: advanced 1476 1477 Notes: 1478 When `N` >= A's communicator size, each subdomain is local -- contained within a single MPI process. 1479 When `N` < size, the subdomains are 'straddling' (rank boundaries) and are no longer local. 1480 The resulting subdomains can be use in `PCGASMSetSubdomains`(pc,n,iss,`NULL`). The overlapping 1481 outer subdomains will be automatically generated from these according to the requested amount of 1482 overlap; this is currently supported only with local subdomains. 1483 1484 Use `PCGASMDestroySubdomains()` to free the array and the list of index sets. 1485 1486 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMDestroySubdomains()` 1487 @*/ 1488 PetscErrorCode PCGASMCreateSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[]) 1489 { 1490 PetscMPIInt size; 1491 1492 PetscFunctionBegin; 1493 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1494 PetscAssertPointer(iis, 4); 1495 1496 PetscCheck(N >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of subdomains must be > 0, N = %" PetscInt_FMT, N); 1497 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 1498 if (N >= size) { 1499 *n = N / size + (N % size); 1500 PetscCall(PCGASMCreateLocalSubdomains(A, *n, iis)); 1501 } else { 1502 PetscCall(PCGASMCreateStraddlingSubdomains(A, N, n, iis)); 1503 } 1504 PetscFunctionReturn(PETSC_SUCCESS); 1505 } 1506 1507 /*@C 1508 PCGASMDestroySubdomains - Destroys the index sets created with 1509 `PCGASMCreateSubdomains()` or `PCGASMCreateSubdomains2D()`. Should be 1510 called after setting subdomains with `PCGASMSetSubdomains()`. 1511 1512 Collective 1513 1514 Input Parameters: 1515 + n - the number of index sets 1516 . iis - the array of inner subdomains 1517 - ois - the array of outer subdomains, can be `NULL` 1518 1519 Level: intermediate 1520 1521 Note: 1522 This is a convenience subroutine that walks each list, 1523 destroys each `IS` on the list, and then frees the list. At the end the 1524 list pointers are set to `NULL`. 1525 1526 .seealso: [](ch_ksp), `PCGASM`, `PCGASMCreateSubdomains()`, `PCGASMSetSubdomains()` 1527 @*/ 1528 PetscErrorCode PCGASMDestroySubdomains(PetscInt n, IS **iis, IS **ois) 1529 { 1530 PetscInt i; 1531 1532 PetscFunctionBegin; 1533 if (n <= 0) PetscFunctionReturn(PETSC_SUCCESS); 1534 if (ois) { 1535 PetscAssertPointer(ois, 3); 1536 if (*ois) { 1537 PetscAssertPointer(*ois, 3); 1538 for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*ois)[i])); 1539 PetscCall(PetscFree((*ois))); 1540 } 1541 } 1542 if (iis) { 1543 PetscAssertPointer(iis, 2); 1544 if (*iis) { 1545 PetscAssertPointer(*iis, 2); 1546 for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*iis)[i])); 1547 PetscCall(PetscFree((*iis))); 1548 } 1549 } 1550 PetscFunctionReturn(PETSC_SUCCESS); 1551 } 1552 1553 #define PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, xleft_loc, ylow_loc, xright_loc, yhigh_loc, n) \ 1554 do { \ 1555 PetscInt first_row = first / M, last_row = last / M + 1; \ 1556 /* \ 1557 Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \ 1558 of the bounding box of the intersection of the subdomain with the local ownership range (local \ 1559 subdomain). \ 1560 Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \ 1561 of the intersection. \ 1562 */ \ 1563 /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \ 1564 *ylow_loc = PetscMax(first_row, ylow); \ 1565 /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1566 *xleft_loc = *ylow_loc == first_row ? PetscMax(first % M, xleft) : xleft; \ 1567 /* yhigh_loc is the grid row above the last local subdomain element */ \ 1568 *yhigh_loc = PetscMin(last_row, yhigh); \ 1569 /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \ 1570 *xright_loc = *yhigh_loc == last_row ? PetscMin(xright, last % M) : xright; \ 1571 /* Now compute the size of the local subdomain n. */ \ 1572 *n = 0; \ 1573 if (*ylow_loc < *yhigh_loc) { \ 1574 PetscInt width = xright - xleft; \ 1575 *n += width * (*yhigh_loc - *ylow_loc - 1); \ 1576 *n += PetscMin(PetscMax(*xright_loc - xleft, 0), width); \ 1577 *n -= PetscMin(PetscMax(*xleft_loc - xleft, 0), width); \ 1578 } \ 1579 } while (0) 1580 1581 /*@C 1582 PCGASMCreateSubdomains2D - Creates the index sets for the `PCGASM` overlapping Schwarz 1583 preconditioner for a two-dimensional problem on a regular grid. 1584 1585 Collective 1586 1587 Input Parameters: 1588 + pc - the preconditioner context 1589 . M - the global number of grid points in the x direction 1590 . N - the global number of grid points in the y direction 1591 . Mdomains - the global number of subdomains in the x direction 1592 . Ndomains - the global number of subdomains in the y direction 1593 . dof - degrees of freedom per node 1594 - overlap - overlap in mesh lines 1595 1596 Output Parameters: 1597 + nsub - the number of local subdomains created 1598 . iis - array of index sets defining inner (nonoverlapping) subdomains 1599 - ois - array of index sets defining outer (overlapping, if overlap > 0) subdomains 1600 1601 Level: advanced 1602 1603 Note: 1604 Use `PCGASMDestroySubdomains()` to free the index sets and the arrays 1605 1606 Fortran Notes: 1607 The `IS` must be declared as an array of length long enough to hold `Nsub` entries 1608 1609 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, `PCGASMSetOverlap()`, `PCASMCreateSubdomains2D()`, 1610 `PCGASMDestroySubdomains()` 1611 @*/ 1612 PetscErrorCode PCGASMCreateSubdomains2D(PC pc, PetscInt M, PetscInt N, PetscInt Mdomains, PetscInt Ndomains, PetscInt dof, PetscInt overlap, PetscInt *nsub, IS **iis, IS **ois) 1613 { 1614 PetscMPIInt size, rank; 1615 PetscInt i, j; 1616 PetscInt maxheight, maxwidth; 1617 PetscInt xstart, xleft, xright, xleft_loc, xright_loc; 1618 PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc; 1619 PetscInt x[2][2], y[2][2], n[2]; 1620 PetscInt first, last; 1621 PetscInt nidx, *idx; 1622 PetscInt ii, jj, s, q, d; 1623 PetscInt k, kk; 1624 PetscMPIInt color; 1625 MPI_Comm comm, subcomm; 1626 IS **xis = NULL, **is = ois, **is_local = iis; 1627 1628 PetscFunctionBegin; 1629 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm)); 1630 PetscCallMPI(MPI_Comm_size(comm, &size)); 1631 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1632 PetscCall(MatGetOwnershipRange(pc->pmat, &first, &last)); 1633 PetscCheck((first % dof) == 0 && (last % dof) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, 1634 "Matrix row partitioning unsuitable for domain decomposition: local row range (%" PetscInt_FMT ",%" PetscInt_FMT ") " 1635 "does not respect the number of degrees of freedom per grid point %" PetscInt_FMT, 1636 first, last, dof); 1637 1638 /* Determine the number of domains with nonzero intersections with the local ownership range. */ 1639 s = 0; 1640 ystart = 0; 1641 for (j = 0; j < Ndomains; ++j) { 1642 maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1643 PetscCheck(maxheight >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the vertical direction for mesh height %" PetscInt_FMT, Ndomains, N); 1644 /* Vertical domain limits with an overlap. */ 1645 ylow = PetscMax(ystart - overlap, 0); 1646 yhigh = PetscMin(ystart + maxheight + overlap, N); 1647 xstart = 0; 1648 for (i = 0; i < Mdomains; ++i) { 1649 maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1650 PetscCheck(maxwidth >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the horizontal direction for mesh width %" PetscInt_FMT, Mdomains, M); 1651 /* Horizontal domain limits with an overlap. */ 1652 xleft = PetscMax(xstart - overlap, 0); 1653 xright = PetscMin(xstart + maxwidth + overlap, M); 1654 /* 1655 Determine whether this subdomain intersects this rank's ownership range of pc->pmat. 1656 */ 1657 PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx)); 1658 if (nidx) ++s; 1659 xstart += maxwidth; 1660 } /* for (i = 0; i < Mdomains; ++i) */ 1661 ystart += maxheight; 1662 } /* for (j = 0; j < Ndomains; ++j) */ 1663 1664 /* Now we can allocate the necessary number of ISs. */ 1665 *nsub = s; 1666 PetscCall(PetscMalloc1(*nsub, is)); 1667 PetscCall(PetscMalloc1(*nsub, is_local)); 1668 s = 0; 1669 ystart = 0; 1670 for (j = 0; j < Ndomains; ++j) { 1671 maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */ 1672 PetscCheck(maxheight >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the vertical direction for mesh height %" PetscInt_FMT, Ndomains, N); 1673 /* Vertical domain limits with an overlap. */ 1674 y[0][0] = PetscMax(ystart - overlap, 0); 1675 y[0][1] = PetscMin(ystart + maxheight + overlap, N); 1676 /* Vertical domain limits without an overlap. */ 1677 y[1][0] = ystart; 1678 y[1][1] = PetscMin(ystart + maxheight, N); 1679 xstart = 0; 1680 for (i = 0; i < Mdomains; ++i) { 1681 maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */ 1682 PetscCheck(maxwidth >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the horizontal direction for mesh width %" PetscInt_FMT, Mdomains, M); 1683 /* Horizontal domain limits with an overlap. */ 1684 x[0][0] = PetscMax(xstart - overlap, 0); 1685 x[0][1] = PetscMin(xstart + maxwidth + overlap, M); 1686 /* Horizontal domain limits without an overlap. */ 1687 x[1][0] = xstart; 1688 x[1][1] = PetscMin(xstart + maxwidth, M); 1689 /* 1690 Determine whether this domain intersects this rank's ownership range of pc->pmat. 1691 Do this twice: first for the domains with overlaps, and once without. 1692 During the first pass create the subcommunicators, and use them on the second pass as well. 1693 */ 1694 for (q = 0; q < 2; ++q) { 1695 PetscBool split = PETSC_FALSE; 1696 /* 1697 domain limits, (xleft, xright) and (ylow, yheigh) are adjusted 1698 according to whether the domain with an overlap or without is considered. 1699 */ 1700 xleft = x[q][0]; 1701 xright = x[q][1]; 1702 ylow = y[q][0]; 1703 yhigh = y[q][1]; 1704 PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx)); 1705 nidx *= dof; 1706 n[q] = nidx; 1707 /* 1708 Based on the counted number of indices in the local domain *with an overlap*, 1709 construct a subcommunicator of all the MPI ranks supporting this domain. 1710 Observe that a domain with an overlap might have nontrivial local support, 1711 while the domain without an overlap might not. Hence, the decision to participate 1712 in the subcommunicator must be based on the domain with an overlap. 1713 */ 1714 if (q == 0) { 1715 if (nidx) color = 1; 1716 else color = MPI_UNDEFINED; 1717 PetscCallMPI(MPI_Comm_split(comm, color, rank, &subcomm)); 1718 split = PETSC_TRUE; 1719 } 1720 /* 1721 Proceed only if the number of local indices *with an overlap* is nonzero. 1722 */ 1723 if (n[0]) { 1724 if (q == 0) xis = is; 1725 if (q == 1) { 1726 /* 1727 The IS for the no-overlap subdomain shares a communicator with the overlapping domain. 1728 Moreover, if the overlap is zero, the two ISs are identical. 1729 */ 1730 if (overlap == 0) { 1731 (*is_local)[s] = (*is)[s]; 1732 PetscCall(PetscObjectReference((PetscObject)(*is)[s])); 1733 continue; 1734 } else { 1735 xis = is_local; 1736 subcomm = ((PetscObject)(*is)[s])->comm; 1737 } 1738 } /* if (q == 1) */ 1739 idx = NULL; 1740 PetscCall(PetscMalloc1(nidx, &idx)); 1741 if (nidx) { 1742 k = 0; 1743 for (jj = ylow_loc; jj < yhigh_loc; ++jj) { 1744 PetscInt x0 = (jj == ylow_loc) ? xleft_loc : xleft; 1745 PetscInt x1 = (jj == yhigh_loc - 1) ? xright_loc : xright; 1746 kk = dof * (M * jj + x0); 1747 for (ii = x0; ii < x1; ++ii) { 1748 for (d = 0; d < dof; ++d) idx[k++] = kk++; 1749 } 1750 } 1751 } 1752 PetscCall(ISCreateGeneral(subcomm, nidx, idx, PETSC_OWN_POINTER, (*xis) + s)); 1753 if (split) PetscCallMPI(MPI_Comm_free(&subcomm)); 1754 } /* if (n[0]) */ 1755 } /* for (q = 0; q < 2; ++q) */ 1756 if (n[0]) ++s; 1757 xstart += maxwidth; 1758 } /* for (i = 0; i < Mdomains; ++i) */ 1759 ystart += maxheight; 1760 } /* for (j = 0; j < Ndomains; ++j) */ 1761 PetscFunctionReturn(PETSC_SUCCESS); 1762 } 1763 1764 /*@C 1765 PCGASMGetSubdomains - Gets the subdomains supported on this MPI rank 1766 for the `PCGASM` additive Schwarz preconditioner. 1767 1768 Not Collective 1769 1770 Input Parameter: 1771 . pc - the preconditioner context 1772 1773 Output Parameters: 1774 + n - the number of subdomains for this MPI rank (default value = 1) 1775 . iis - the index sets that define the inner subdomains (without overlap) supported on this rank (can be `NULL`) 1776 - ois - the index sets that define the outer subdomains (with overlap) supported on this rank (can be `NULL`) 1777 1778 Level: advanced 1779 1780 Notes: 1781 The user is responsible for destroying the `IS`s and freeing the returned arrays, this can be done with 1782 `PCGASMDestroySubdomains()` 1783 1784 The `IS` numbering is in the parallel, global numbering of the vector. 1785 1786 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, `PCGASMCreateSubdomains2D()`, 1787 `PCGASMSetSubdomains()`, `PCGASMGetSubmatrices()`, `PCGASMDestroySubdomains()` 1788 @*/ 1789 PetscErrorCode PCGASMGetSubdomains(PC pc, PetscInt *n, IS *iis[], IS *ois[]) 1790 { 1791 PC_GASM *osm; 1792 PetscBool match; 1793 PetscInt i; 1794 1795 PetscFunctionBegin; 1796 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1797 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1798 PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Incorrect object type: expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1799 osm = (PC_GASM *)pc->data; 1800 if (n) *n = osm->n; 1801 if (iis) PetscCall(PetscMalloc1(osm->n, iis)); 1802 if (ois) PetscCall(PetscMalloc1(osm->n, ois)); 1803 if (iis || ois) { 1804 for (i = 0; i < osm->n; ++i) { 1805 if (iis) (*iis)[i] = osm->iis[i]; 1806 if (ois) (*ois)[i] = osm->ois[i]; 1807 } 1808 } 1809 PetscFunctionReturn(PETSC_SUCCESS); 1810 } 1811 1812 /*@C 1813 PCGASMGetSubmatrices - Gets the local submatrices (for this MPI rank 1814 only) for the `PCGASM` additive Schwarz preconditioner. 1815 1816 Not Collective 1817 1818 Input Parameter: 1819 . pc - the preconditioner context 1820 1821 Output Parameters: 1822 + n - the number of matrices for this MPI rank (default value = 1) 1823 - mat - the matrices 1824 1825 Level: advanced 1826 1827 Note: 1828 Matrices returned by this routine have the same communicators as the index sets (`IS`) 1829 used to define subdomains in `PCGASMSetSubdomains()` 1830 1831 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, 1832 `PCGASMCreateSubdomains2D()`, `PCGASMSetSubdomains()`, `PCGASMGetSubdomains()` 1833 @*/ 1834 PetscErrorCode PCGASMGetSubmatrices(PC pc, PetscInt *n, Mat *mat[]) 1835 { 1836 PC_GASM *osm; 1837 PetscBool match; 1838 1839 PetscFunctionBegin; 1840 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1841 PetscAssertPointer(n, 2); 1842 if (mat) PetscAssertPointer(mat, 3); 1843 PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call after KSPSetUp() or PCSetUp()."); 1844 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1845 PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name); 1846 osm = (PC_GASM *)pc->data; 1847 if (n) *n = osm->n; 1848 if (mat) *mat = osm->pmat; 1849 PetscFunctionReturn(PETSC_SUCCESS); 1850 } 1851 1852 /*@ 1853 PCGASMSetUseDMSubdomains - Indicates whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible for `PCGASM` 1854 1855 Logically Collective 1856 1857 Input Parameters: 1858 + pc - the preconditioner 1859 - flg - boolean indicating whether to use subdomains defined by the `DM` 1860 1861 Options Database Key: 1862 + -pc_gasm_dm_subdomains - configure subdomains 1863 . -pc_gasm_overlap - set overlap 1864 - -pc_gasm_total_subdomains - set number of subdomains 1865 1866 Level: intermediate 1867 1868 Note: 1869 `PCGASMSetSubdomains()`, `PCGASMSetTotalSubdomains()` or `PCGASMSetOverlap()` take precedence over `PCGASMSetUseDMSubdomains()`, 1870 so setting `PCGASMSetSubdomains()` with nontrivial subdomain ISs or any of `PCGASMSetTotalSubdomains()` and `PCGASMSetOverlap()` 1871 automatically turns the latter off. 1872 1873 .seealso: [](ch_ksp), `PCGASM`, `PCGASMGetUseDMSubdomains()`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()` 1874 `PCGASMCreateSubdomains2D()` 1875 @*/ 1876 PetscErrorCode PCGASMSetUseDMSubdomains(PC pc, PetscBool flg) 1877 { 1878 PC_GASM *osm = (PC_GASM *)pc->data; 1879 PetscBool match; 1880 1881 PetscFunctionBegin; 1882 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1883 PetscValidLogicalCollectiveBool(pc, flg, 2); 1884 PetscCheck(!pc->setupcalled, ((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Not for a setup PC."); 1885 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1886 if (match) { 1887 if (!osm->user_subdomains && osm->N == PETSC_DETERMINE && osm->overlap < 0) osm->dm_subdomains = flg; 1888 } 1889 PetscFunctionReturn(PETSC_SUCCESS); 1890 } 1891 1892 /*@ 1893 PCGASMGetUseDMSubdomains - Returns flag indicating whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible with `PCGASM` 1894 1895 Not Collective 1896 1897 Input Parameter: 1898 . pc - the preconditioner 1899 1900 Output Parameter: 1901 . flg - boolean indicating whether to use subdomains defined by the `DM` 1902 1903 Level: intermediate 1904 1905 .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetUseDMSubdomains()`, `PCGASMSetOverlap()` 1906 `PCGASMCreateSubdomains2D()` 1907 @*/ 1908 PetscErrorCode PCGASMGetUseDMSubdomains(PC pc, PetscBool *flg) 1909 { 1910 PC_GASM *osm = (PC_GASM *)pc->data; 1911 PetscBool match; 1912 1913 PetscFunctionBegin; 1914 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1915 PetscAssertPointer(flg, 2); 1916 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match)); 1917 if (match) { 1918 if (flg) *flg = osm->dm_subdomains; 1919 } 1920 PetscFunctionReturn(PETSC_SUCCESS); 1921 } 1922