Chi-Tech
lagC_03_buildsparsity.cc
Go to the documentation of this file.
2
4
5#include "chi_runtime.h"
6#include "chi_log.h"
7#include "chi_mpi.h"
8
9#include <algorithm>
10
11#define sc_int64 static_cast<int64_t>
12
14{
15
16// ###################################################################
17/**Builds the sparsity pattern for a Continuous Finite Element Method.*/
19 std::vector<int64_t>& nodal_nnz_in_diag,
20 std::vector<int64_t>& nodal_nnz_off_diag,
21 const chi_math::UnknownManager& unknown_manager) const
22{
23 //**************************************** DEFINE UTILITIES
24
25 //=================================== Dof-handler
26 /**Utility mappings*/
27 struct DOFHandler
28 {
29 const int64_t local_block_start = 0;
30 const int64_t local_block_end = 0;
31 const std::vector<uint64_t>& locI_block_addr;
32
33 DOFHandler(int64_t block_start,
34 int64_t block_end,
35 const std::vector<uint64_t>& in_locJ_block_address)
36 : local_block_start(block_start),
37 local_block_end(block_end),
38 locI_block_addr(in_locJ_block_address)
39 {
40 }
41
42 bool IsMapLocal(int64_t ir) const
43 {
44 return (ir >= local_block_start and ir < local_block_end);
45 }
46
47 int64_t MapIRLocal(int64_t ir) const { return ir - local_block_start; }
48
49 int64_t GetLocFromIR(int64_t ir)
50 {
51 int64_t locI =
52 std::upper_bound(locI_block_addr.begin(), locI_block_addr.end(), ir) -
53 locI_block_addr.begin() - 1;
54 return locI;
55 }
56
57 } dof_handler(sc_int64(local_block_address_),
60
61 // Writes a message on ir error
62 auto IR_MAP_ERROR = []()
63 {
64 Chi::log.LogAllError() << "PWL-MapCFEMDOF: ir Mapping error node ";
65 Chi::Exit(EXIT_FAILURE);
66 };
67
68 // Writes a message on jr error
69 auto JR_MAP_ERROR = []()
70 {
71 Chi::log.LogAllError() << "PWL-MapCFEMDOF: jr Mapping error node ";
72 Chi::Exit(EXIT_FAILURE);
73 };
74
75 // Checks whether an integer is already in a vector
76 auto IS_VALUE_IN_VECTOR = [](const std::vector<int64_t>& vec, int64_t val)
77 {
78 bool already_there = false;
79 for (auto check_val : vec)
80 if (check_val == val)
81 {
82 already_there = true;
83 break;
84 }
85 return already_there;
86 };
87
88 //**************************************** END OF UTILITIES
89
90 //======================================== Build local sparsity pattern
91 Chi::log.Log0Verbose1() << "Building local sparsity pattern.";
92 std::vector<std::vector<int64_t>> nodal_connections(local_base_block_size_);
93
94 nodal_nnz_in_diag.clear();
95 nodal_nnz_off_diag.clear();
96
97 nodal_nnz_in_diag.resize(local_base_block_size_, 0);
98 nodal_nnz_off_diag.resize(local_base_block_size_, 0);
99
100 for (auto& cell : ref_grid_.local_cells)
101 {
102 const auto& cell_mapping = GetCellMapping(cell);
103 for (unsigned int i = 0; i < cell_mapping.NumNodes(); ++i)
104 {
105 const int64_t ir = MapDOF(cell, i);
106 if (ir < 0) IR_MAP_ERROR();
107
108 if (dof_handler.IsMapLocal(ir))
109 {
110 const int64_t il = dof_handler.MapIRLocal(ir);
111 std::vector<int64_t>& node_links = nodal_connections[il];
112
113 for (unsigned int j = 0; j < cell_mapping.NumNodes(); ++j)
114 {
115 const int64_t jr = MapDOF(cell, j);
116 if (jr < 0) JR_MAP_ERROR();
117
118 if (IS_VALUE_IN_VECTOR(node_links, jr)) continue;
119
120 node_links.push_back(jr);
121 if (dof_handler.IsMapLocal(jr)) nodal_nnz_in_diag[il] += 1;
122 else
123 nodal_nnz_off_diag[il] += 1;
124 } // for j
125 } // if i local
126 } // for i
127 } // for cell
128
129 //======================================== Build non-local sparsity pattern
130 Chi::log.Log0Verbose1() << "Building non-local sparsity pattern.";
131
132 // In this process we build a list
133 // of ir-nodes that are not local. Each ir-node needs to
134 // be furnished with the jr-nodes it links to.
135
136 typedef std::pair<int64_t, std::vector<int64_t>> ROWJLINKS;
137 std::vector<ROWJLINKS> ir_links;
138
139 for (auto& cell : ref_grid_.local_cells)
140 {
141 const auto& cell_mapping = GetCellMapping(cell);
142
143 for (unsigned int i = 0; i < cell_mapping.NumNodes(); ++i)
144 {
145 const int64_t ir = MapDOF(cell, i);
146 if (ir < 0) IR_MAP_ERROR();
147
148 if (not dof_handler.IsMapLocal(ir))
149 {
150 ROWJLINKS new_ir_link;
151 ROWJLINKS* cur_ir_link = &new_ir_link;
152
153 //============================= Check if ir already there
154 bool ir_already_there = false;
155 for (auto& ir_link : ir_links)
156 if (ir == ir_link.first)
157 {
158 ir_already_there = true;
159 cur_ir_link = &ir_link;
160 break;
161 }
162
163 //============================= Now add links
164 auto& node_links = cur_ir_link->second;
165 for (unsigned int j = 0; j < cell_mapping.NumNodes(); ++j)
166 {
167 const int64_t jr = MapDOF(cell, j);
168 if (jr < 0) JR_MAP_ERROR();
169
170 if (IS_VALUE_IN_VECTOR(node_links, jr)) continue;
171 else
172 node_links.push_back(jr);
173 } // for j
174
175 //============================= If its not add it
176 if (not ir_already_there)
177 {
178 cur_ir_link->first = ir;
179 ir_links.push_back(*cur_ir_link);
180 }
181
182 } // if i not local
183 } // for i
184 } // for cell
185
186 //======================================== Build communication structure
187 Chi::log.Log0Verbose1() << "Building communication structure.";
188
189 //=================================== Step 1
190 // We now serialize the non-local data
191 std::vector<std::vector<int64_t>> locI_serialized(Chi::mpi.process_count);
192
193 for (const auto& ir_linkage : ir_links)
194 {
195 const int64_t locI = dof_handler.GetLocFromIR(ir_linkage.first);
196
197 locI_serialized[locI].push_back(
198 sc_int64(ir_linkage.second.size())); // row cols amount
199 locI_serialized[locI].push_back(sc_int64(ir_linkage.first)); // row num
200 for (int64_t jr : ir_linkage.second)
201 locI_serialized[locI].push_back(jr); // col num
202 }
203
204 //=================================== Step 2
205 // Establish the size of the serialized data
206 // to send to each location and communicate
207 // to get receive count.
208 std::vector<int> sendcount(Chi::mpi.process_count, 0);
209 std::vector<int> recvcount(Chi::mpi.process_count, 0);
210 int locI = 0;
211 for (const auto& locI_data : locI_serialized)
212 {
213 sendcount[locI] = static_cast<int>(locI_data.size());
214
215 if (Chi::mpi.location_id == 0)
217 << "To send to " << locI << " = " << sendcount[locI];
218
219 ++locI;
220 }
221
222 MPI_Alltoall(
223 sendcount.data(), 1, MPI_INT, recvcount.data(), 1, MPI_INT, Chi::mpi.comm);
224
225 //=================================== Step 3
226 // We now establish send displacements and
227 // receive displacements.
228 std::vector<int> send_displs(Chi::mpi.process_count, 0);
229 std::vector<int> recv_displs(Chi::mpi.process_count, 0);
230
231 int send_displ_c = 0;
232 int recv_displ_c = 0;
233
234 int c = 0;
235 for (int send_count : sendcount)
236 {
237 send_displs[c++] = send_displ_c;
238 send_displ_c += send_count;
239 }
240 c = 0;
241 for (int recv_count : recvcount)
242 {
243 recv_displs[c++] = recv_displ_c;
244 recv_displ_c += recv_count;
245 }
246
247 //======================================== Communicate data
248 Chi::log.Log0Verbose1() << "Communicating non-local rows.";
249
250 // We now initialize the buffers and
251 // communicate the data
252 std::vector<int64_t> sendbuf;
253 std::vector<int64_t> recvbuf;
254
255 sendbuf.reserve(send_displ_c);
256 recvbuf.resize(recv_displ_c, 0);
257
258 for (const auto& serial_block : locI_serialized)
259 for (int64_t data_val : serial_block)
260 sendbuf.push_back(data_val);
261
262 MPI_Alltoallv(sendbuf.data(),
263 sendcount.data(),
264 send_displs.data(),
265 MPI_INT64_T,
266 recvbuf.data(),
267 recvcount.data(),
268 recv_displs.data(),
269 MPI_INT64_T,
270 Chi::mpi.comm);
271
272 //======================================== Deserialze data
273 Chi::log.Log0Verbose1() << "Deserialize data.";
274
275 std::vector<ROWJLINKS> foreign_ir_links;
276
277 for (size_t k = 0; k < recvbuf.size();)
278 {
279 const int64_t num_values = recvbuf[k++];
280 const int64_t ir = recvbuf[k++];
281
282 ROWJLINKS new_links;
283 new_links.first = ir;
284 new_links.second.reserve(num_values);
285 for (int i = 0; i < num_values; ++i)
286 new_links.second.push_back(recvbuf[k++]);
287
288 foreign_ir_links.push_back(std::move(new_links));
289 }
290
291 //======================================== Adding to sparsity pattern
292 for (const auto& ir_linkage : foreign_ir_links)
293 {
294 const int64_t ir = ir_linkage.first;
295
296 if (not dof_handler.IsMapLocal(ir)) IR_MAP_ERROR();
297
298 int64_t il = dof_handler.MapIRLocal(ir);
299 std::vector<int64_t>& node_links = nodal_connections[il];
300
301 for (int64_t jr : ir_linkage.second)
302 {
303 if (IS_VALUE_IN_VECTOR(node_links, jr)) continue;
304
305 node_links.push_back(jr);
306 if (dof_handler.IsMapLocal(jr)) nodal_nnz_in_diag[il] += 1;
307 else
308 nodal_nnz_off_diag[il] += 1;
309 }
310 }
311
313
314 //======================================== Spacing according to unknown
315 // manager
316 auto backup_nnz_in_diag = nodal_nnz_in_diag;
317 auto backup_nnz_off_diag = nodal_nnz_off_diag;
318
319 unsigned int N = unknown_manager.GetTotalUnknownStructureSize();
320
321 nodal_nnz_in_diag.clear();
322 nodal_nnz_off_diag.clear();
323
324 nodal_nnz_in_diag.resize(local_base_block_size_ * N, 0);
325 nodal_nnz_off_diag.resize(local_base_block_size_ * N, 0);
326
328 {
329 int ir = -1;
330 for (int i = 0; i < local_base_block_size_; ++i)
331 {
332 for (int j = 0; j < N; ++j)
333 {
334 ++ir;
335 nodal_nnz_in_diag[ir] = backup_nnz_in_diag[i];
336 nodal_nnz_off_diag[ir] = backup_nnz_off_diag[i];
337 } // for j
338 } // for i
339 }
340 else if (unknown_manager.dof_storage_type_ ==
342 {
343 int ir = -1;
344 for (int j = 0; j < N; ++j)
345 {
346 for (int i = 0; i < local_base_block_size_; ++i)
347 {
348 ++ir;
349 nodal_nnz_in_diag[ir] = backup_nnz_in_diag[i];
350 nodal_nnz_off_diag[ir] = backup_nnz_off_diag[i];
351 } // for i
352 } // for j
353 }
354}
355
356} // namespace chi_math::spatial_discretization
static void Exit(int error_code)
Definition: chi_runtime.cc:342
static chi::ChiLog & log
Definition: chi_runtime.h:81
static chi::MPI_Info & mpi
Definition: chi_runtime.h:78
LogStream LogAllVerbose1()
Definition: chi_log.h:241
LogStream LogAllError()
Definition: chi_log.h:239
LogStream Log0Verbose1()
Definition: chi_log.h:234
const MPI_Comm & comm
MPI communicator.
Definition: mpi_info.h:28
void Barrier() const
Definition: mpi_info.cc:38
const CellMapping & GetCellMapping(const chi_mesh::Cell &cell) const
std::vector< uint64_t > locJ_block_address_
const chi_mesh::MeshContinuum & ref_grid_
UnknownStorageType dof_storage_type_
unsigned int GetTotalUnknownStructureSize() const
void BuildSparsityPattern(std::vector< int64_t > &nodal_nnz_in_diag, std::vector< int64_t > &nodal_nnz_off_diag, const UnknownManager &unknown_manager) const override
int64_t MapDOF(const chi_mesh::Cell &cell, unsigned int node, const UnknownManager &unknown_manager, unsigned int unknown_id, unsigned int component) const override
LocalCellHandler local_cells
#define sc_int64