1 #define PETSCDM_DLL 2 #include <petsc/private/dmnetworkimpl.h> /*I "petscdmnetwork.h" I*/ 3 #include <petsc/private/vecimpl.h> 4 5 PetscErrorCode DMSetFromOptions_Network(DM dm, PetscOptionItems *PetscOptionsObject) 6 { 7 PetscFunctionBegin; 8 PetscOptionsHeadBegin(PetscOptionsObject, "DMNetwork Options"); 9 PetscOptionsHeadEnd(); 10 PetscFunctionReturn(0); 11 } 12 13 /* External function declarations here */ 14 extern PetscErrorCode DMCreateMatrix_Network(DM, Mat *); 15 extern PetscErrorCode DMDestroy_Network(DM); 16 extern PetscErrorCode DMView_Network(DM, PetscViewer); 17 extern PetscErrorCode DMGlobalToLocalBegin_Network(DM, Vec, InsertMode, Vec); 18 extern PetscErrorCode DMGlobalToLocalEnd_Network(DM, Vec, InsertMode, Vec); 19 extern PetscErrorCode DMLocalToGlobalBegin_Network(DM, Vec, InsertMode, Vec); 20 extern PetscErrorCode DMLocalToGlobalEnd_Network(DM, Vec, InsertMode, Vec); 21 extern PetscErrorCode DMSetUp_Network(DM); 22 extern PetscErrorCode DMClone_Network(DM, DM *); 23 24 static PetscErrorCode VecArrayPrint_private(PetscViewer viewer, PetscInt n, const PetscScalar *xv) 25 { 26 PetscInt i; 27 28 PetscFunctionBegin; 29 for (i = 0; i < n; i++) { 30 #if defined(PETSC_USE_COMPLEX) 31 if (PetscImaginaryPart(xv[i]) > 0.0) { 32 PetscCall(PetscViewerASCIIPrintf(viewer, " %g + %g i\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i]))); 33 } else if (PetscImaginaryPart(xv[i]) < 0.0) { 34 PetscCall(PetscViewerASCIIPrintf(viewer, " %g - %g i\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i]))); 35 } else PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)PetscRealPart(xv[i]))); 36 #else 37 PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)xv[i])); 38 #endif 39 } 40 PetscFunctionReturn(0); 41 } 42 43 static PetscErrorCode VecView_Network_Seq(DM networkdm, Vec X, PetscViewer viewer) 44 { 45 PetscInt e, v, Start, End, offset, nvar, id; 46 const PetscScalar *xv; 47 48 PetscFunctionBegin; 49 PetscCall(VecGetArrayRead(X, &xv)); 50 51 /* iterate over edges */ 52 PetscCall(DMNetworkGetEdgeRange(networkdm, &Start, &End)); 53 for (e = Start; e < End; e++) { 54 PetscCall(DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar)); 55 if (!nvar) continue; 56 57 PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset)); 58 PetscCall(DMNetworkGetGlobalEdgeIndex(networkdm, e, &id)); 59 60 PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id)); 61 PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset)); 62 } 63 64 /* iterate over vertices */ 65 PetscCall(DMNetworkGetVertexRange(networkdm, &Start, &End)); 66 for (v = Start; v < End; v++) { 67 PetscCall(DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar)); 68 if (!nvar) continue; 69 70 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset)); 71 PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, v, &id)); 72 73 PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id)); 74 PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset)); 75 } 76 PetscCall(PetscViewerFlush(viewer)); 77 PetscCall(VecRestoreArrayRead(X, &xv)); 78 PetscFunctionReturn(0); 79 } 80 81 static PetscErrorCode VecView_Network_MPI(DM networkdm, Vec X, PetscViewer viewer) 82 { 83 PetscInt i, e, v, eStart, eEnd, vStart, vEnd, offset, nvar, len_loc, len, k; 84 const PetscScalar *xv; 85 MPI_Comm comm; 86 PetscMPIInt size, rank, tag = ((PetscObject)viewer)->tag; 87 Vec localX; 88 PetscBool ghostvtex; 89 PetscScalar *values; 90 PetscInt j, ne, nv, id; 91 MPI_Status status; 92 93 PetscFunctionBegin; 94 PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm)); 95 PetscCallMPI(MPI_Comm_size(comm, &size)); 96 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 97 98 PetscCall(DMGetLocalVector(networkdm, &localX)); 99 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX)); 100 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX)); 101 PetscCall(VecGetArrayRead(localX, &xv)); 102 103 PetscCall(VecGetLocalSize(localX, &len_loc)); 104 105 PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd)); 106 PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd)); 107 len_loc += 2 * (1 + eEnd - eStart + vEnd - vStart); 108 109 /* values = [nedges, nvertices; id, nvar, xedge; ...; id, nvars, xvertex;...], to be sent to proc[0] */ 110 PetscCallMPI(MPI_Allreduce(&len_loc, &len, 1, MPIU_INT, MPI_MAX, comm)); 111 PetscCall(PetscCalloc1(len, &values)); 112 113 if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", rank)); 114 115 /* iterate over edges */ 116 k = 2; 117 for (e = eStart; e < eEnd; e++) { 118 PetscCall(DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar)); 119 if (!nvar) continue; 120 121 PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset)); 122 PetscCall(DMNetworkGetGlobalEdgeIndex(networkdm, e, &id)); 123 124 if (rank == 0) { /* print its own entries */ 125 PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id)); 126 PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset)); 127 } else { 128 values[0] += 1; /* number of edges */ 129 values[k++] = id; 130 values[k++] = nvar; 131 for (i = offset; i < offset + nvar; i++) values[k++] = xv[i]; 132 } 133 } 134 135 /* iterate over vertices */ 136 for (v = vStart; v < vEnd; v++) { 137 PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex)); 138 if (ghostvtex) continue; 139 PetscCall(DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar)); 140 if (!nvar) continue; 141 142 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset)); 143 PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, v, &id)); 144 145 if (rank == 0) { 146 PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id)); 147 PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset)); 148 } else { 149 values[1] += 1; /* number of vertices */ 150 values[k++] = id; 151 values[k++] = nvar; 152 for (i = offset; i < offset + nvar; i++) values[k++] = xv[i]; 153 } 154 } 155 156 if (rank == 0) { 157 /* proc[0] receives and prints messages */ 158 for (j = 1; j < size; j++) { 159 PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%" PetscInt_FMT "]\n", j)); 160 161 PetscCallMPI(MPI_Recv(values, (PetscMPIInt)len, MPIU_SCALAR, j, tag, comm, &status)); 162 163 ne = (PetscInt)PetscAbsScalar(values[0]); 164 nv = (PetscInt)PetscAbsScalar(values[1]); 165 166 /* print received edges */ 167 k = 2; 168 for (i = 0; i < ne; i++) { 169 id = (PetscInt)PetscAbsScalar(values[k++]); 170 nvar = (PetscInt)PetscAbsScalar(values[k++]); 171 PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id)); 172 PetscCall(VecArrayPrint_private(viewer, nvar, values + k)); 173 k += nvar; 174 } 175 176 /* print received vertices */ 177 for (i = 0; i < nv; i++) { 178 id = (PetscInt)PetscAbsScalar(values[k++]); 179 nvar = (PetscInt)PetscAbsScalar(values[k++]); 180 PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id)); 181 PetscCall(VecArrayPrint_private(viewer, nvar, values + k)); 182 k += nvar; 183 } 184 } 185 } else { 186 /* sends values to proc[0] */ 187 PetscCallMPI(MPI_Send((void *)values, k, MPIU_SCALAR, 0, tag, comm)); 188 } 189 190 PetscCall(PetscFree(values)); 191 PetscCall(VecRestoreArrayRead(localX, &xv)); 192 PetscCall(DMRestoreLocalVector(networkdm, &localX)); 193 PetscFunctionReturn(0); 194 } 195 196 PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer); 197 198 PetscErrorCode VecView_Network(Vec v, PetscViewer viewer) 199 { 200 DM dm; 201 PetscBool isseq; 202 PetscBool iascii; 203 204 PetscFunctionBegin; 205 PetscCall(VecGetDM(v, &dm)); 206 PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM"); 207 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 208 PetscCall(PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq)); 209 210 /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */ 211 if (iascii) { 212 if (isseq) PetscCall(VecView_Network_Seq(dm, v, viewer)); 213 else PetscCall(VecView_Network_MPI(dm, v, viewer)); 214 } else { 215 if (isseq) PetscCall(VecView_Seq(v, viewer)); 216 else PetscCall(VecView_MPI(v, viewer)); 217 } 218 PetscFunctionReturn(0); 219 } 220 221 static PetscErrorCode DMCreateGlobalVector_Network(DM dm, Vec *vec) 222 { 223 DM_Network *network = (DM_Network *)dm->data; 224 225 PetscFunctionBegin; 226 PetscCall(DMCreateGlobalVector(network->plex, vec)); 227 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Network)); 228 PetscCall(VecSetDM(*vec, dm)); 229 PetscFunctionReturn(0); 230 } 231 232 static PetscErrorCode DMCreateLocalVector_Network(DM dm, Vec *vec) 233 { 234 DM_Network *network = (DM_Network *)dm->data; 235 236 PetscFunctionBegin; 237 PetscCall(DMCreateLocalVector(network->plex, vec)); 238 PetscCall(VecSetDM(*vec, dm)); 239 PetscFunctionReturn(0); 240 } 241 242 PetscErrorCode DMNetworkInitializeToDefault_NonShared(DM dm) 243 { 244 DM_Network *network = (DM_Network *)dm->data; 245 246 PetscFunctionBegin; 247 network->Je = NULL; 248 network->Jv = NULL; 249 network->Jvptr = NULL; 250 network->userEdgeJacobian = PETSC_FALSE; 251 network->userVertexJacobian = PETSC_FALSE; 252 253 network->vertex.DofSection = NULL; 254 network->vertex.GlobalDofSection = NULL; 255 network->vertex.mapping = NULL; 256 network->vertex.sf = NULL; 257 258 network->edge.DofSection = NULL; 259 network->edge.GlobalDofSection = NULL; 260 network->edge.mapping = NULL; 261 network->edge.sf = NULL; 262 263 network->DataSection = NULL; 264 network->DofSection = NULL; 265 network->GlobalDofSection = NULL; 266 network->componentsetup = PETSC_FALSE; 267 268 network->plex = NULL; 269 270 network->component = NULL; 271 network->ncomponent = 0; 272 273 network->header = NULL; 274 network->cvalue = NULL; 275 network->componentdataarray = NULL; 276 277 network->max_comps_registered = DMNETWORK_MAX_COMP_REGISTERED_DEFAULT; /* return to default */ 278 PetscFunctionReturn(0); 279 } 280 /* Default values for the parameters in DMNetwork */ 281 PetscErrorCode DMNetworkInitializeToDefault(DM dm) 282 { 283 DM_Network *network = (DM_Network *)dm->data; 284 DMNetworkCloneShared cloneshared = network->cloneshared; 285 286 PetscFunctionBegin; 287 PetscCall(DMNetworkInitializeToDefault_NonShared(dm)); 288 /* Default values for shared data */ 289 cloneshared->refct = 1; 290 cloneshared->NVertices = 0; 291 cloneshared->NEdges = 0; 292 cloneshared->nVertices = 0; 293 cloneshared->nEdges = 0; 294 cloneshared->nsubnet = 0; 295 cloneshared->pStart = -1; 296 cloneshared->pEnd = -1; 297 cloneshared->vStart = -1; 298 cloneshared->vEnd = -1; 299 cloneshared->eStart = -1; 300 cloneshared->eEnd = -1; 301 cloneshared->vltog = NULL; 302 cloneshared->distributecalled = PETSC_FALSE; 303 304 cloneshared->subnet = NULL; 305 cloneshared->subnetvtx = NULL; 306 cloneshared->subnetedge = NULL; 307 cloneshared->svtx = NULL; 308 cloneshared->nsvtx = 0; 309 cloneshared->Nsvtx = 0; 310 cloneshared->svertices = NULL; 311 cloneshared->sedgelist = NULL; 312 cloneshared->svtable = NULL; 313 PetscFunctionReturn(0); 314 } 315 316 PetscErrorCode DMInitialize_Network(DM dm) 317 { 318 PetscFunctionBegin; 319 PetscCall(DMSetDimension(dm, 1)); 320 dm->ops->view = DMView_Network; 321 dm->ops->setfromoptions = DMSetFromOptions_Network; 322 dm->ops->clone = DMClone_Network; 323 dm->ops->setup = DMSetUp_Network; 324 dm->ops->createglobalvector = DMCreateGlobalVector_Network; 325 dm->ops->createlocalvector = DMCreateLocalVector_Network; 326 dm->ops->getlocaltoglobalmapping = NULL; 327 dm->ops->createfieldis = NULL; 328 dm->ops->createcoordinatedm = NULL; 329 dm->ops->getcoloring = NULL; 330 dm->ops->creatematrix = DMCreateMatrix_Network; 331 dm->ops->createinterpolation = NULL; 332 dm->ops->createinjection = NULL; 333 dm->ops->refine = NULL; 334 dm->ops->coarsen = NULL; 335 dm->ops->refinehierarchy = NULL; 336 dm->ops->coarsenhierarchy = NULL; 337 dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network; 338 dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network; 339 dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network; 340 dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network; 341 dm->ops->destroy = DMDestroy_Network; 342 dm->ops->createsubdm = NULL; 343 dm->ops->locatepoints = NULL; 344 PetscFunctionReturn(0); 345 } 346 /* 347 copies over the subnetid and index portions of the DMNetworkComponentHeader from orignal dm to the newdm 348 */ 349 static PetscErrorCode DMNetworkCopyHeaderTopological(DM dm, DM newdm) 350 { 351 DM_Network *network = (DM_Network *)dm->data, *newnetwork = (DM_Network *)newdm->data; 352 PetscInt p, i, np, index, subnetid; 353 354 PetscFunctionBegin; 355 np = network->cloneshared->pEnd - network->cloneshared->pStart; 356 PetscCall(PetscCalloc2(np, &newnetwork->header, np, &newnetwork->cvalue)); 357 for (i = 0; i < np; i++) { 358 p = i + network->cloneshared->pStart; 359 PetscCall(DMNetworkGetSubnetID(dm, p, &subnetid)); 360 PetscCall(DMNetworkGetIndex(dm, p, &index)); 361 newnetwork->header[i].index = index; 362 newnetwork->header[i].subnetid = subnetid; 363 newnetwork->header[i].size = NULL; 364 newnetwork->header[i].key = NULL; 365 newnetwork->header[i].offset = NULL; 366 newnetwork->header[i].nvar = NULL; 367 newnetwork->header[i].offsetvarrel = NULL; 368 newnetwork->header[i].ndata = 0; 369 newnetwork->header[i].maxcomps = DMNETWORK_MAX_COMP_AT_POINT_DEFAULT; 370 newnetwork->header[i].hsize = sizeof(struct _p_DMNetworkComponentHeader) / sizeof(sizeof(DMNetworkComponentGenericDataType)); 371 } 372 PetscFunctionReturn(0); 373 } 374 375 PetscErrorCode DMClone_Network(DM dm, DM *newdm) 376 { 377 DM_Network *network = (DM_Network *)dm->data, *newnetwork = NULL; 378 379 PetscFunctionBegin; 380 network->cloneshared->refct++; 381 PetscCall(PetscNew(&newnetwork)); 382 (*newdm)->data = newnetwork; 383 PetscCall(DMNetworkInitializeToDefault_NonShared(*newdm)); 384 newnetwork->cloneshared = network->cloneshared; /* Share all data that can be cloneshared */ 385 PetscCall(DMClone(network->plex, &newnetwork->plex)); 386 PetscCall(DMNetworkCopyHeaderTopological(dm, *newdm)); 387 PetscCall(DMNetworkInitializeNonTopological(*newdm)); /* initialize all non-topological data to the state after DMNetworkLayoutSetUp as been called */ 388 PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMNETWORK)); 389 PetscCall(DMInitialize_Network(*newdm)); 390 PetscFunctionReturn(0); 391 } 392 393 /*MC 394 DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object 395 DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of 396 the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary. 397 This is specified by a PetscSection object. Ownership in the global representation is determined by 398 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 399 400 Level: intermediate 401 402 .seealso: `DMType`, `DMNetworkCreate()`, `DMCreate()`, `DMSetType()` 403 M*/ 404 405 PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm) 406 { 407 DM_Network *network; 408 409 PetscFunctionBegin; 410 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 411 PetscCall(PetscNew(&network)); 412 PetscCall(PetscNew(&network->cloneshared)); 413 dm->data = network; 414 415 PetscCall(DMNetworkInitializeToDefault(dm)); 416 PetscCall(DMInitialize_Network(dm)); 417 PetscFunctionReturn(0); 418 } 419 420 /*@ 421 DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network. 422 423 Collective 424 425 Input Parameter: 426 . comm - The communicator for the DMNetwork object 427 428 Output Parameter: 429 . network - The DMNetwork object 430 431 Level: beginner 432 433 @*/ 434 PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network) 435 { 436 PetscFunctionBegin; 437 PetscValidPointer(network, 2); 438 PetscCall(DMCreate(comm, network)); 439 PetscCall(DMSetType(*network, DMNETWORK)); 440 PetscFunctionReturn(0); 441 } 442