1 #include <petsc/private/sfimpl.h> /*I "petscsf.h" I*/ 2 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 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 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 @*/ 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 PetscCall(PetscObjectGetComm((PetscObject)sf, &comm)); 130 PetscCall(GetBoundingBox_Internal(nroots, dim, rootcoords, bbox)); 131 PetscCallMPI(MPI_Comm_size(comm, &size)); 132 PetscCall(PetscMalloc1(size * 2 * dim, &bboxes)); 133 PetscCall(PetscMPIIntCast(2 * dim, &msize)); 134 PetscCallMPI(MPI_Allgather(bbox, msize, MPIU_REAL, bboxes, msize, MPIU_REAL, comm)); 135 PetscCall(GetBoundingBox_Internal(nleaves, dim, leafcoords, bbox)); 136 PetscCall(PetscMalloc1(size, &root_sizes)); 137 PetscCallMPI(MPI_Allgather(&nroots, 1, MPIU_INT, root_sizes, 1, MPIU_INT, comm)); 138 139 PetscCall(PetscMalloc2(size, &ranks_needed, size + 1, &root_starts)); 140 root_starts[0] = 0; 141 num_ranks = 0; 142 for (PetscMPIInt r = 0; r < size; r++) { 143 if (IntersectBoundingBox_Internal(dim, bbox, &bboxes[2 * dim * r], tol)) { 144 ranks_needed[num_ranks++] = r; 145 root_starts[num_ranks] = root_starts[num_ranks - 1] + root_sizes[r]; 146 } 147 } 148 PetscCall(PetscFree(root_sizes)); 149 PetscCall(PetscMalloc1(root_starts[num_ranks], &premote)); 150 for (PetscInt i = 0; i < num_ranks; i++) { 151 for (PetscInt j = root_starts[i]; j < root_starts[i + 1]; j++) { 152 premote[j].rank = ranks_needed[i]; 153 premote[j].index = j - root_starts[i]; 154 } 155 } 156 PetscCall(PetscSFCreate(comm, &psf)); 157 PetscCall(PetscSFSetGraph(psf, nroots, root_starts[num_ranks], NULL, PETSC_USE_POINTER, premote, PETSC_USE_POINTER)); 158 PetscCall(PetscMalloc1(root_starts[num_ranks] * dim, &target_coords)); 159 PetscCall(PetscMPIIntCast(dim, &msize)); 160 PetscCallMPI(MPI_Type_contiguous(msize, MPIU_REAL, &unit)); 161 PetscCallMPI(MPI_Type_commit(&unit)); 162 PetscCall(PetscSFBcastBegin(psf, unit, rootcoords, target_coords, MPI_REPLACE)); 163 PetscCall(PetscSFBcastEnd(psf, unit, rootcoords, target_coords, MPI_REPLACE)); 164 PetscCallMPI(MPI_Type_free(&unit)); 165 PetscCall(PetscSFDestroy(&psf)); 166 167 // Condense targets to only those that lie within our bounding box 168 PetscInt num_targets = 0; 169 for (PetscInt i = 0; i < root_starts[num_ranks]; i++) { 170 if (InBoundingBox_Internal(dim, &target_coords[i * dim], bbox, tol)) { 171 premote[num_targets] = premote[i]; 172 for (PetscInt d = 0; d < dim; d++) target_coords[num_targets * dim + d] = target_coords[i * dim + d]; 173 num_targets++; 174 } 175 } 176 PetscCall(PetscFree(bboxes)); 177 PetscCall(PetscFree2(ranks_needed, root_starts)); 178 179 PetscCall(PetscMalloc1(nleaves, &lremote)); 180 PetscKDTree tree; 181 PetscCount *indices; 182 PetscReal *distances; 183 184 PetscCall(PetscKDTreeCreate(num_targets, dim, target_coords, PETSC_USE_POINTER, PETSC_DETERMINE, &tree)); 185 PetscCall(PetscMalloc2(nleaves, &indices, nleaves, &distances)); 186 PetscCall(PetscKDTreeQueryPointsNearestNeighbor(tree, nleaves, leafcoords, tol, indices, distances)); 187 for (PetscInt i = 0; i < nleaves; i++) { 188 if (distances[i] < tol) { 189 lremote[i] = premote[indices[i]]; 190 } else { 191 switch (dim) { 192 case 1: 193 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "No target found for leaf coordinate %g", (double)leafcoords[i * dim + 0]); 194 case 2: 195 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]); 196 case 3: 197 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]); 198 } 199 } 200 } 201 PetscCall(PetscFree2(indices, distances)); 202 PetscCall(PetscKDTreeDestroy(&tree)); 203 PetscCall(PetscFree(premote)); 204 PetscCall(PetscFree(target_coords)); 205 PetscCall(PetscSFSetGraph(sf, nroots, nleaves, NULL, PETSC_USE_POINTER, lremote, PETSC_OWN_POINTER)); 206 PetscFunctionReturn(PETSC_SUCCESS); 207 } 208