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) { 114 PetscCall(PetscViewerASCIIPrintf(viewer,"Process [%d]\n",rank)); 115 } 116 117 /* iterate over edges */ 118 k = 2; 119 for (e=eStart; e<eEnd; e++) { 120 PetscCall(DMNetworkGetComponent(networkdm,e,ALL_COMPONENTS,NULL,NULL,&nvar)); 121 if (!nvar) continue; 122 123 PetscCall(DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&offset)); 124 PetscCall(DMNetworkGetGlobalEdgeIndex(networkdm,e,&id)); 125 126 if (rank == 0) { /* print its own entries */ 127 PetscCall(PetscViewerASCIIPrintf(viewer," Edge %" PetscInt_FMT ":\n",id)); 128 PetscCall(VecArrayPrint_private(viewer,nvar,xv+offset)); 129 } else { 130 values[0] += 1; /* number of edges */ 131 values[k++] = id; 132 values[k++] = nvar; 133 for (i=offset; i< offset+nvar; i++) values[k++] = xv[i]; 134 } 135 } 136 137 /* iterate over vertices */ 138 for (v=vStart; v<vEnd; v++) { 139 PetscCall(DMNetworkIsGhostVertex(networkdm,v,&ghostvtex)); 140 if (ghostvtex) continue; 141 PetscCall(DMNetworkGetComponent(networkdm,v,ALL_COMPONENTS,NULL,NULL,&nvar)); 142 if (!nvar) continue; 143 144 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&offset)); 145 PetscCall(DMNetworkGetGlobalVertexIndex(networkdm,v,&id)); 146 147 if (rank == 0) { 148 PetscCall(PetscViewerASCIIPrintf(viewer," Vertex %" PetscInt_FMT ":\n",id)); 149 PetscCall(VecArrayPrint_private(viewer,nvar,xv+offset)); 150 } else { 151 values[1] += 1; /* number of vertices */ 152 values[k++] = id; 153 values[k++] = nvar; 154 for (i=offset; i< offset+nvar; i++) values[k++] = xv[i]; 155 } 156 } 157 158 if (rank == 0) { 159 /* proc[0] receives and prints messages */ 160 for (j=1; j<size; j++) { 161 PetscCall(PetscViewerASCIIPrintf(viewer,"Process [%" PetscInt_FMT "]\n",j)); 162 163 PetscCallMPI(MPI_Recv(values,(PetscMPIInt)len,MPIU_SCALAR,j,tag,comm,&status)); 164 165 ne = (PetscInt)PetscAbsScalar(values[0]); 166 nv = (PetscInt)PetscAbsScalar(values[1]); 167 168 /* print received edges */ 169 k = 2; 170 for (i=0; i<ne; i++) { 171 id = (PetscInt)PetscAbsScalar(values[k++]); 172 nvar = (PetscInt)PetscAbsScalar(values[k++]); 173 PetscCall(PetscViewerASCIIPrintf(viewer," Edge %" PetscInt_FMT ":\n",id)); 174 PetscCall(VecArrayPrint_private(viewer,nvar,values+k)); 175 k += nvar; 176 } 177 178 /* print received vertices */ 179 for (i=0; i<nv; i++) { 180 id = (PetscInt)PetscAbsScalar(values[k++]); 181 nvar = (PetscInt)PetscAbsScalar(values[k++]); 182 PetscCall(PetscViewerASCIIPrintf(viewer," Vertex %" PetscInt_FMT ":\n",id)); 183 PetscCall(VecArrayPrint_private(viewer,nvar,values+k)); 184 k += nvar; 185 } 186 } 187 } else { 188 /* sends values to proc[0] */ 189 PetscCallMPI(MPI_Send((void*)values,k,MPIU_SCALAR,0,tag,comm)); 190 } 191 192 PetscCall(PetscFree(values)); 193 PetscCall(VecRestoreArrayRead(localX,&xv)); 194 PetscCall(DMRestoreLocalVector(networkdm,&localX)); 195 PetscFunctionReturn(0); 196 } 197 198 PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer); 199 200 PetscErrorCode VecView_Network(Vec v,PetscViewer viewer) 201 { 202 DM dm; 203 PetscBool isseq; 204 PetscBool iascii; 205 206 PetscFunctionBegin; 207 PetscCall(VecGetDM(v,&dm)); 208 PetscCheck(dm,PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_WRONG,"Vector not generated from a DM"); 209 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 210 PetscCall(PetscObjectTypeCompare((PetscObject)v,VECSEQ,&isseq)); 211 212 /* Use VecView_Network if the viewer is ASCII; use VecView_Seq/MPI for other viewer formats */ 213 if (iascii) { 214 if (isseq) PetscCall(VecView_Network_Seq(dm,v,viewer)); 215 else PetscCall(VecView_Network_MPI(dm,v,viewer)); 216 } else { 217 if (isseq) PetscCall(VecView_Seq(v,viewer)); 218 else PetscCall(VecView_MPI(v,viewer)); 219 } 220 PetscFunctionReturn(0); 221 } 222 223 static PetscErrorCode DMCreateGlobalVector_Network(DM dm,Vec *vec) 224 { 225 DM_Network *network = (DM_Network*) dm->data; 226 227 PetscFunctionBegin; 228 PetscCall(DMCreateGlobalVector(network->plex,vec)); 229 PetscCall(VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Network)); 230 PetscCall(VecSetDM(*vec,dm)); 231 PetscFunctionReturn(0); 232 } 233 234 static PetscErrorCode DMCreateLocalVector_Network(DM dm,Vec *vec) 235 { 236 DM_Network *network = (DM_Network*) dm->data; 237 238 PetscFunctionBegin; 239 PetscCall(DMCreateLocalVector(network->plex,vec)); 240 PetscCall(VecSetDM(*vec,dm)); 241 PetscFunctionReturn(0); 242 } 243 244 PetscErrorCode DMNetworkInitializeToDefault_NonShared(DM dm) 245 { 246 DM_Network *network = (DM_Network*) dm->data; 247 248 PetscFunctionBegin; 249 network->Je = NULL; 250 network->Jv = NULL; 251 network->Jvptr = NULL; 252 network->userEdgeJacobian = PETSC_FALSE; 253 network->userVertexJacobian = PETSC_FALSE; 254 255 network->vertex.DofSection = NULL; 256 network->vertex.GlobalDofSection = NULL; 257 network->vertex.mapping = NULL; 258 network->vertex.sf = NULL; 259 260 network->edge.DofSection = NULL; 261 network->edge.GlobalDofSection = NULL; 262 network->edge.mapping = NULL; 263 network->edge.sf = NULL; 264 265 network->DataSection = NULL; 266 network->DofSection = NULL; 267 network->GlobalDofSection = NULL; 268 network->componentsetup = PETSC_FALSE; 269 270 network->plex = NULL; 271 272 network->component = NULL; 273 network->ncomponent = 0; 274 275 network->header = NULL; 276 network->cvalue = NULL; 277 network->componentdataarray = NULL; 278 279 network->max_comps_registered = DMNETWORK_MAX_COMP_REGISTERED_DEFAULT; /* return to default */ 280 PetscFunctionReturn(0); 281 } 282 /* Default values for the parameters in DMNetwork */ 283 PetscErrorCode DMNetworkInitializeToDefault(DM dm) 284 { 285 DM_Network *network = (DM_Network*) dm->data; 286 DMNetworkCloneShared cloneshared = network->cloneshared; 287 288 PetscFunctionBegin; 289 PetscCall(DMNetworkInitializeToDefault_NonShared(dm)); 290 /* Default values for shared data */ 291 cloneshared->refct = 1; 292 cloneshared->NVertices = 0; 293 cloneshared->NEdges = 0; 294 cloneshared->nVertices = 0; 295 cloneshared->nEdges = 0; 296 cloneshared->nsubnet = 0; 297 cloneshared->pStart = -1; 298 cloneshared->pEnd = -1; 299 cloneshared->vStart = -1; 300 cloneshared->vEnd = -1; 301 cloneshared->eStart = -1; 302 cloneshared->eEnd = -1; 303 cloneshared->vltog = NULL; 304 cloneshared->distributecalled = PETSC_FALSE; 305 306 cloneshared->subnet = NULL; 307 cloneshared->subnetvtx = NULL; 308 cloneshared->subnetedge = NULL; 309 cloneshared->svtx = NULL; 310 cloneshared->nsvtx = 0; 311 cloneshared->Nsvtx = 0; 312 cloneshared->svertices = NULL; 313 cloneshared->sedgelist = NULL; 314 cloneshared->svtable = NULL; 315 PetscFunctionReturn(0); 316 } 317 318 PetscErrorCode DMInitialize_Network(DM dm) 319 { 320 PetscFunctionBegin; 321 PetscCall(DMSetDimension(dm,1)); 322 dm->ops->view = DMView_Network; 323 dm->ops->setfromoptions = DMSetFromOptions_Network; 324 dm->ops->clone = DMClone_Network; 325 dm->ops->setup = DMSetUp_Network; 326 dm->ops->createglobalvector = DMCreateGlobalVector_Network; 327 dm->ops->createlocalvector = DMCreateLocalVector_Network; 328 dm->ops->getlocaltoglobalmapping = NULL; 329 dm->ops->createfieldis = NULL; 330 dm->ops->createcoordinatedm = NULL; 331 dm->ops->getcoloring = NULL; 332 dm->ops->creatematrix = DMCreateMatrix_Network; 333 dm->ops->createinterpolation = NULL; 334 dm->ops->createinjection = NULL; 335 dm->ops->refine = NULL; 336 dm->ops->coarsen = NULL; 337 dm->ops->refinehierarchy = NULL; 338 dm->ops->coarsenhierarchy = NULL; 339 dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network; 340 dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network; 341 dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network; 342 dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network; 343 dm->ops->destroy = DMDestroy_Network; 344 dm->ops->createsubdm = NULL; 345 dm->ops->locatepoints = NULL; 346 PetscFunctionReturn(0); 347 } 348 /* 349 copies over the subnetid and index portions of the DMNetworkComponentHeader from orignal dm to the newdm 350 */ 351 static PetscErrorCode DMNetworkCopyHeaderTopological(DM dm, DM newdm) 352 { 353 DM_Network *network = (DM_Network *) dm->data, *newnetwork = (DM_Network *) newdm->data; 354 PetscInt p,i,np,index,subnetid; 355 356 PetscFunctionBegin; 357 np = network->cloneshared->pEnd - network->cloneshared->pStart; 358 PetscCall(PetscCalloc2(np,&newnetwork->header,np,&newnetwork->cvalue)); 359 for (i = 0; i < np; i++ ) { 360 p = i + network->cloneshared->pStart; 361 PetscCall(DMNetworkGetSubnetID(dm,p,&subnetid)); 362 PetscCall(DMNetworkGetIndex(dm,p,&index)); 363 newnetwork->header[i].index = index; 364 newnetwork->header[i].subnetid = subnetid; 365 newnetwork->header[i].size = NULL; 366 newnetwork->header[i].key = NULL; 367 newnetwork->header[i].offset = NULL; 368 newnetwork->header[i].nvar = NULL; 369 newnetwork->header[i].offsetvarrel = NULL; 370 newnetwork->header[i].ndata = 0; 371 newnetwork->header[i].maxcomps = DMNETWORK_MAX_COMP_AT_POINT_DEFAULT; 372 newnetwork->header[i].hsize = sizeof(struct _p_DMNetworkComponentHeader)/sizeof(sizeof(DMNetworkComponentGenericDataType)); 373 } 374 PetscFunctionReturn(0); 375 } 376 377 PetscErrorCode DMClone_Network(DM dm, DM *newdm) 378 { 379 DM_Network *network = (DM_Network *) dm->data, *newnetwork = NULL; 380 381 PetscFunctionBegin; 382 network->cloneshared->refct++; 383 PetscCall(PetscNewLog(*newdm,&newnetwork)); 384 (*newdm)->data = newnetwork; 385 PetscCall(DMNetworkInitializeToDefault_NonShared(*newdm)); 386 newnetwork->cloneshared = network->cloneshared; /* Share all data that can be cloneshared */ 387 PetscCall(DMClone(network->plex,&newnetwork->plex)); 388 PetscCall(DMNetworkCopyHeaderTopological(dm, *newdm)); 389 PetscCall(DMNetworkInitializeNonTopological(*newdm)); /* initialize all non-topological data to the state after DMNetworkLayoutSetUp as been called */ 390 PetscCall(PetscObjectChangeTypeName((PetscObject) *newdm, DMNETWORK)); 391 PetscCall(DMInitialize_Network(*newdm)); 392 PetscFunctionReturn(0); 393 } 394 395 /*MC 396 DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object 397 DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of 398 the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary. 399 This is specified by a PetscSection object. Ownership in the global representation is determined by 400 ownership of the underlying DMPlex points. This is specified by another PetscSection object. 401 402 Level: intermediate 403 404 .seealso: `DMType`, `DMNetworkCreate()`, `DMCreate()`, `DMSetType()` 405 M*/ 406 407 PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm) 408 { 409 DM_Network *network; 410 411 PetscFunctionBegin; 412 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 413 PetscCall(PetscNewLog(dm,&network)); 414 PetscCall(PetscNewLog(dm,&network->cloneshared)); 415 dm->data = network; 416 417 PetscCall(DMNetworkInitializeToDefault(dm)); 418 PetscCall(DMInitialize_Network(dm)); 419 PetscFunctionReturn(0); 420 } 421 422 /*@ 423 DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network. 424 425 Collective 426 427 Input Parameter: 428 . comm - The communicator for the DMNetwork object 429 430 Output Parameter: 431 . network - The DMNetwork object 432 433 Level: beginner 434 435 @*/ 436 PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network) 437 { 438 PetscFunctionBegin; 439 PetscValidPointer(network,2); 440 PetscCall(DMCreate(comm, network)); 441 PetscCall(DMSetType(*network, DMNETWORK)); 442 PetscFunctionReturn(0); 443 } 444