1 #include <petsc/private/dmnetworkimpl.h> /*I "petscdmnetwork.h" I*/
2 #include <petsc/private/vecimpl.h>
3
DMSetFromOptions_Network(DM dm,PetscOptionItems PetscOptionsObject)4 static PetscErrorCode DMSetFromOptions_Network(DM dm, PetscOptionItems PetscOptionsObject)
5 {
6 PetscFunctionBegin;
7 PetscOptionsHeadBegin(PetscOptionsObject, "DMNetwork Options");
8 PetscOptionsHeadEnd();
9 PetscFunctionReturn(PETSC_SUCCESS);
10 }
11
12 /* External function declarations here */
13 extern PetscErrorCode DMCreateMatrix_Network(DM, Mat *);
14 extern PetscErrorCode DMDestroy_Network(DM);
15 extern PetscErrorCode DMView_Network(DM, PetscViewer);
16 extern PetscErrorCode DMGlobalToLocalBegin_Network(DM, Vec, InsertMode, Vec);
17 extern PetscErrorCode DMGlobalToLocalEnd_Network(DM, Vec, InsertMode, Vec);
18 extern PetscErrorCode DMLocalToGlobalBegin_Network(DM, Vec, InsertMode, Vec);
19 extern PetscErrorCode DMLocalToGlobalEnd_Network(DM, Vec, InsertMode, Vec);
20 extern PetscErrorCode DMSetUp_Network(DM);
21 extern PetscErrorCode DMClone_Network(DM, DM *);
22
VecArrayPrint_private(PetscViewer viewer,PetscInt n,const PetscScalar * xv)23 static PetscErrorCode VecArrayPrint_private(PetscViewer viewer, PetscInt n, const PetscScalar *xv)
24 {
25 PetscInt i;
26
27 PetscFunctionBegin;
28 for (i = 0; i < n; i++) {
29 #if defined(PETSC_USE_COMPLEX)
30 if (PetscImaginaryPart(xv[i]) > 0.0) {
31 PetscCall(PetscViewerASCIIPrintf(viewer, " %g + %g i\n", (double)PetscRealPart(xv[i]), (double)PetscImaginaryPart(xv[i])));
32 } else if (PetscImaginaryPart(xv[i]) < 0.0) {
33 PetscCall(PetscViewerASCIIPrintf(viewer, " %g - %g i\n", (double)PetscRealPart(xv[i]), -(double)PetscImaginaryPart(xv[i])));
34 } else PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)PetscRealPart(xv[i])));
35 #else
36 PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)xv[i]));
37 #endif
38 }
39 PetscFunctionReturn(PETSC_SUCCESS);
40 }
41
VecView_Network_Seq(DM networkdm,Vec X,PetscViewer viewer)42 static PetscErrorCode VecView_Network_Seq(DM networkdm, Vec X, PetscViewer viewer)
43 {
44 PetscInt e, v, Start, End, offset, nvar, id;
45 const PetscScalar *xv;
46
47 PetscFunctionBegin;
48 PetscCall(VecGetArrayRead(X, &xv));
49
50 /* iterate over edges */
51 PetscCall(DMNetworkGetEdgeRange(networkdm, &Start, &End));
52 for (e = Start; e < End; e++) {
53 PetscCall(DMNetworkGetComponent(networkdm, e, ALL_COMPONENTS, NULL, NULL, &nvar));
54 if (!nvar) continue;
55
56 PetscCall(DMNetworkGetLocalVecOffset(networkdm, e, ALL_COMPONENTS, &offset));
57 PetscCall(DMNetworkGetGlobalEdgeIndex(networkdm, e, &id));
58
59 PetscCall(PetscViewerASCIIPrintf(viewer, " Edge %" PetscInt_FMT ":\n", id));
60 PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset));
61 }
62
63 /* iterate over vertices */
64 PetscCall(DMNetworkGetVertexRange(networkdm, &Start, &End));
65 for (v = Start; v < End; v++) {
66 PetscCall(DMNetworkGetComponent(networkdm, v, ALL_COMPONENTS, NULL, NULL, &nvar));
67 if (!nvar) continue;
68
69 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, ALL_COMPONENTS, &offset));
70 PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, v, &id));
71
72 PetscCall(PetscViewerASCIIPrintf(viewer, " Vertex %" PetscInt_FMT ":\n", id));
73 PetscCall(VecArrayPrint_private(viewer, nvar, xv + offset));
74 }
75 PetscCall(PetscViewerFlush(viewer));
76 PetscCall(VecRestoreArrayRead(X, &xv));
77 PetscFunctionReturn(PETSC_SUCCESS);
78 }
79
VecView_Network_MPI(DM networkdm,Vec X,PetscViewer viewer)80 static PetscErrorCode VecView_Network_MPI(DM networkdm, Vec X, PetscViewer viewer)
81 {
82 PetscInt i, e, v, eStart, eEnd, vStart, vEnd, offset, nvar, len_loc, k;
83 const PetscScalar *xv;
84 MPI_Comm comm;
85 PetscMPIInt size, rank, tag = ((PetscObject)viewer)->tag, len;
86 Vec localX;
87 PetscBool ghostvtex;
88 PetscScalar *values;
89 PetscInt ne, nv, id;
90 MPI_Status status;
91
92 PetscFunctionBegin;
93 PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
94 PetscCallMPI(MPI_Comm_size(comm, &size));
95 PetscCallMPI(MPI_Comm_rank(comm, &rank));
96
97 PetscCall(DMGetLocalVector(networkdm, &localX));
98 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
99 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
100 PetscCall(VecGetArrayRead(localX, &xv));
101
102 PetscCall(VecGetLocalSize(localX, &len_loc));
103
104 PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
105 PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
106 len_loc += 2 * (1 + eEnd - eStart + vEnd - vStart);
107
108 /* values = [nedges, nvertices; id, nvar, xedge; ...; id, nvars, xvertex;...], to be sent to proc[0] */
109 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &len_loc, 1, MPIU_INT, MPI_MAX, comm));
110 PetscCall(PetscMPIIntCast(len_loc, &len));
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 (PetscMPIInt j = 1; j < size; j++) {
159 PetscCall(PetscViewerASCIIPrintf(viewer, "Process [%d]\n", j));
160
161 PetscCallMPI(MPI_Recv(values, 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(MPIU_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(PETSC_SUCCESS);
194 }
195
196 PETSC_SINGLE_LIBRARY_INTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
197
VecView_Network(Vec v,PetscViewer viewer)198 static PetscErrorCode VecView_Network(Vec v, PetscViewer viewer)
199 {
200 DM dm;
201 PetscBool isseq;
202 PetscBool isascii;
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, &isascii));
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 (isascii) {
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(PETSC_SUCCESS);
219 }
220
DMCreateGlobalVector_Network(DM dm,Vec * vec)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, (PetscErrorCodeFn *)VecView_Network));
228 PetscCall(VecSetDM(*vec, dm));
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231
DMCreateLocalVector_Network(DM dm,Vec * vec)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(PETSC_SUCCESS);
240 }
241
DMNetworkInitializeToDefault_NonShared(DM dm)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(PETSC_SUCCESS);
279 }
280 /* Default values for the parameters in DMNetwork */
DMNetworkInitializeToDefault(DM dm)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(PETSC_SUCCESS);
314 }
315
DMInitialize_Network(DM dm)316 static 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 = DMCreateCoordinateDM_Network;
329 dm->ops->createcellcoordinatedm = NULL;
330 dm->ops->getcoloring = NULL;
331 dm->ops->creatematrix = DMCreateMatrix_Network;
332 dm->ops->createinterpolation = NULL;
333 dm->ops->createinjection = NULL;
334 dm->ops->refine = NULL;
335 dm->ops->coarsen = NULL;
336 dm->ops->refinehierarchy = NULL;
337 dm->ops->coarsenhierarchy = NULL;
338 dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Network;
339 dm->ops->globaltolocalend = DMGlobalToLocalEnd_Network;
340 dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Network;
341 dm->ops->localtoglobalend = DMLocalToGlobalEnd_Network;
342 dm->ops->destroy = DMDestroy_Network;
343 dm->ops->createsubdm = NULL;
344 dm->ops->locatepoints = NULL;
345 PetscFunctionReturn(PETSC_SUCCESS);
346 }
347 /*
348 copies over the subnetid and index portions of the DMNetworkComponentHeader from original dm to the newdm
349 */
DMNetworkCopyHeaderTopological(DM dm,DM newdm)350 static PetscErrorCode DMNetworkCopyHeaderTopological(DM dm, DM newdm)
351 {
352 DM_Network *network = (DM_Network *)dm->data, *newnetwork = (DM_Network *)newdm->data;
353 PetscInt p, i, np, index, subnetid;
354
355 PetscFunctionBegin;
356 np = network->cloneshared->pEnd - network->cloneshared->pStart;
357 PetscCall(PetscCalloc2(np, &newnetwork->header, np, &newnetwork->cvalue));
358 for (i = 0; i < np; i++) {
359 p = i + network->cloneshared->pStart;
360 PetscCall(DMNetworkGetSubnetID(dm, p, &subnetid));
361 PetscCall(DMNetworkGetIndex(dm, p, &index));
362 newnetwork->header[i].index = index;
363 newnetwork->header[i].subnetid = subnetid;
364 newnetwork->header[i].size = NULL;
365 newnetwork->header[i].key = NULL;
366 newnetwork->header[i].offset = NULL;
367 newnetwork->header[i].nvar = NULL;
368 newnetwork->header[i].offsetvarrel = NULL;
369 newnetwork->header[i].ndata = 0;
370 newnetwork->header[i].maxcomps = DMNETWORK_MAX_COMP_AT_POINT_DEFAULT;
371 newnetwork->header[i].hsize = sizeof(struct _p_DMNetworkComponentHeader) / sizeof(sizeof(DMNetworkComponentGenericDataType));
372 }
373 PetscFunctionReturn(PETSC_SUCCESS);
374 }
375
DMClone_Network(DM dm,DM * newdm)376 PetscErrorCode DMClone_Network(DM dm, DM *newdm)
377 {
378 DM_Network *network = (DM_Network *)dm->data, *newnetwork = NULL;
379
380 PetscFunctionBegin;
381 network->cloneshared->refct++;
382 PetscCall(PetscNew(&newnetwork));
383 (*newdm)->data = newnetwork;
384 PetscCall(DMNetworkInitializeToDefault_NonShared(*newdm));
385 newnetwork->cloneshared = network->cloneshared; /* Share all data that can be cloneshared */
386
387 PetscCheck(network->plex, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_NULL, "Must call DMNetworkLayoutSetUp() first");
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(PETSC_SUCCESS);
394 }
395
396 /* Developer Note: Be aware that the plex inside of the network does not have a coordinate plex.
397 */
DMCreateCoordinateDM_Network(DM dm,DM * cdm)398 PetscErrorCode DMCreateCoordinateDM_Network(DM dm, DM *cdm)
399 {
400 DM_Network *newnetwork = NULL;
401 PetscInt Nf;
402 const char *prefix;
403
404 PetscFunctionBegin;
405 PetscCall(DMClone(dm, cdm));
406 newnetwork = (DM_Network *)(*cdm)->data;
407 PetscCall(DMGetNumFields(newnetwork->plex, &Nf));
408 PetscCall(DMSetNumFields(*cdm, Nf)); /* consistency with the coordinate plex */
409 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
410 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*cdm, prefix));
411 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*cdm, "cdm_"));
412 PetscFunctionReturn(PETSC_SUCCESS);
413 }
414
415 /*MC
416 DMNETWORK = "network" - A DM object that encapsulates an unstructured network. The implementation is based on the DM object
417 DMPlex that manages unstructured grids. Distributed networks use a non-overlapping partitioning of
418 the edges. In the local representation, Vecs contain all unknowns in the interior and shared boundary.
419 This is specified by a PetscSection object. Ownership in the global representation is determined by
420 ownership of the underlying DMPlex points. This is specified by another PetscSection object.
421
422 Level: intermediate
423
424 .seealso: `DMType`, `DMNetworkCreate()`, `DMCreate()`, `DMSetType()`
425 M*/
426
DMCreate_Network(DM dm)427 PETSC_EXTERN PetscErrorCode DMCreate_Network(DM dm)
428 {
429 DM_Network *network;
430
431 PetscFunctionBegin;
432 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
433 PetscCall(PetscNew(&network));
434 PetscCall(PetscNew(&network->cloneshared));
435 dm->data = network;
436
437 PetscCall(DMNetworkInitializeToDefault(dm));
438 PetscCall(DMInitialize_Network(dm));
439 PetscFunctionReturn(PETSC_SUCCESS);
440 }
441
442 /*@
443 DMNetworkCreate - Creates a DMNetwork object, which encapsulates an unstructured network.
444
445 Collective
446
447 Input Parameter:
448 . comm - The communicator for the DMNetwork object
449
450 Output Parameter:
451 . network - The DMNetwork object
452
453 Level: beginner
454
455 .seealso: `DMCreate()`
456 @*/
DMNetworkCreate(MPI_Comm comm,DM * network)457 PetscErrorCode DMNetworkCreate(MPI_Comm comm, DM *network)
458 {
459 PetscFunctionBegin;
460 PetscAssertPointer(network, 2);
461 PetscCall(DMCreate(comm, network));
462 PetscCall(DMSetType(*network, DMNETWORK));
463 PetscFunctionReturn(PETSC_SUCCESS);
464 }
465