1 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/ 2 #include <petscds.h> 3 4 // Greatest common divisor of two nonnegative integers 5 PetscInt PetscGCD(PetscInt a, PetscInt b) { 6 while (b != 0) { 7 PetscInt tmp = b; 8 b = a % b; 9 a = tmp; 10 } 11 return a; 12 } 13 14 PetscErrorCode DMCreateGlobalVector_Section_Private(DM dm,Vec *vec) 15 { 16 PetscSection gSection; 17 PetscInt localSize, bs, blockSize = -1, pStart, pEnd, p; 18 PetscErrorCode ierr; 19 PetscInt in[2],out[2]; 20 21 PetscFunctionBegin; 22 ierr = DMGetGlobalSection(dm, &gSection);CHKERRQ(ierr); 23 ierr = PetscSectionGetChart(gSection, &pStart, &pEnd);CHKERRQ(ierr); 24 for (p = pStart; p < pEnd; ++p) { 25 PetscInt dof, cdof; 26 27 ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr); 28 ierr = PetscSectionGetConstraintDof(gSection, p, &cdof);CHKERRQ(ierr); 29 30 if (dof - cdof > 0) { 31 if (blockSize < 0) { 32 /* set blockSize */ 33 blockSize = dof-cdof; 34 } else { 35 blockSize = PetscGCD(dof - cdof, blockSize); 36 } 37 } 38 } 39 40 in[0] = blockSize < 0 ? PETSC_MIN_INT : -blockSize; 41 in[1] = blockSize; 42 ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));CHKERRMPI(ierr); 43 /* -out[0] = min(blockSize), out[1] = max(blockSize) */ 44 if (-out[0] == out[1]) { 45 bs = out[1]; 46 } else bs = 1; 47 48 if (bs < 0) { /* Everyone was empty */ 49 blockSize = 1; 50 bs = 1; 51 } 52 53 ierr = PetscSectionGetConstrainedStorageSize(gSection, &localSize);CHKERRQ(ierr); 54 PetscCheckFalse(localSize%blockSize,PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %d and local storage size %d", blockSize, localSize); 55 ierr = VecCreate(PetscObjectComm((PetscObject)dm), vec);CHKERRQ(ierr); 56 ierr = VecSetSizes(*vec, localSize, PETSC_DETERMINE);CHKERRQ(ierr); 57 ierr = VecSetBlockSize(*vec, bs);CHKERRQ(ierr); 58 ierr = VecSetType(*vec,dm->vectype);CHKERRQ(ierr); 59 ierr = VecSetDM(*vec, dm);CHKERRQ(ierr); 60 /* ierr = VecSetLocalToGlobalMapping(*vec, dm->ltogmap);CHKERRQ(ierr); */ 61 PetscFunctionReturn(0); 62 } 63 64 PetscErrorCode DMCreateLocalVector_Section_Private(DM dm,Vec *vec) 65 { 66 PetscSection section; 67 PetscInt localSize, blockSize = -1, pStart, pEnd, p; 68 PetscErrorCode ierr; 69 70 PetscFunctionBegin; 71 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 72 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 73 for (p = pStart; p < pEnd; ++p) { 74 PetscInt dof; 75 76 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 77 if ((blockSize < 0) && (dof > 0)) blockSize = dof; 78 if (dof > 0) blockSize = PetscGCD(dof, blockSize); 79 } 80 ierr = PetscSectionGetStorageSize(section, &localSize);CHKERRQ(ierr); 81 ierr = VecCreate(PETSC_COMM_SELF, vec);CHKERRQ(ierr); 82 ierr = VecSetSizes(*vec, localSize, localSize);CHKERRQ(ierr); 83 ierr = VecSetBlockSize(*vec, blockSize);CHKERRQ(ierr); 84 ierr = VecSetType(*vec,dm->vectype);CHKERRQ(ierr); 85 ierr = VecSetDM(*vec, dm);CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 89 /*@C 90 DMCreateSectionSubDM - Returns an IS and subDM+subSection encapsulating a subproblem defined by the fields in a PetscSection in the DM. 91 92 Not collective 93 94 Input Parameters: 95 + dm - The DM object 96 . numFields - The number of fields in this subproblem 97 - fields - The field numbers of the selected fields 98 99 Output Parameters: 100 + is - The global indices for the subproblem 101 - subdm - The DM for the subproblem, which must already have be cloned from dm 102 103 Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes, 104 such as Plex and Forest. 105 106 Level: intermediate 107 108 .seealso DMCreateSubDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView() 109 @*/ 110 PetscErrorCode DMCreateSectionSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm) 111 { 112 PetscSection section, sectionGlobal; 113 PetscInt *subIndices; 114 PetscInt subSize = 0, subOff = 0, Nf, f, pStart, pEnd, p; 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 if (!numFields) PetscFunctionReturn(0); 119 ierr = DMGetLocalSection(dm, §ion);CHKERRQ(ierr); 120 ierr = DMGetGlobalSection(dm, §ionGlobal);CHKERRQ(ierr); 121 PetscCheckFalse(!section,PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 122 PetscCheckFalse(!sectionGlobal,PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 123 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 124 PetscCheckFalse(numFields > Nf,PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of DM fields %d", numFields, Nf); 125 if (is) { 126 PetscInt bs, bsLocal[2], bsMinMax[2]; 127 128 for (f = 0, bs = 0; f < numFields; ++f) { 129 PetscInt Nc; 130 131 ierr = PetscSectionGetFieldComponents(section, fields[f], &Nc);CHKERRQ(ierr); 132 bs += Nc; 133 } 134 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 135 for (p = pStart; p < pEnd; ++p) { 136 PetscInt gdof, pSubSize = 0; 137 138 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 139 if (gdof > 0) { 140 for (f = 0; f < numFields; ++f) { 141 PetscInt fdof, fcdof; 142 143 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr); 144 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr); 145 pSubSize += fdof-fcdof; 146 } 147 subSize += pSubSize; 148 if (pSubSize && bs != pSubSize) { 149 /* Layout does not admit a pointwise block size */ 150 bs = 1; 151 } 152 } 153 } 154 /* Must have same blocksize on all procs (some might have no points) */ 155 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs; 156 ierr = PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);CHKERRQ(ierr); 157 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;} 158 else {bs = bsMinMax[0];} 159 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr); 160 for (p = pStart; p < pEnd; ++p) { 161 PetscInt gdof, goff; 162 163 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 164 if (gdof > 0) { 165 ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr); 166 for (f = 0; f < numFields; ++f) { 167 PetscInt fdof, fcdof, fc, f2, poff = 0; 168 169 /* Can get rid of this loop by storing field information in the global section */ 170 for (f2 = 0; f2 < fields[f]; ++f2) { 171 ierr = PetscSectionGetFieldDof(section, p, f2, &fdof);CHKERRQ(ierr); 172 ierr = PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);CHKERRQ(ierr); 173 poff += fdof-fcdof; 174 } 175 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr); 176 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr); 177 for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) { 178 subIndices[subOff] = goff+poff+fc; 179 } 180 } 181 } 182 } 183 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);CHKERRQ(ierr); 184 if (bs > 1) { 185 /* We need to check that the block size does not come from non-contiguous fields */ 186 PetscInt i, j, set = 1; 187 for (i = 0; i < subSize; i += bs) { 188 for (j = 0; j < bs; ++j) { 189 if (subIndices[i+j] != subIndices[i]+j) {set = 0; break;} 190 } 191 } 192 if (set) {ierr = ISSetBlockSize(*is, bs);CHKERRQ(ierr);} 193 } 194 } 195 if (subdm) { 196 PetscSection subsection; 197 PetscBool haveNull = PETSC_FALSE; 198 PetscInt f, nf = 0, of = 0; 199 200 ierr = PetscSectionCreateSubsection(section, numFields, fields, &subsection);CHKERRQ(ierr); 201 ierr = DMSetLocalSection(*subdm, subsection);CHKERRQ(ierr); 202 ierr = PetscSectionDestroy(&subsection);CHKERRQ(ierr); 203 for (f = 0; f < numFields; ++f) { 204 (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]]; 205 if ((*subdm)->nullspaceConstructors[f]) { 206 haveNull = PETSC_TRUE; 207 nf = f; 208 of = fields[f]; 209 } 210 } 211 if (dm->probs) { 212 ierr = DMSetNumFields(*subdm, numFields);CHKERRQ(ierr); 213 for (f = 0; f < numFields; ++f) { 214 PetscObject disc; 215 216 ierr = DMGetField(dm, fields[f], NULL, &disc);CHKERRQ(ierr); 217 ierr = DMSetField(*subdm, f, NULL, disc);CHKERRQ(ierr); 218 } 219 ierr = DMCreateDS(*subdm);CHKERRQ(ierr); 220 if (numFields == 1 && is) { 221 PetscObject disc, space, pmat; 222 223 ierr = DMGetField(*subdm, 0, NULL, &disc);CHKERRQ(ierr); 224 ierr = PetscObjectQuery(disc, "nullspace", &space);CHKERRQ(ierr); 225 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nullspace", space);CHKERRQ(ierr);} 226 ierr = PetscObjectQuery(disc, "nearnullspace", &space);CHKERRQ(ierr); 227 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nearnullspace", space);CHKERRQ(ierr);} 228 ierr = PetscObjectQuery(disc, "pmat", &pmat);CHKERRQ(ierr); 229 if (pmat) {ierr = PetscObjectCompose((PetscObject) *is, "pmat", pmat);CHKERRQ(ierr);} 230 } 231 /* Check if DSes record their DM fields */ 232 if (dm->probs[0].fields) { 233 PetscInt d, e; 234 235 for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) { 236 const PetscInt Nf = dm->probs[d].ds->Nf; 237 const PetscInt *fld; 238 PetscInt f, g; 239 240 ierr = ISGetIndices(dm->probs[d].fields, &fld);CHKERRQ(ierr); 241 for (f = 0; f < Nf; ++f) { 242 for (g = 0; g < numFields; ++g) if (fld[f] == fields[g]) break; 243 if (g < numFields) break; 244 } 245 ierr = ISRestoreIndices(dm->probs[d].fields, &fld);CHKERRQ(ierr); 246 if (f == Nf) continue; 247 ierr = PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds);CHKERRQ(ierr); 248 ierr = PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds);CHKERRQ(ierr); 249 /* Translate DM fields to DS fields */ 250 { 251 IS infields, dsfields; 252 const PetscInt *fld, *ofld; 253 PetscInt *fidx; 254 PetscInt onf, nf, f, g; 255 256 ierr = ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields);CHKERRQ(ierr); 257 ierr = ISIntersect(infields, dm->probs[d].fields, &dsfields);CHKERRQ(ierr); 258 ierr = ISDestroy(&infields);CHKERRQ(ierr); 259 ierr = ISGetLocalSize(dsfields, &nf);CHKERRQ(ierr); 260 PetscCheckFalse(!nf,PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields"); 261 ierr = ISGetIndices(dsfields, &fld);CHKERRQ(ierr); 262 ierr = ISGetLocalSize(dm->probs[d].fields, &onf);CHKERRQ(ierr); 263 ierr = ISGetIndices(dm->probs[d].fields, &ofld);CHKERRQ(ierr); 264 ierr = PetscMalloc1(nf, &fidx);CHKERRQ(ierr); 265 for (f = 0, g = 0; f < onf && g < nf; ++f) { 266 if (ofld[f] == fld[g]) fidx[g++] = f; 267 } 268 ierr = ISRestoreIndices(dm->probs[d].fields, &ofld);CHKERRQ(ierr); 269 ierr = ISRestoreIndices(dsfields, &fld);CHKERRQ(ierr); 270 ierr = ISDestroy(&dsfields);CHKERRQ(ierr); 271 ierr = PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);CHKERRQ(ierr); 272 ierr = PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds);CHKERRQ(ierr); 273 ierr = PetscFree(fidx);CHKERRQ(ierr); 274 } 275 ++e; 276 } 277 } else { 278 ierr = PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds);CHKERRQ(ierr); 279 ierr = PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds);CHKERRQ(ierr); 280 ierr = PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);CHKERRQ(ierr); 281 ierr = PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds);CHKERRQ(ierr); 282 } 283 } 284 if (haveNull && is) { 285 MatNullSpace nullSpace; 286 287 ierr = (*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace);CHKERRQ(ierr); 288 ierr = PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr); 289 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr); 290 } 291 if (dm->coarseMesh) { 292 ierr = DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh);CHKERRQ(ierr); 293 } 294 } 295 PetscFunctionReturn(0); 296 } 297 298 /*@C 299 DMCreateSectionSuperDM - Returns an arrays of ISes and DM+Section encapsulating a superproblem defined by the DM+Sections passed in. 300 301 Not collective 302 303 Input Parameters: 304 + dms - The DM objects 305 - len - The number of DMs 306 307 Output Parameters: 308 + is - The global indices for the subproblem, or NULL 309 - superdm - The DM for the superproblem, which must already have be cloned 310 311 Note: This handles all information in the DM class and the PetscSection. This is used as the basis for creating subDMs in specialized classes, 312 such as Plex and Forest. 313 314 Level: intermediate 315 316 .seealso DMCreateSuperDM(), DMGetLocalSection(), DMPlexSetMigrationSF(), DMView() 317 @*/ 318 PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm) 319 { 320 MPI_Comm comm; 321 PetscSection supersection, *sections, *sectionGlobals; 322 PetscInt *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i; 323 PetscBool haveNull = PETSC_FALSE; 324 PetscErrorCode ierr; 325 326 PetscFunctionBegin; 327 ierr = PetscObjectGetComm((PetscObject)dms[0], &comm);CHKERRQ(ierr); 328 /* Pull out local and global sections */ 329 ierr = PetscMalloc3(len, &Nfs, len, §ions, len, §ionGlobals);CHKERRQ(ierr); 330 for (i = 0 ; i < len; ++i) { 331 ierr = DMGetLocalSection(dms[i], §ions[i]);CHKERRQ(ierr); 332 ierr = DMGetGlobalSection(dms[i], §ionGlobals[i]);CHKERRQ(ierr); 333 PetscCheckFalse(!sections[i],comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields"); 334 PetscCheckFalse(!sectionGlobals[i],comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields"); 335 ierr = PetscSectionGetNumFields(sections[i], &Nfs[i]);CHKERRQ(ierr); 336 Nf += Nfs[i]; 337 } 338 /* Create the supersection */ 339 ierr = PetscSectionCreateSupersection(sections, len, &supersection);CHKERRQ(ierr); 340 ierr = DMSetLocalSection(*superdm, supersection);CHKERRQ(ierr); 341 /* Create ISes */ 342 if (is) { 343 PetscSection supersectionGlobal; 344 PetscInt bs = -1, startf = 0; 345 346 ierr = PetscMalloc1(len, is);CHKERRQ(ierr); 347 ierr = DMGetGlobalSection(*superdm, &supersectionGlobal);CHKERRQ(ierr); 348 for (i = 0 ; i < len; startf += Nfs[i], ++i) { 349 PetscInt *subIndices; 350 PetscInt subSize, subOff, pStart, pEnd, p, start, end, dummy; 351 352 ierr = PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd);CHKERRQ(ierr); 353 ierr = PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize);CHKERRQ(ierr); 354 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr); 355 for (p = pStart, subOff = 0; p < pEnd; ++p) { 356 PetscInt gdof, gcdof, gtdof, d; 357 358 ierr = PetscSectionGetDof(sectionGlobals[i], p, &gdof);CHKERRQ(ierr); 359 ierr = PetscSectionGetConstraintDof(sections[i], p, &gcdof);CHKERRQ(ierr); 360 gtdof = gdof - gcdof; 361 if (gdof > 0 && gtdof) { 362 if (bs < 0) {bs = gtdof;} 363 else if (bs != gtdof) {bs = 1;} 364 ierr = DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy);CHKERRQ(ierr); 365 ierr = DMGetGlobalFieldOffset_Private(*superdm, p, startf+Nfs[i]-1, &dummy, &end);CHKERRQ(ierr); 366 PetscCheckFalse(end-start != gtdof,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid number of global dofs %D != %D for dm %D on point %D", end-start, gtdof, i, p); 367 for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d; 368 } 369 } 370 ierr = ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]);CHKERRQ(ierr); 371 /* Must have same blocksize on all procs (some might have no points) */ 372 { 373 PetscInt bs = -1, bsLocal[2], bsMinMax[2]; 374 375 bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs; 376 ierr = PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax);CHKERRQ(ierr); 377 if (bsMinMax[0] != bsMinMax[1]) {bs = 1;} 378 else {bs = bsMinMax[0];} 379 ierr = ISSetBlockSize((*is)[i], bs);CHKERRQ(ierr); 380 } 381 } 382 } 383 /* Preserve discretizations */ 384 if (len && dms[0]->probs) { 385 ierr = DMSetNumFields(*superdm, Nf);CHKERRQ(ierr); 386 for (i = 0, supf = 0; i < len; ++i) { 387 for (f = 0; f < Nfs[i]; ++f, ++supf) { 388 PetscObject disc; 389 390 ierr = DMGetField(dms[i], f, NULL, &disc);CHKERRQ(ierr); 391 ierr = DMSetField(*superdm, supf, NULL, disc);CHKERRQ(ierr); 392 } 393 } 394 ierr = DMCreateDS(*superdm);CHKERRQ(ierr); 395 } 396 /* Preserve nullspaces */ 397 for (i = 0, supf = 0; i < len; ++i) { 398 for (f = 0; f < Nfs[i]; ++f, ++supf) { 399 (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f]; 400 if ((*superdm)->nullspaceConstructors[supf]) { 401 haveNull = PETSC_TRUE; 402 nullf = supf; 403 oldf = f; 404 } 405 } 406 } 407 /* Attach nullspace to IS */ 408 if (haveNull && is) { 409 MatNullSpace nullSpace; 410 411 ierr = (*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace);CHKERRQ(ierr); 412 ierr = PetscObjectCompose((PetscObject) (*is)[nullf], "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr); 413 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr); 414 } 415 ierr = PetscSectionDestroy(&supersection);CHKERRQ(ierr); 416 ierr = PetscFree3(Nfs, sections, sectionGlobals);CHKERRQ(ierr); 417 PetscFunctionReturn(0); 418 } 419