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