1 #include <petsc/private/sfimpl.h> /*I "petscsf.h" I*/
2
GetBoundingBox_Internal(PetscInt npoints,PetscInt dim,const PetscReal * coords,PetscReal * bbox)3 static PetscErrorCode GetBoundingBox_Internal(PetscInt npoints, PetscInt dim, const PetscReal *coords, PetscReal *bbox)
4 {
5 PetscFunctionBegin;
6 for (PetscInt d = 0; d < dim; d++) {
7 bbox[0 * dim + d] = PETSC_MAX_REAL;
8 bbox[1 * dim + d] = PETSC_MIN_REAL;
9 }
10 for (PetscInt i = 0; i < npoints; i++) {
11 for (PetscInt d = 0; d < dim; d++) {
12 bbox[0 * dim + d] = PetscMin(bbox[0 * dim + d], coords[i * dim + d]);
13 bbox[1 * dim + d] = PetscMax(bbox[1 * dim + d], coords[i * dim + d]);
14 }
15 }
16 PetscFunctionReturn(PETSC_SUCCESS);
17 }
18
IntersectBoundingBox_Internal(PetscInt dim,const PetscReal * a,const PetscReal * b,PetscReal tol)19 static PetscBool IntersectBoundingBox_Internal(PetscInt dim, const PetscReal *a, const PetscReal *b, PetscReal tol)
20 {
21 for (PetscInt d = 0; d < dim; d++) {
22 if (a[1 * dim + d] + tol < b[0 * dim + d] || b[1 * dim + d] + tol < a[0 * dim + d]) return PETSC_FALSE;
23 }
24 return PETSC_TRUE;
25 }
26
InBoundingBox_Internal(PetscInt dim,const PetscReal * x,const PetscReal * bbox,PetscReal tol)27 static PetscBool InBoundingBox_Internal(PetscInt dim, const PetscReal *x, const PetscReal *bbox, PetscReal tol)
28 {
29 for (PetscInt d = 0; d < dim; d++) {
30 if (x[d] + tol < bbox[0 * dim + d] || bbox[1 * dim + d] + tol < x[d]) return PETSC_FALSE;
31 }
32 return PETSC_TRUE;
33 }
34
35 /*@
36 PetscSFSetGraphFromCoordinates - Create SF by fuzzy matching leaf coordinates to root coordinates
37
38 Collective
39
40 Input Parameters:
41 + sf - PetscSF to set graph on
42 . nroots - number of root coordinates
43 . nleaves - number of leaf coordinates
44 . dim - spatial dimension of coordinates
45 . tol - positive tolerance for matching
46 . rootcoords - array of root coordinates in which root i component d is [i*dim+d]
47 - leafcoords - array of root coordinates in which leaf i component d is [i*dim+d]
48
49 Notes:
50 The tolerance typically represents the rounding error incurred by numerically computing coordinates via
51 possibly-different procedures. Passing anything from `PETSC_SMALL` to `100 * PETSC_MACHINE_EPSILON` is appropriate for
52 most use cases.
53
54 Example:
55 As a motivating example, consider fluid flow in the x direction with y (distance from a wall). The spanwise direction,
56 z, has periodic boundary conditions and needs some spanwise length to allow turbulent structures to develop. The
57 distribution is stationary with respect to z, so you want to average turbulence variables (like Reynolds stress) over
58 the z direction. It is complicated in a 3D simulation with arbitrary partitioner to uniquely number the nodes or
59 quadrature point coordinates to average these quantities into a 2D plane where they will be visualized, but it's easy
60 to compute the projection of each 3D point into the 2D plane.
61
62 Suppose a 2D target mesh and 3D source mesh (logically an extrusion of the 2D, though perhaps not created in that way)
63 are distributed independently on a communicator. Each rank passes its 2D target points as root coordinates and the 2D
64 projection of its 3D source points as leaf coordinates. Calling `PetscSFReduceBegin()`/`PetscSFReduceEnd()` on the
65 result will sum data from the 3D sources to the 2D targets.
66
67 As a concrete example, consider three MPI ranks with targets (roots)
68 .vb
69 Rank 0: (0, 0), (0, 1)
70 Rank 1: (0.1, 0), (0.1, 1)
71 Rank 2: (0.2, 0), (0.2, 1)
72 .ve
73 Note that targets must be uniquely owned. Suppose also that we identify the following leaf coordinates (perhaps via projection from a 3D space).
74 .vb
75 Rank 0: (0, 0), (0.1, 0), (0, 1), (0.1, 1)
76 Rank 1: (0, 0), (0.1, 0), (0.2, 0), (0, 1), (0.1, 1)
77 Rank 2: (0.1, 0), (0.2, 0), (0.1, 1), (0.2, 1)
78 .ve
79 Leaf coordinates may be repeated, both on a rank and between ranks. This example yields the following `PetscSF` capable of reducing from sources to targets.
80 .vb
81 Roots by rank
82 [0] 0: 0.0000e+00 0.0000e+00 0.0000e+00 1.0000e+00
83 [1] 0: 1.0000e-01 0.0000e+00 1.0000e-01 1.0000e+00
84 [2] 0: 2.0000e-01 0.0000e+00 2.0000e-01 1.0000e+00
85 Leaves by rank
86 [0] 0: 0.0000e+00 0.0000e+00 1.0000e-01 0.0000e+00 0.0000e+00
87 [0] 5: 1.0000e+00 1.0000e-01 1.0000e+00
88 [1] 0: 0.0000e+00 0.0000e+00 1.0000e-01 0.0000e+00 2.0000e-01
89 [1] 5: 0.0000e+00 0.0000e+00 1.0000e+00 1.0000e-01 1.0000e+00
90 [1] 10: 2.0000e-01 1.0000e+00
91 [2] 0: 1.0000e-01 0.0000e+00 2.0000e-01 0.0000e+00 1.0000e-01
92 [2] 5: 1.0000e+00 2.0000e-01 1.0000e+00
93 PetscSF Object: 3 MPI processes
94 type: basic
95 [0] Number of roots=2, leaves=4, remote ranks=2
96 [0] 0 <- (0,0)
97 [0] 1 <- (1,0)
98 [0] 2 <- (0,1)
99 [0] 3 <- (1,1)
100 [1] Number of roots=2, leaves=6, remote ranks=3
101 [1] 0 <- (0,0)
102 [1] 1 <- (1,0)
103 [1] 2 <- (2,0)
104 [1] 3 <- (0,1)
105 [1] 4 <- (1,1)
106 [1] 5 <- (2,1)
107 [2] Number of roots=2, leaves=4, remote ranks=2
108 [2] 0 <- (1,0)
109 [2] 1 <- (2,0)
110 [2] 2 <- (1,1)
111 [2] 3 <- (2,1)
112 .ve
113
114 Level: advanced
115
116 .seealso: `PetscSFCreate()`, `PetscSFSetGraph()`, `PetscSFCreateByMatchingIndices()`
117 @*/
PetscSFSetGraphFromCoordinates(PetscSF sf,PetscInt nroots,PetscInt nleaves,PetscInt dim,PetscReal tol,const PetscReal rootcoords[],const PetscReal leafcoords[])118 PetscErrorCode PetscSFSetGraphFromCoordinates(PetscSF sf, PetscInt nroots, PetscInt nleaves, PetscInt dim, PetscReal tol, const PetscReal rootcoords[], const PetscReal leafcoords[])
119 {
120 PetscReal bbox[6], *bboxes, *target_coords;
121 PetscMPIInt size, *ranks_needed, num_ranks, msize;
122 PetscInt *root_sizes, *root_starts;
123 PetscSFNode *premote, *lremote;
124 PetscSF psf;
125 MPI_Datatype unit;
126 MPI_Comm comm;
127
128 PetscFunctionBegin;
129 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 1);
130 PetscValidLogicalCollectiveInt(sf, dim, 4);
131 PetscValidLogicalCollectiveReal(sf, tol, 5);
132 if (nroots != 0) PetscAssertPointer(rootcoords, 6);
133 if (nleaves != 0) PetscAssertPointer(leafcoords, 7);
134 PetscCall(PetscObjectGetComm((PetscObject)sf, &comm));
135 PetscCall(GetBoundingBox_Internal(nroots, dim, rootcoords, bbox));
136 PetscCallMPI(MPI_Comm_size(comm, &size));
137 PetscCall(PetscMalloc1(size * 2 * dim, &bboxes));
138 PetscCall(PetscMPIIntCast(2 * dim, &msize));
139 PetscCallMPI(MPI_Allgather(bbox, msize, MPIU_REAL, bboxes, msize, MPIU_REAL, comm));
140 PetscCall(GetBoundingBox_Internal(nleaves, dim, leafcoords, bbox));
141 PetscCall(PetscMalloc1(size, &root_sizes));
142 PetscCallMPI(MPI_Allgather(&nroots, 1, MPIU_INT, root_sizes, 1, MPIU_INT, comm));
143
144 PetscCall(PetscMalloc2(size, &ranks_needed, size + 1, &root_starts));
145 root_starts[0] = 0;
146 num_ranks = 0;
147 for (PetscMPIInt r = 0; r < size; r++) {
148 if (IntersectBoundingBox_Internal(dim, bbox, &bboxes[2 * dim * r], tol)) {
149 ranks_needed[num_ranks++] = r;
150 root_starts[num_ranks] = root_starts[num_ranks - 1] + root_sizes[r];
151 }
152 }
153 PetscCall(PetscFree(root_sizes));
154 PetscCall(PetscMalloc1(root_starts[num_ranks], &premote));
155 for (PetscInt i = 0; i < num_ranks; i++) {
156 for (PetscInt j = root_starts[i]; j < root_starts[i + 1]; j++) {
157 premote[j].rank = ranks_needed[i];
158 premote[j].index = j - root_starts[i];
159 }
160 }
161 PetscCall(PetscSFCreate(comm, &psf));
162 PetscCall(PetscSFSetGraph(psf, nroots, root_starts[num_ranks], NULL, PETSC_USE_POINTER, premote, PETSC_USE_POINTER));
163 PetscCall(PetscMalloc1(root_starts[num_ranks] * dim, &target_coords));
164 PetscCall(PetscMPIIntCast(dim, &msize));
165 PetscCallMPI(MPI_Type_contiguous(msize, MPIU_REAL, &unit));
166 PetscCallMPI(MPI_Type_commit(&unit));
167 PetscCall(PetscSFBcastBegin(psf, unit, rootcoords, target_coords, MPI_REPLACE));
168 PetscCall(PetscSFBcastEnd(psf, unit, rootcoords, target_coords, MPI_REPLACE));
169 PetscCallMPI(MPI_Type_free(&unit));
170 PetscCall(PetscSFDestroy(&psf));
171
172 // Condense targets to only those that lie within our bounding box
173 PetscInt num_targets = 0;
174 for (PetscInt i = 0; i < root_starts[num_ranks]; i++) {
175 if (InBoundingBox_Internal(dim, &target_coords[i * dim], bbox, tol)) {
176 premote[num_targets] = premote[i];
177 for (PetscInt d = 0; d < dim; d++) target_coords[num_targets * dim + d] = target_coords[i * dim + d];
178 num_targets++;
179 }
180 }
181 PetscCall(PetscFree(bboxes));
182 PetscCall(PetscFree2(ranks_needed, root_starts));
183
184 PetscCall(PetscMalloc1(nleaves, &lremote));
185 PetscKDTree tree;
186 PetscCount *indices;
187 PetscReal *distances;
188
189 PetscCall(PetscKDTreeCreate(num_targets, dim, target_coords, PETSC_USE_POINTER, PETSC_DETERMINE, &tree));
190 PetscCall(PetscMalloc2(nleaves, &indices, nleaves, &distances));
191 PetscCall(PetscKDTreeQueryPointsNearestNeighbor(tree, nleaves, leafcoords, tol, indices, distances));
192 for (PetscInt i = 0; i < nleaves; i++) {
193 if (distances[i] < tol) {
194 lremote[i] = premote[indices[i]];
195 } else {
196 switch (dim) {
197 case 1:
198 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "No target found for leaf coordinate %g", (double)leafcoords[i * dim + 0]);
199 case 2:
200 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "No target found for leaf coordinate (%g, %g)", (double)leafcoords[i * dim + 0], (double)leafcoords[i * dim + 1]);
201 case 3:
202 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "No target found for leaf coordinate (%g, %g, %g)", (double)leafcoords[i * dim + 0], (double)leafcoords[i * dim + 1], (double)leafcoords[i * dim + 2]);
203 }
204 }
205 }
206 PetscCall(PetscFree2(indices, distances));
207 PetscCall(PetscKDTreeDestroy(&tree));
208 PetscCall(PetscFree(premote));
209 PetscCall(PetscFree(target_coords));
210 PetscCall(PetscSFSetGraph(sf, nroots, nleaves, NULL, PETSC_USE_POINTER, lremote, PETSC_OWN_POINTER));
211 PetscFunctionReturn(PETSC_SUCCESS);
212 }
213