Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #include "libmesh/libmesh_common.h"
20 : #include "libmesh/parallel.h"
21 :
22 :
23 : // Local Include
24 : #include "libmesh/libmesh_version.h"
25 : #include "libmesh/system.h"
26 : #include "libmesh/mesh_base.h"
27 : #include "libmesh/elem.h"
28 : #include "libmesh/xdr_cxx.h"
29 : #include "libmesh/numeric_vector.h"
30 : #include "libmesh/dof_map.h"
31 :
32 :
33 : // C++ Includes
34 : #include <memory>
35 : #include <numeric> // for std::partial_sum
36 : #include <set>
37 :
38 :
39 : // Anonymous namespace for implementation details.
40 : namespace {
41 :
42 : using libMesh::DofObject;
43 : using libMesh::Number;
44 : using libMesh::cast_int;
45 :
46 : // Comments:
47 : // ---------
48 : // - The max_io_blksize governs how many nodes or elements will be
49 : // treated as a single block when performing parallel IO on large
50 : // systems.
51 : // - This parameter only loosely affects the size of the actual IO
52 : // buffer as this depends on the number of components a given
53 : // variable has for the nodes/elements in the block.
54 : // - When reading/writing each processor uses an ID map which is
55 : // 3*io_blksize*sizeof(dof_id_type) bytes long, so with unsigned int
56 : // and // io_blksize=256000 we would expect that buffer alone to be
57 : // ~3Mb.
58 : // - In general, an increase in max_io_blksize should increase the
59 : // efficiency of large parallel read/writes by reducing the number
60 : // of MPI messages at the expense of memory.
61 : // - If the library exhausts memory during IO you might reduce this
62 : // parameter.
63 :
64 : const std::size_t max_io_blksize = 256000;
65 :
66 : /**
67 : *
68 : */
69 : template <typename InValType>
70 : class ThreadedIO
71 : {
72 : private:
73 : libMesh::Xdr & _io;
74 : std::vector<InValType> & _data;
75 :
76 : public:
77 45988 : ThreadedIO (libMesh::Xdr & io, std::vector<InValType> & data) :
78 37624 : _io(io),
79 43604 : _data(data)
80 4182 : {}
81 :
82 41806 : void operator()()
83 : {
84 41806 : if (_data.empty()) return;
85 20945 : _io.data_stream (_data.data(), cast_int<unsigned int>(_data.size()));
86 : }
87 : };
88 : }
89 :
90 :
91 : namespace libMesh
92 : {
93 :
94 :
95 : // ------------------------------------------------------------
96 : // System class implementation
97 11283 : void System::read_header (Xdr & io,
98 : std::string_view version,
99 : const bool read_header_in,
100 : const bool read_additional_data,
101 : const bool read_legacy_format)
102 : {
103 : // This method implements the input of a
104 : // System object, embedded in the output of
105 : // an EquationSystems<T_sys>. This warrants some
106 : // documentation. The output file essentially
107 : // consists of 5 sections:
108 : //
109 : // for this system
110 : //
111 : // 5.) The number of variables in the system (unsigned int)
112 : //
113 : // for each variable in the system
114 : //
115 : // 6.) The name of the variable (string)
116 : //
117 : // 6.1.) Variable subdomains
118 : //
119 : // 7.) Combined in an FEType:
120 : // - The approximation order(s) of the variable
121 : // (Order Enum, cast to int/s)
122 : // - The finite element family/ies of the variable
123 : // (FEFamily Enum, cast to int/s)
124 : //
125 : // end variable loop
126 : //
127 : // 8.) The number of additional vectors (unsigned int),
128 : //
129 : // for each additional vector in the system object
130 : //
131 : // 9.) the name of the additional vector (string)
132 : //
133 : // end system
134 324 : libmesh_assert (io.reading());
135 :
136 : // Possibly clear data structures and start from scratch.
137 11283 : if (read_header_in)
138 919 : this->clear ();
139 :
140 : // Figure out if we need to read infinite element information.
141 : // This will be true if the version string contains " with infinite elements"
142 : const bool read_ifem_info =
143 11315 : Utility::contains(version, " with infinite elements") ||
144 11315 : libMesh::on_command_line ("--read-ifem-systems");
145 :
146 :
147 : {
148 : // 5.)
149 : // Read the number of variables in the system
150 11283 : unsigned int nv=0;
151 11607 : if (this->processor_id() == 0)
152 1782 : io.data (nv);
153 11283 : this->comm().broadcast(nv);
154 :
155 11283 : _written_var_indices.clear();
156 11283 : _written_var_indices.resize(nv, 0);
157 :
158 22992 : for (unsigned int var=0; var<nv; var++)
159 : {
160 : // 6.)
161 : // Read the name of the var-th variable
162 672 : std::string var_name;
163 12045 : if (this->processor_id() == 0)
164 1854 : io.data (var_name);
165 11709 : this->comm().broadcast(var_name);
166 :
167 : // 6.1.)
168 672 : std::set<subdomain_id_type> domains;
169 11709 : if (io.version() >= LIBMESH_VERSION_ID(0,7,2))
170 : {
171 672 : std::vector<subdomain_id_type> domain_array;
172 12045 : if (this->processor_id() == 0)
173 1854 : io.data (domain_array);
174 11709 : for (const auto & id : domain_array)
175 0 : domains.insert(id);
176 : }
177 11709 : this->comm().broadcast(domains);
178 :
179 : // 7.)
180 : // Read the approximation order(s) of the var-th variable
181 11709 : int order=0;
182 12045 : if (this->processor_id() == 0)
183 1854 : io.data (order);
184 11709 : this->comm().broadcast(order);
185 :
186 :
187 : // do the same for infinite element radial_order
188 11709 : int rad_order=0;
189 11709 : if (read_ifem_info)
190 : {
191 974 : if (this->processor_id() == 0)
192 334 : io.data(rad_order);
193 654 : this->comm().broadcast(rad_order);
194 : }
195 :
196 : // Read the finite element type of the var-th variable
197 11709 : int fam=0;
198 12045 : if (this->processor_id() == 0)
199 1854 : io.data (fam);
200 11709 : this->comm().broadcast(fam);
201 336 : FEType type;
202 11709 : type.order = static_cast<Order>(order);
203 11709 : type.family = static_cast<FEFamily>(fam);
204 :
205 : // Check for incompatibilities. The shape function indexing was
206 : // changed for the monomial and xyz finite element families to
207 : // simplify extension to arbitrary p. The consequence is that
208 : // old restart files will not be read correctly. This is expected
209 : // to be an unlikely occurrence, but catch it anyway.
210 11709 : if (read_legacy_format)
211 0 : if ((type.family == MONOMIAL || type.family == XYZ) &&
212 0 : ((type.order.get_order() > 2 && this->get_mesh().mesh_dimension() == 2) ||
213 0 : (type.order.get_order() > 1 && this->get_mesh().mesh_dimension() == 3)))
214 : {
215 0 : libmesh_here();
216 0 : libMesh::out << "*****************************************************************\n"
217 0 : << "* WARNING: reading a potentially incompatible restart file!!! *\n"
218 0 : << "* contact libmesh-users@lists.sourceforge.net for more details *\n"
219 0 : << "*****************************************************************"
220 0 : << std::endl;
221 : }
222 :
223 : // Read additional information for infinite elements
224 11709 : int radial_fam=0;
225 11709 : int i_map=0;
226 11709 : if (read_ifem_info)
227 : {
228 974 : if (this->processor_id() == 0)
229 334 : io.data (radial_fam);
230 654 : this->comm().broadcast(radial_fam);
231 974 : if (this->processor_id() == 0)
232 334 : io.data (i_map);
233 654 : this->comm().broadcast(i_map);
234 : }
235 :
236 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
237 :
238 686 : type.radial_order = static_cast<Order>(rad_order);
239 686 : type.radial_family = static_cast<FEFamily>(radial_fam);
240 686 : type.inf_map = static_cast<InfMapType>(i_map);
241 :
242 : #endif
243 :
244 11709 : if (read_header_in)
245 : {
246 1061 : if (domains.empty())
247 1061 : _written_var_indices[var] = this->add_variable (var_name, type);
248 : else
249 0 : _written_var_indices[var] = this->add_variable (var_name, type, &domains);
250 : }
251 : else
252 10648 : _written_var_indices[var] = this->variable_number(var_name);
253 : }
254 : }
255 :
256 : // 8.)
257 : // Read the number of additional vectors.
258 11283 : unsigned int nvecs=0;
259 11607 : if (this->processor_id() == 0)
260 1782 : io.data (nvecs);
261 11283 : this->comm().broadcast(nvecs);
262 :
263 : // If nvecs > 0, this means that write_additional_data
264 : // was true when this file was written. We will need to
265 : // make use of this fact later.
266 11283 : this->_additional_data_written = nvecs;
267 :
268 83383 : for (unsigned int vec=0; vec<nvecs; vec++)
269 : {
270 : // 9.)
271 : // Read the name of the vec-th additional vector
272 4120 : std::string vec_name;
273 74160 : if (this->processor_id() == 0)
274 11330 : io.data (vec_name);
275 72100 : this->comm().broadcast(vec_name);
276 72100 : if (io.version() >= LIBMESH_VERSION_ID(1,7,0))
277 : {
278 72100 : int vec_projection = 0;
279 74160 : if (this->processor_id() == 0)
280 11330 : io.data (vec_projection);
281 72100 : this->comm().broadcast(vec_projection);
282 : int vec_type;
283 74160 : if (this->processor_id() == 0)
284 11330 : io.data (vec_type);
285 72100 : this->comm().broadcast(vec_type);
286 :
287 72100 : if (read_additional_data)
288 74160 : this->add_vector(vec_name, bool(vec_projection), ParallelType(vec_type));
289 : }
290 0 : else if (read_additional_data)
291 : // Systems now can handle adding post-initialization vectors
292 : // libmesh_assert(this->_can_add_vectors);
293 : // Some systems may have added their own vectors already
294 : // libmesh_assert_equal_to (this->_vectors.count(vec_name), 0);
295 0 : this->add_vector(vec_name);
296 : }
297 11283 : }
298 :
299 :
300 :
301 : template <typename InValType>
302 0 : void System::read_parallel_data (Xdr & io,
303 : const bool read_additional_data)
304 : {
305 : /**
306 : * This method implements the output of the vectors
307 : * contained in this System object, embedded in the
308 : * output of an EquationSystems<T_sys>.
309 : *
310 : * 9.) The global solution vector, re-ordered to be node-major
311 : * (More on this later.)
312 : *
313 : * for each additional vector in the object
314 : *
315 : * 10.) The global additional vector, re-ordered to be
316 : * node-major (More on this later.)
317 : *
318 : * Note that the actual IO is handled through the Xdr class
319 : * (to be renamed later?) which provides a uniform interface to
320 : * both the XDR (eXternal Data Representation) interface and standard
321 : * ASCII output. Thus this one section of code will read XDR or ASCII
322 : * files with no changes.
323 : */
324 : // PerfLog pl("IO Performance",false);
325 : // pl.push("read_parallel_data");
326 0 : [[maybe_unused]] dof_id_type total_read_size = 0;
327 :
328 0 : libmesh_assert (io.reading());
329 0 : libmesh_assert (io.is_open());
330 :
331 : // build the ordered nodes and element maps.
332 : // when writing/reading parallel files we need to iterate
333 : // over our nodes/elements in order of increasing global id().
334 : // however, this is not guaranteed to be ordering we obtain
335 : // by using the node_iterators/element_iterators directly.
336 : // so build a set, sorted by id(), that provides the ordering.
337 : // further, for memory economy build the set but then transfer
338 : // its contents to vectors, which will be sorted.
339 0 : std::vector<const DofObject *> ordered_nodes, ordered_elements;
340 : {
341 : std::set<const DofObject *, CompareDofObjectsByID>
342 0 : ordered_nodes_set (this->get_mesh().local_nodes_begin(),
343 0 : this->get_mesh().local_nodes_end());
344 :
345 0 : ordered_nodes.insert(ordered_nodes.end(),
346 : ordered_nodes_set.begin(),
347 : ordered_nodes_set.end());
348 : }
349 : {
350 : std::set<const DofObject *, CompareDofObjectsByID>
351 0 : ordered_elements_set (this->get_mesh().local_elements_begin(),
352 0 : this->get_mesh().local_elements_end());
353 :
354 0 : ordered_elements.insert(ordered_elements.end(),
355 : ordered_elements_set.begin(),
356 : ordered_elements_set.end());
357 : }
358 :
359 : // std::vector<Number> io_buffer;
360 0 : std::vector<InValType> io_buffer;
361 :
362 : // 9.)
363 : //
364 : // Actually read the solution components
365 : // for the ith system to disk
366 0 : io.data(io_buffer);
367 :
368 0 : total_read_size += cast_int<dof_id_type>(io_buffer.size());
369 :
370 0 : const unsigned int sys_num = this->number();
371 0 : const unsigned int nv = cast_int<unsigned int>
372 0 : (this->_written_var_indices.size());
373 0 : libmesh_assert_less_equal (nv, this->n_vars());
374 :
375 0 : dof_id_type cnt=0;
376 :
377 : // Loop over each non-SCALAR variable and each node, and read out the value.
378 0 : for (unsigned int data_var=0; data_var<nv; data_var++)
379 : {
380 0 : const unsigned int var = _written_var_indices[data_var];
381 0 : if (this->variable(var).type().family != SCALAR)
382 : {
383 : // First read the node DOF values
384 0 : for (const auto & node : ordered_nodes)
385 0 : for (auto comp : make_range(node->n_comp(sys_num,var)))
386 : {
387 0 : libmesh_assert_not_equal_to (node->dof_number(sys_num, var, comp),
388 : DofObject::invalid_id);
389 0 : libmesh_assert_less (cnt, io_buffer.size());
390 0 : this->solution->set(node->dof_number(sys_num, var, comp), io_buffer[cnt++]);
391 : }
392 :
393 : // Then read the element DOF values
394 0 : for (const auto & elem : ordered_elements)
395 0 : for (auto comp : make_range(elem->n_comp(sys_num,var)))
396 : {
397 0 : libmesh_assert_not_equal_to (elem->dof_number(sys_num, var, comp),
398 : DofObject::invalid_id);
399 0 : libmesh_assert_less (cnt, io_buffer.size());
400 0 : this->solution->set(elem->dof_number(sys_num, var, comp), io_buffer[cnt++]);
401 : }
402 : }
403 : }
404 :
405 : // Finally, read the SCALAR variables on the last processor
406 0 : for (unsigned int data_var=0; data_var<nv; data_var++)
407 : {
408 0 : const unsigned int var = _written_var_indices[data_var];
409 0 : if (this->variable(var).type().family == SCALAR)
410 : {
411 0 : if (this->processor_id() == (this->n_processors()-1))
412 : {
413 0 : const DofMap & dof_map = this->get_dof_map();
414 0 : std::vector<dof_id_type> SCALAR_dofs;
415 0 : dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
416 :
417 0 : for (auto dof : SCALAR_dofs)
418 0 : this->solution->set(dof, io_buffer[cnt++]);
419 : }
420 : }
421 : }
422 :
423 : // And we're done setting solution entries
424 0 : this->solution->close();
425 :
426 : // For each additional vector, simply go through the list.
427 : // ONLY attempt to do this IF additional data was actually
428 : // written to the file for this system (controlled by the
429 : // _additional_data_written flag).
430 0 : if (this->_additional_data_written)
431 : {
432 0 : const std::size_t nvecs = this->_vectors.size();
433 :
434 : // If the number of additional vectors written is non-zero, and
435 : // the number of additional vectors we have is non-zero, and
436 : // they don't match, then something is wrong and we can't be
437 : // sure we're reading data into the correct places.
438 0 : if (read_additional_data && nvecs &&
439 0 : nvecs != this->_additional_data_written)
440 0 : libmesh_error_msg
441 : ("Additional vectors in file do not match system");
442 :
443 0 : auto pos = _vectors.begin();
444 :
445 0 : for (std::size_t i = 0; i != this->_additional_data_written; ++i)
446 : {
447 0 : cnt=0;
448 0 : io_buffer.clear();
449 :
450 : // 10.)
451 : //
452 : // Actually read the additional vector components
453 : // for the ith system from disk
454 0 : io.data(io_buffer);
455 :
456 0 : total_read_size += cast_int<dof_id_type>(io_buffer.size());
457 :
458 : // If read_additional_data==true and we have additional vectors,
459 : // then we will keep this vector data; otherwise we are going to
460 : // throw it away.
461 0 : if (read_additional_data && nvecs)
462 : {
463 : // Loop over each non-SCALAR variable and each node, and read out the value.
464 0 : for (unsigned int data_var=0; data_var<nv; data_var++)
465 : {
466 0 : const unsigned int var = _written_var_indices[data_var];
467 0 : if (this->variable(var).type().family != SCALAR)
468 : {
469 : // First read the node DOF values
470 0 : for (const auto & node : ordered_nodes)
471 0 : for (auto comp : make_range(node->n_comp(sys_num,var)))
472 : {
473 0 : libmesh_assert_not_equal_to (node->dof_number(sys_num, var, comp),
474 : DofObject::invalid_id);
475 0 : libmesh_assert_less (cnt, io_buffer.size());
476 0 : pos->second->set(node->dof_number(sys_num, var, comp), io_buffer[cnt++]);
477 : }
478 :
479 : // Then read the element DOF values
480 0 : for (const auto & elem : ordered_elements)
481 0 : for (auto comp : make_range(elem->n_comp(sys_num,var)))
482 : {
483 0 : libmesh_assert_not_equal_to (elem->dof_number(sys_num, var, comp),
484 : DofObject::invalid_id);
485 0 : libmesh_assert_less (cnt, io_buffer.size());
486 0 : pos->second->set(elem->dof_number(sys_num, var, comp), io_buffer[cnt++]);
487 : }
488 : }
489 : }
490 :
491 : // Finally, read the SCALAR variables on the last processor
492 0 : for (unsigned int data_var=0; data_var<nv; data_var++)
493 : {
494 0 : const unsigned int var = _written_var_indices[data_var];
495 0 : if (this->variable(var).type().family == SCALAR)
496 : {
497 0 : if (this->processor_id() == (this->n_processors()-1))
498 : {
499 0 : const DofMap & dof_map = this->get_dof_map();
500 0 : std::vector<dof_id_type> SCALAR_dofs;
501 0 : dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
502 :
503 0 : for (auto dof : SCALAR_dofs)
504 0 : pos->second->set(dof, io_buffer[cnt++]);
505 : }
506 : }
507 : }
508 :
509 : // And we're done setting entries for this variable
510 0 : pos->second->close();
511 : }
512 :
513 : // If we've got vectors then we need to be iterating through
514 : // those too
515 0 : if (pos != this->_vectors.end())
516 0 : ++pos;
517 : }
518 : }
519 :
520 : // const Real
521 : // dt = pl.get_elapsed_time(),
522 : // rate = total_read_size*sizeof(Number)/dt;
523 :
524 : // libMesh::err << "Read " << total_read_size << " \"Number\" values\n"
525 : // << " Elapsed time = " << dt << '\n'
526 : // << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
527 :
528 : // pl.pop("read_parallel_data");
529 0 : }
530 :
531 :
532 : template <typename InValType>
533 10999 : void System::read_serialized_data (Xdr & io,
534 : const bool read_additional_data)
535 : {
536 : // This method implements the input of the vectors
537 : // contained in this System object, embedded in the
538 : // output of an EquationSystems<T_sys>.
539 : //
540 : // 10.) The global solution vector, re-ordered to be node-major
541 : // (More on this later.)
542 : //
543 : // for each additional vector in the object
544 : //
545 : // 11.) The global additional vector, re-ordered to be
546 : // node-major (More on this later.)
547 316 : parallel_object_only();
548 632 : std::string comment;
549 :
550 : // PerfLog pl("IO Performance",false);
551 : // pl.push("read_serialized_data");
552 : // std::size_t total_read_size = 0;
553 :
554 : // 10.)
555 : // Read the global solution vector
556 : {
557 : // total_read_size +=
558 10999 : this->read_serialized_vector<InValType>(io, this->solution.get());
559 :
560 : // get the comment
561 11315 : if (this->processor_id() == 0)
562 1734 : io.comment (comment);
563 : }
564 :
565 : // 11.)
566 : // Only read additional vectors if data is available, and only use
567 : // that data to fill our vectors if the user requested it.
568 10999 : if (this->_additional_data_written)
569 : {
570 304 : const std::size_t nvecs = this->_vectors.size();
571 :
572 : // If the number of additional vectors written is non-zero, and
573 : // the number of additional vectors we have is non-zero, and
574 : // they don't match, then we can't read additional vectors
575 : // and be sure we're reading data into the correct places.
576 10640 : if (read_additional_data && nvecs &&
577 10640 : nvecs != this->_additional_data_written)
578 0 : libmesh_error_msg
579 : ("Additional vectors in file do not match system");
580 :
581 304 : auto pos = _vectors.begin();
582 :
583 82740 : for (std::size_t i = 0; i != this->_additional_data_written; ++i)
584 : {
585 : // Read data, but only put it into a vector if we've been
586 : // asked to and if we have a corresponding vector to read.
587 :
588 : // total_read_size +=
589 8240 : this->read_serialized_vector<InValType>
590 138020 : (io, (read_additional_data && nvecs) ? pos->second.get() : nullptr);
591 :
592 : // get the comment
593 74160 : if (this->processor_id() == 0)
594 11330 : io.comment (comment);
595 :
596 :
597 : // If we've got vectors then we need to be iterating through
598 : // those too
599 72100 : if (pos != this->_vectors.end())
600 2060 : ++pos;
601 : }
602 : }
603 :
604 : // const Real
605 : // dt = pl.get_elapsed_time(),
606 : // rate = total_read_size*sizeof(Number)/dt;
607 :
608 : // libMesh::out << "Read " << total_read_size << " \"Number\" values\n"
609 : // << " Elapsed time = " << dt << '\n'
610 : // << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
611 :
612 : // pl.pop("read_serialized_data");
613 10999 : }
614 :
615 :
616 :
617 : template <typename iterator_type, typename InValType>
618 166766 : std::size_t System::read_serialized_blocked_dof_objects (const dof_id_type n_objs,
619 : const iterator_type begin,
620 : const iterator_type end,
621 : const InValType ,
622 : Xdr & io,
623 : const std::vector<NumericVector<Number> *> & vecs,
624 : const unsigned int var_to_read) const
625 : {
626 : //-------------------------------------------------------
627 : // General order: (IO format 0.7.4 & greater)
628 : //
629 : // for (objects ...)
630 : // for (vecs ....)
631 : // for (vars ....)
632 : // for (comps ...)
633 : //
634 : // where objects are nodes or elements, sorted to be
635 : // partition independent,
636 : // vecs are one or more *identically distributed* solution
637 : // coefficient vectors, vars are one or more variables
638 : // to write, and comps are all the components for said
639 : // vars on the object.
640 :
641 : // variables to read. Unless specified otherwise, defaults to _written_var_indices.
642 171534 : std::vector<unsigned int> vars_to_read (_written_var_indices);
643 :
644 166766 : if (var_to_read != libMesh::invalid_uint)
645 0 : vars_to_read.assign({var_to_read});
646 :
647 : const unsigned int
648 9536 : sys_num = this->number(),
649 9536 : num_vecs = cast_int<unsigned int>(vecs.size());
650 : const dof_id_type
651 171550 : io_blksize = cast_int<dof_id_type>(std::min(max_io_blksize, static_cast<std::size_t>(n_objs))),
652 166766 : num_blks = cast_int<unsigned int>(std::ceil(static_cast<double>(n_objs)/
653 157230 : static_cast<double>(io_blksize)));
654 :
655 4768 : libmesh_assert_less_equal (_written_var_indices.size(), this->n_vars());
656 :
657 4768 : std::size_t n_read_values=0;
658 :
659 176302 : std::vector<std::vector<dof_id_type>> xfer_ids(num_blks); // The global IDs and # of components for the local objects in all blocks
660 176302 : std::vector<std::vector<Number>> recv_vals(num_blks); // The raw values for the local objects in all blocks
661 : std::vector<Parallel::Request>
662 181070 : id_requests(num_blks), val_requests(num_blks);
663 : std::vector<Parallel::MessageTag>
664 181070 : id_tags(num_blks), val_tags(num_blks);
665 :
666 : // ------------------------------------------------------
667 : // First pass - count the number of objects in each block
668 : // traverse all the objects and figure out which block they
669 : // will ultimately live in.
670 : std::vector<std::size_t>
671 171534 : xfer_ids_size (num_blks,0),
672 171534 : recv_vals_size (num_blks,0);
673 :
674 :
675 32632620 : for (iterator_type it=begin; it!=end; ++it)
676 : {
677 : const dof_id_type
678 32465854 : id = (*it)->id(),
679 32465854 : block = id/io_blksize;
680 :
681 2951215 : libmesh_assert_less (block, num_blks);
682 :
683 35417069 : xfer_ids_size[block] += 2; // for each object, we send its id, as well as the total number of components for all variables
684 :
685 2951215 : dof_id_type n_comp_tot=0;
686 65180504 : for (const auto & var : vars_to_read)
687 32714650 : n_comp_tot += (*it)->n_comp(sys_num, var); // for each variable, we will receive the nonzero components
688 :
689 35417069 : recv_vals_size[block] += n_comp_tot*num_vecs;
690 : }
691 :
692 : // knowing the recv_vals_size[block] for each processor allows
693 : // us to sum them and find the global size for each block.
694 171534 : std::vector<std::size_t> tot_vals_size(recv_vals_size);
695 166766 : this->comm().sum (tot_vals_size);
696 :
697 :
698 : //------------------------------------------
699 : // Collect the ids & number of values needed
700 : // for all local objects, binning them into
701 : // 'blocks' that will be sent to processor 0
702 333532 : for (dof_id_type blk=0; blk<num_blks; blk++)
703 : {
704 : // Each processor should build up its transfer buffers for its
705 : // local objects in [first_object,last_object).
706 : const dof_id_type
707 166766 : first_object = blk*io_blksize,
708 166766 : last_object = std::min(cast_int<dof_id_type>((blk+1)*io_blksize), n_objs);
709 :
710 : // convenience
711 9552 : std::vector<dof_id_type> & ids (xfer_ids[blk]);
712 9536 : std::vector<Number> & vals (recv_vals[blk]);
713 :
714 : // we now know the number of values we will store for each block,
715 : // so we can do efficient preallocation
716 171534 : ids.clear(); /**/ ids.reserve (xfer_ids_size[blk]);
717 171534 : vals.resize(recv_vals_size[blk]);
718 :
719 : #ifdef DEBUG
720 9536 : std::unordered_set<dof_id_type> seen_ids;
721 : #endif
722 :
723 171534 : if (recv_vals_size[blk] != 0) // only if there are nonzero values to receive
724 26031865 : for (iterator_type it=begin; it!=end; ++it)
725 : {
726 14287886 : dof_id_type id = (*it)->id();
727 : #ifdef DEBUG
728 : // Any renumbering tricks should not have given us any
729 : // duplicate ids.
730 1298479 : libmesh_assert(!seen_ids.count(id));
731 1298479 : seen_ids.insert(id);
732 : #endif
733 :
734 14287886 : if ((id >= first_object) && // object in [first_object,last_object)
735 1298479 : (id < last_object))
736 : {
737 14287886 : ids.push_back(id);
738 :
739 1298479 : unsigned int n_comp_tot=0;
740 :
741 28871020 : for (const auto & var : vars_to_read)
742 14583134 : n_comp_tot += (*it)->n_comp(sys_num, var);
743 :
744 14287886 : ids.push_back (n_comp_tot*num_vecs);
745 : }
746 : }
747 :
748 : #ifdef LIBMESH_HAVE_MPI
749 166766 : id_tags[blk] = this->comm().get_unique_tag(100*num_blks + blk);
750 166766 : val_tags[blk] = this->comm().get_unique_tag(200*num_blks + blk);
751 :
752 : // nonblocking send the data for this block
753 171534 : this->comm().send (0, ids, id_requests[blk], id_tags[blk]);
754 :
755 : // Go ahead and post the receive too
756 171534 : this->comm().receive (0, vals, val_requests[blk], val_tags[blk]);
757 : #endif
758 : }
759 :
760 : //---------------------------------------------------
761 : // Here processor 0 will read and distribute the data.
762 : // We have to do this block-wise to ensure that we
763 : // do not exhaust memory on processor 0.
764 :
765 : // give these variables scope outside the block to avoid reallocation
766 181070 : std::vector<std::vector<dof_id_type>> recv_ids (this->n_processors());
767 181070 : std::vector<std::vector<Number>> send_vals (this->n_processors());
768 181070 : std::vector<Parallel::Request> reply_requests (this->n_processors());
769 9536 : std::vector<unsigned int> obj_val_offsets; // map to traverse entry-wise rather than processor-wise
770 9536 : std::vector<Number> input_vals; // The input buffer for the current block
771 4768 : std::vector<InValType> input_vals_tmp; // The input buffer for the current block
772 :
773 333532 : for (dof_id_type blk=0; blk<num_blks; blk++)
774 : {
775 : // Each processor should build up its transfer buffers for its
776 : // local objects in [first_object,last_object).
777 : const dof_id_type
778 166766 : first_object = blk*io_blksize,
779 166766 : last_object = std::min(cast_int<dof_id_type>((blk+1)*io_blksize), n_objs),
780 166766 : n_objects_blk = last_object - first_object;
781 :
782 : // Processor 0 has a special job. It needs to gather the requested indices
783 : // in [first_object,last_object) from all processors, read the data from
784 : // disk, and reply
785 171534 : if (this->processor_id() == 0)
786 : {
787 : // we know the input buffer size for this block and can begin reading it now
788 28608 : input_vals.resize(tot_vals_size[blk]);
789 28608 : input_vals_tmp.resize(tot_vals_size[blk]);
790 :
791 : // a ThreadedIO object to perform asynchronous file IO
792 2384 : ThreadedIO<InValType> threaded_io(io, input_vals_tmp);
793 28608 : Threads::Thread async_io(threaded_io);
794 :
795 : // offset array. this will define where each object's values
796 : // map into the actual input_vals buffer. this must get
797 : // 0-initialized because 0-component objects are not actually sent
798 26224 : obj_val_offsets.resize (n_objects_blk); /**/ std::fill (obj_val_offsets.begin(), obj_val_offsets.end(), 0);
799 28608 : recv_vals_size.resize(this->n_processors()); // reuse this to count how many values are going to each processor
800 :
801 : #ifndef NDEBUG
802 2384 : std::size_t n_vals_blk = 0;
803 : #endif
804 :
805 : // loop over all processors and process their index request
806 192990 : for (processor_id_type comm_step=0, tnp=this->n_processors(); comm_step != tnp; ++comm_step)
807 : {
808 : #ifdef LIBMESH_HAVE_MPI
809 : // blocking receive indices for this block, imposing no particular order on processor
810 171534 : Parallel::Status id_status (this->comm().probe (Parallel::any_source, id_tags[blk]));
811 166766 : std::vector<dof_id_type> & ids (recv_ids[id_status.source()]);
812 9536 : std::size_t & n_vals_proc (recv_vals_size[id_status.source()]);
813 171534 : this->comm().receive (id_status.source(), ids, id_tags[blk]);
814 : #else
815 : // straight copy without MPI
816 : std::vector<dof_id_type> & ids (recv_ids[0]);
817 : std::size_t & n_vals_proc (recv_vals_size[0]);
818 : ids = xfer_ids[blk];
819 : #endif
820 :
821 166766 : n_vals_proc = 0;
822 :
823 : // note its possible we didn't receive values for objects in
824 : // this block if they have no components allocated.
825 14454652 : for (std::size_t idx=0, sz=ids.size(); idx<sz; idx+=2)
826 : {
827 : const dof_id_type
828 14287886 : local_idx = ids[idx+0]-first_object,
829 14287886 : n_vals_tot_allvecs = ids[idx+1];
830 :
831 1298479 : libmesh_assert_less (local_idx, n_objects_blk);
832 :
833 14287886 : obj_val_offsets[local_idx] = n_vals_tot_allvecs;
834 14287886 : n_vals_proc += n_vals_tot_allvecs;
835 : }
836 :
837 : #ifndef NDEBUG
838 4768 : n_vals_blk += n_vals_proc;
839 : #endif
840 : }
841 :
842 : // We need the offsets into the input_vals vector for each object.
843 : // fortunately, this is simply the partial sum of the total number
844 : // of components for each object
845 23840 : std::partial_sum(obj_val_offsets.begin(), obj_val_offsets.end(),
846 : obj_val_offsets.begin());
847 :
848 2384 : libmesh_assert_equal_to (n_vals_blk, obj_val_offsets.back());
849 2384 : libmesh_assert_equal_to (n_vals_blk, tot_vals_size[blk]);
850 :
851 : // Wait for read completion
852 26224 : async_io.join();
853 : // now copy the values back to the main vector for transfer
854 15661729 : for (auto i_val : index_range(input_vals))
855 18458375 : input_vals[i_val] = input_vals_tmp[i_val];
856 :
857 26224 : n_read_values += input_vals.size();
858 :
859 : // pack data replies for each processor
860 192990 : for (auto proc : make_range(this->n_processors()))
861 : {
862 166766 : const std::vector<dof_id_type> & ids (recv_ids[proc]);
863 9536 : std::vector<Number> & vals (send_vals[proc]);
864 9536 : const std::size_t & n_vals_proc (recv_vals_size[proc]);
865 :
866 166766 : vals.clear(); /**/ vals.reserve(n_vals_proc);
867 :
868 14454652 : for (std::size_t idx=0, sz=ids.size(); idx<sz; idx+=2)
869 : {
870 : const dof_id_type
871 14287886 : local_idx = ids[idx+0]-first_object,
872 15586365 : n_vals_tot_allvecs = ids[idx+1];
873 :
874 1298479 : std::vector<Number>::const_iterator in_vals(input_vals.begin());
875 14287886 : if (local_idx != 0)
876 15572036 : std::advance (in_vals, obj_val_offsets[local_idx-1]);
877 :
878 29923391 : for (unsigned int val=0; val<n_vals_tot_allvecs; val++, ++in_vals)
879 : {
880 1411435 : libmesh_assert (in_vals != input_vals.end());
881 : //libMesh::out << "*in_vals=" << *in_vals << '\n';
882 15635505 : vals.push_back(*in_vals);
883 : }
884 : }
885 :
886 : #ifdef LIBMESH_HAVE_MPI
887 : // send the relevant values to this processor
888 171534 : this->comm().send (proc, vals, reply_requests[proc], val_tags[blk]);
889 : #else
890 : recv_vals[blk] = vals;
891 : #endif
892 : }
893 : } // end processor 0 read/reply
894 :
895 : // all processors complete the (already posted) read for this block
896 : {
897 9552 : Parallel::wait (val_requests[blk]);
898 :
899 9536 : const std::vector<Number> & vals (recv_vals[blk]);
900 9536 : std::vector<Number>::const_iterator val_it(vals.begin());
901 :
902 166766 : if (!recv_vals[blk].empty()) // nonzero values to receive
903 26031865 : for (iterator_type it=begin; it!=end; ++it)
904 27277293 : if (((*it)->id() >= first_object) && // object in [first_object,last_object)
905 14287886 : ((*it)->id() < last_object))
906 : // unpack & set the values
907 30370184 : for (auto & vec : vecs)
908 34573284 : for (const auto & var : vars_to_read)
909 : {
910 18490986 : const unsigned int n_comp = (*it)->n_comp(sys_num, var);
911 :
912 34126491 : for (unsigned int comp=0; comp<n_comp; comp++, ++val_it)
913 : {
914 15635505 : const dof_id_type dof_index = (*it)->dof_number (sys_num, var, comp);
915 1411435 : libmesh_assert (val_it != vals.end());
916 15635505 : if (vec)
917 : {
918 1411435 : libmesh_assert_greater_equal (dof_index, vec->first_local_index());
919 1411435 : libmesh_assert_less (dof_index, vec->last_local_index());
920 : //libMesh::out << "dof_index, *val_it = \t" << dof_index << ", " << *val_it << '\n';
921 15635505 : vec->set (dof_index, *val_it);
922 : }
923 : }
924 : }
925 : }
926 :
927 : // processor 0 needs to make sure all replies have been handed off
928 171534 : if (this->processor_id () == 0)
929 26224 : Parallel::wait(reply_requests);
930 : }
931 :
932 166766 : Parallel::wait(id_requests);
933 :
934 171534 : return n_read_values;
935 314460 : }
936 :
937 :
938 :
939 840 : unsigned int System::read_SCALAR_dofs (const unsigned int var,
940 : Xdr & io,
941 : NumericVector<Number> * vec) const
942 : {
943 24 : unsigned int n_assigned_vals = 0; // the number of values assigned, this will be returned.
944 :
945 : // Processor 0 will read the block from the buffer stream and send it to the last processor
946 840 : const unsigned int n_SCALAR_dofs = this->variable(var).type().order.get_order();
947 840 : std::vector<Number> input_buffer(n_SCALAR_dofs);
948 864 : if (this->processor_id() == 0)
949 132 : io.data_stream(input_buffer.data(), n_SCALAR_dofs);
950 :
951 : #ifdef LIBMESH_HAVE_MPI
952 864 : if (this->n_processors() > 1)
953 : {
954 876 : const Parallel::MessageTag val_tag = this->comm().get_unique_tag();
955 :
956 : // Post the receive on the last processor
957 852 : if (this->processor_id() == (this->n_processors()-1))
958 120 : this->comm().receive(0, input_buffer, val_tag);
959 :
960 : // Send the data to processor 0
961 852 : if (this->processor_id() == 0)
962 120 : this->comm().send(this->n_processors()-1, input_buffer, val_tag);
963 780 : }
964 : #endif
965 :
966 : // Finally, set the SCALAR values
967 864 : if (this->processor_id() == (this->n_processors()-1))
968 : {
969 12 : const DofMap & dof_map = this->get_dof_map();
970 24 : std::vector<dof_id_type> SCALAR_dofs;
971 132 : dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
972 :
973 264 : for (auto i : index_range(SCALAR_dofs))
974 : {
975 132 : if (vec)
976 156 : vec->set (SCALAR_dofs[i], input_buffer[i]);
977 132 : ++n_assigned_vals;
978 : }
979 : }
980 :
981 864 : return n_assigned_vals;
982 : }
983 :
984 :
985 : template <typename InValType>
986 83099 : numeric_index_type System::read_serialized_vector (Xdr & io,
987 : NumericVector<Number> * vec)
988 : {
989 2376 : parallel_object_only();
990 :
991 : #ifndef NDEBUG
992 : // In parallel we better be reading a parallel vector -- if not
993 : // we will not set all of its components below!!
994 2376 : if (this->n_processors() > 1 && vec)
995 : {
996 2376 : libmesh_assert (vec->type() == PARALLEL ||
997 : vec->type() == GHOSTED);
998 : }
999 : #endif
1000 :
1001 2376 : libmesh_assert (io.reading());
1002 :
1003 : // vector length
1004 83099 : unsigned int vector_length=0; // FIXME? size_t would break binary compatibility...
1005 : #ifndef NDEBUG
1006 2376 : std::size_t n_assigned_vals=0;
1007 : #endif
1008 :
1009 : // Get the buffer size
1010 85475 : if (this->processor_id() == 0)
1011 13064 : io.data(vector_length, "# vector length");
1012 83099 : this->comm().broadcast(vector_length);
1013 :
1014 2376 : const unsigned int nv = cast_int<unsigned int>
1015 4752 : (this->_written_var_indices.size());
1016 : const dof_id_type
1017 83099 : n_nodes = this->get_mesh().n_nodes(),
1018 83099 : n_elem = this->get_mesh().n_elem();
1019 :
1020 2376 : libmesh_assert_less_equal (nv, this->n_vars());
1021 :
1022 : // for newer versions, read variables node/elem major
1023 83099 : if (io.version() >= LIBMESH_VERSION_ID(0,7,4))
1024 : {
1025 : //---------------------------------
1026 : // Collect the values for all nodes
1027 : #ifndef NDEBUG
1028 4752 : n_assigned_vals +=
1029 : #endif
1030 163822 : this->read_serialized_blocked_dof_objects (n_nodes,
1031 83099 : this->get_mesh().local_nodes_begin(),
1032 83099 : this->get_mesh().local_nodes_end(),
1033 : InValType(),
1034 : io,
1035 78347 : std::vector<NumericVector<Number> *> (1,vec));
1036 :
1037 :
1038 : //------------------------------------
1039 : // Collect the values for all elements
1040 : #ifndef NDEBUG
1041 4752 : n_assigned_vals +=
1042 : #endif
1043 166198 : this->read_serialized_blocked_dof_objects (n_elem,
1044 83099 : this->get_mesh().local_elements_begin(),
1045 83099 : this->get_mesh().local_elements_end(),
1046 : InValType(),
1047 : io,
1048 156694 : std::vector<NumericVector<Number> *> (1,vec));
1049 : }
1050 :
1051 : // for older versions, read variables var-major
1052 : else
1053 : {
1054 : // Loop over each variable in the system, and then each node/element in the mesh.
1055 0 : for (unsigned int data_var=0; data_var<nv; data_var++)
1056 : {
1057 0 : const unsigned int var = _written_var_indices[data_var];
1058 0 : if (this->variable(var).type().family != SCALAR)
1059 : {
1060 : //---------------------------------
1061 : // Collect the values for all nodes
1062 : #ifndef NDEBUG
1063 0 : n_assigned_vals +=
1064 : #endif
1065 0 : this->read_serialized_blocked_dof_objects (n_nodes,
1066 0 : this->get_mesh().local_nodes_begin(),
1067 0 : this->get_mesh().local_nodes_end(),
1068 : InValType(),
1069 : io,
1070 0 : std::vector<NumericVector<Number> *> (1,vec),
1071 : var);
1072 :
1073 :
1074 : //------------------------------------
1075 : // Collect the values for all elements
1076 : #ifndef NDEBUG
1077 0 : n_assigned_vals +=
1078 : #endif
1079 0 : this->read_serialized_blocked_dof_objects (n_elem,
1080 0 : this->get_mesh().local_elements_begin(),
1081 0 : this->get_mesh().local_elements_end(),
1082 : InValType(),
1083 : io,
1084 0 : std::vector<NumericVector<Number> *> (1,vec),
1085 : var);
1086 : } // end variable loop
1087 : }
1088 : }
1089 :
1090 : //-------------------------------------------
1091 : // Finally loop over all the SCALAR variables
1092 167180 : for (unsigned int data_var=0; data_var<nv; data_var++)
1093 : {
1094 84081 : const unsigned int var = _written_var_indices[data_var];
1095 84081 : if (this->variable(var).type().family == SCALAR)
1096 : {
1097 : #ifndef NDEBUG
1098 24 : n_assigned_vals +=
1099 : #endif
1100 840 : this->read_SCALAR_dofs (var, io, vec);
1101 : }
1102 : }
1103 :
1104 83099 : if (vec)
1105 83099 : vec->close();
1106 :
1107 : #ifndef NDEBUG
1108 2376 : this->comm().sum (n_assigned_vals);
1109 2376 : libmesh_assert_equal_to (n_assigned_vals, vector_length);
1110 : #endif
1111 :
1112 83099 : return vector_length;
1113 : }
1114 :
1115 :
1116 :
1117 1390 : void System::write_header (Xdr & io,
1118 : std::string_view /* version is currently unused */,
1119 : const bool write_additional_data) const
1120 : {
1121 : /**
1122 : * This method implements the output of a
1123 : * System object, embedded in the output of
1124 : * an EquationSystems<T_sys>. This warrants some
1125 : * documentation. The output of this part
1126 : * consists of 5 sections:
1127 : *
1128 : * for this system
1129 : *
1130 : * 5.) The number of variables in the system (unsigned int)
1131 : *
1132 : * for each variable in the system
1133 : *
1134 : * 6.) The name of the variable (string)
1135 : *
1136 : * 6.1.) subdomain where the variable lives
1137 : *
1138 : * 7.) Combined in an FEType:
1139 : * - The approximation order(s) of the variable
1140 : * (Order Enum, cast to int/s)
1141 : * - The finite element family/ies of the variable
1142 : * (FEFamily Enum, cast to int/s)
1143 : *
1144 : * end variable loop
1145 : *
1146 : * 8.) The number of additional vectors (unsigned int),
1147 : *
1148 : * for each additional vector in the system object
1149 : *
1150 : * 9.) the name of the additional vector (string)
1151 : *
1152 : * end system
1153 : */
1154 127 : libmesh_assert (io.writing());
1155 :
1156 :
1157 : // Only write the header information
1158 : // if we are processor 0.
1159 1517 : if (this->get_mesh().processor_id() != 0)
1160 0 : return;
1161 :
1162 254 : std::string comment;
1163 :
1164 : // 5.)
1165 : // Write the number of variables in the system
1166 :
1167 : {
1168 : // set up the comment
1169 127 : comment = "# No. of Variables in System \"";
1170 127 : comment += this->name();
1171 127 : comment += "\"";
1172 :
1173 1390 : unsigned int nv = this->n_vars();
1174 1390 : io.data (nv, comment);
1175 : }
1176 :
1177 :
1178 2852 : for (auto var : make_range(this->n_vars()))
1179 : {
1180 : // 6.)
1181 : // Write the name of the var-th variable
1182 : {
1183 : // set up the comment
1184 133 : comment = "# Name, Variable No. ";
1185 2658 : comment += std::to_string(var);
1186 133 : comment += ", System \"";
1187 133 : comment += this->name();
1188 133 : comment += "\"";
1189 :
1190 1595 : std::string var_name = this->variable_name(var);
1191 1462 : io.data (var_name, comment);
1192 : }
1193 :
1194 : // 6.1.) Variable subdomains
1195 : {
1196 : // set up the comment
1197 133 : comment = "# Subdomains, Variable \"";
1198 1462 : comment += this->variable_name(var);
1199 133 : comment += "\", System \"";
1200 133 : comment += this->name();
1201 133 : comment += "\"";
1202 :
1203 1462 : const std::set<subdomain_id_type> & domains = this->variable(var).active_subdomains();
1204 266 : std::vector<subdomain_id_type> domain_array;
1205 1329 : domain_array.assign(domains.begin(), domains.end());
1206 1462 : io.data (domain_array, comment);
1207 : }
1208 :
1209 : // 7.)
1210 : // Write the approximation order of the var-th variable
1211 : // in this system
1212 : {
1213 : // set up the comment
1214 133 : comment = "# Approximation Order, Variable \"";
1215 1462 : comment += this->variable_name(var);
1216 133 : comment += "\", System \"";
1217 133 : comment += this->name();
1218 133 : comment += "\"";
1219 :
1220 1462 : int order = static_cast<int>(this->variable_type(var).order);
1221 1462 : io.data (order, comment);
1222 : }
1223 :
1224 :
1225 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1226 :
1227 : // do the same for radial_order
1228 : {
1229 133 : comment = "# Radial Approximation Order, Variable \"";
1230 282 : comment += this->variable_name(var);
1231 133 : comment += "\", System \"";
1232 133 : comment += this->name();
1233 133 : comment += "\"";
1234 :
1235 282 : int rad_order = static_cast<int>(this->variable_type(var).radial_order);
1236 282 : io.data (rad_order, comment);
1237 : }
1238 :
1239 : #endif
1240 :
1241 : // Write the Finite Element type of the var-th variable
1242 : // in this System
1243 : {
1244 : // set up the comment
1245 133 : comment = "# FE Family, Variable \"";
1246 1462 : comment += this->variable_name(var);
1247 133 : comment += "\", System \"";
1248 133 : comment += this->name();
1249 133 : comment += "\"";
1250 :
1251 1462 : const FEType & type = this->variable_type(var);
1252 1462 : int fam = static_cast<int>(type.family);
1253 1462 : io.data (fam, comment);
1254 :
1255 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1256 :
1257 133 : comment = "# Radial FE Family, Variable \"";
1258 282 : comment += this->variable_name(var);
1259 133 : comment += "\", System \"";
1260 133 : comment += this->name();
1261 133 : comment += "\"";
1262 :
1263 282 : int radial_fam = static_cast<int>(type.radial_family);
1264 282 : io.data (radial_fam, comment);
1265 :
1266 133 : comment = "# Infinite Mapping Type, Variable \"";
1267 282 : comment += this->variable_name(var);
1268 133 : comment += "\", System \"";
1269 133 : comment += this->name();
1270 133 : comment += "\"";
1271 :
1272 282 : int i_map = static_cast<int>(type.inf_map);
1273 282 : io.data (i_map, comment);
1274 : #endif
1275 : }
1276 : } // end of the variable loop
1277 :
1278 : // 8.)
1279 : // Write the number of additional vectors in the System.
1280 : // If write_additional_data==false, then write zero for
1281 : // the number of additional vectors.
1282 : {
1283 : {
1284 : // set up the comment
1285 127 : comment = "# No. of Additional Vectors, System \"";
1286 127 : comment += this->name();
1287 127 : comment += "\"";
1288 :
1289 2550 : unsigned int nvecs = write_additional_data ? this->n_vectors () : 0;
1290 1390 : io.data (nvecs, comment);
1291 : }
1292 :
1293 1390 : if (write_additional_data)
1294 : {
1295 116 : unsigned int cnt=0;
1296 9768 : for (const auto & [vec_name, vec] : _vectors)
1297 : {
1298 : // 9.)
1299 : // write the name of the cnt-th additional vector
1300 9264 : const std::string dth_vector = std::to_string(cnt++)+"th vector";
1301 9264 : comment = "# Name of " + dth_vector;
1302 8492 : std::string nonconst_vec_name = vec_name; // Stupid XDR API
1303 :
1304 8492 : io.data (nonconst_vec_name, comment);
1305 8492 : int vec_projection = _vector_projections.at(vec_name);
1306 16212 : comment = "# Whether to do projections for " + dth_vector;
1307 8492 : io.data (vec_projection, comment);
1308 8492 : int vec_type = vec->type();
1309 16212 : comment = "# Parallel type of " + dth_vector;
1310 8492 : io.data (vec_type, comment);
1311 : }
1312 : }
1313 : }
1314 : }
1315 :
1316 :
1317 :
1318 0 : void System::write_parallel_data (Xdr & io,
1319 : const bool write_additional_data) const
1320 : {
1321 : /**
1322 : * This method implements the output of the vectors
1323 : * contained in this System object, embedded in the
1324 : * output of an EquationSystems<T_sys>.
1325 : *
1326 : * 9.) The global solution vector, re-ordered to be node-major
1327 : * (More on this later.)
1328 : *
1329 : * for each additional vector in the object
1330 : *
1331 : * 10.) The global additional vector, re-ordered to be
1332 : * node-major (More on this later.)
1333 : *
1334 : * Note that the actual IO is handled through the Xdr class
1335 : * (to be renamed later?) which provides a uniform interface to
1336 : * both the XDR (eXternal Data Representation) interface and standard
1337 : * ASCII output. Thus this one section of code will read XDR or ASCII
1338 : * files with no changes.
1339 : */
1340 : // PerfLog pl("IO Performance",false);
1341 : // pl.push("write_parallel_data");
1342 : // std::size_t total_written_size = 0;
1343 :
1344 0 : std::string comment;
1345 :
1346 0 : libmesh_assert (io.writing());
1347 :
1348 0 : std::vector<Number> io_buffer; io_buffer.reserve(this->solution->local_size());
1349 :
1350 : // build the ordered nodes and element maps.
1351 : // when writing/reading parallel files we need to iterate
1352 : // over our nodes/elements in order of increasing global id().
1353 : // however, this is not guaranteed to be ordering we obtain
1354 : // by using the node_iterators/element_iterators directly.
1355 : // so build a set, sorted by id(), that provides the ordering.
1356 : // further, for memory economy build the set but then transfer
1357 : // its contents to vectors, which will be sorted.
1358 0 : std::vector<const DofObject *> ordered_nodes, ordered_elements;
1359 : {
1360 : std::set<const DofObject *, CompareDofObjectsByID>
1361 0 : ordered_nodes_set (this->get_mesh().local_nodes_begin(),
1362 0 : this->get_mesh().local_nodes_end());
1363 :
1364 0 : ordered_nodes.insert(ordered_nodes.end(),
1365 : ordered_nodes_set.begin(),
1366 0 : ordered_nodes_set.end());
1367 : }
1368 : {
1369 : std::set<const DofObject *, CompareDofObjectsByID>
1370 0 : ordered_elements_set (this->get_mesh().local_elements_begin(),
1371 0 : this->get_mesh().local_elements_end());
1372 :
1373 0 : ordered_elements.insert(ordered_elements.end(),
1374 : ordered_elements_set.begin(),
1375 0 : ordered_elements_set.end());
1376 : }
1377 :
1378 0 : const unsigned int sys_num = this->number();
1379 0 : const unsigned int nv = this->n_vars();
1380 :
1381 : // Loop over each non-SCALAR variable and each node, and write out the value.
1382 0 : for (unsigned int var=0; var<nv; var++)
1383 0 : if (this->variable(var).type().family != SCALAR)
1384 : {
1385 : // First write the node DOF values
1386 0 : for (const auto & node : ordered_nodes)
1387 0 : for (auto comp : make_range(node->n_comp(sys_num,var)))
1388 : {
1389 0 : libmesh_assert_not_equal_to (node->dof_number(sys_num, var, comp),
1390 : DofObject::invalid_id);
1391 :
1392 0 : io_buffer.push_back((*this->solution)(node->dof_number(sys_num, var, comp)));
1393 : }
1394 :
1395 : // Then write the element DOF values
1396 0 : for (const auto & elem : ordered_elements)
1397 0 : for (auto comp : make_range(elem->n_comp(sys_num,var)))
1398 : {
1399 0 : libmesh_assert_not_equal_to (elem->dof_number(sys_num, var, comp),
1400 : DofObject::invalid_id);
1401 :
1402 0 : io_buffer.push_back((*this->solution)(elem->dof_number(sys_num, var, comp)));
1403 : }
1404 : }
1405 :
1406 : // Finally, write the SCALAR data on the last processor
1407 0 : for (auto var : make_range(this->n_vars()))
1408 0 : if (this->variable(var).type().family == SCALAR)
1409 : {
1410 0 : if (this->processor_id() == (this->n_processors()-1))
1411 : {
1412 0 : const DofMap & dof_map = this->get_dof_map();
1413 0 : std::vector<dof_id_type> SCALAR_dofs;
1414 0 : dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
1415 :
1416 0 : for (auto dof : SCALAR_dofs)
1417 0 : io_buffer.push_back((*this->solution)(dof));
1418 : }
1419 : }
1420 :
1421 : // 9.)
1422 : //
1423 : // Actually write the reordered solution vector
1424 : // for the ith system to disk
1425 :
1426 : // set up the comment
1427 : {
1428 0 : comment = "# System \"";
1429 0 : comment += this->name();
1430 0 : comment += "\" Solution Vector";
1431 : }
1432 :
1433 0 : io.data (io_buffer, comment);
1434 :
1435 : // total_written_size += io_buffer.size();
1436 :
1437 : // Only write additional vectors if wanted
1438 0 : if (write_additional_data)
1439 : {
1440 0 : for (auto & [vec_name, vec] : _vectors)
1441 : {
1442 0 : io_buffer.clear();
1443 0 : io_buffer.reserve(vec->local_size());
1444 :
1445 : // Loop over each non-SCALAR variable and each node, and write out the value.
1446 0 : for (unsigned int var=0; var<nv; var++)
1447 0 : if (this->variable(var).type().family != SCALAR)
1448 : {
1449 : // First write the node DOF values
1450 0 : for (const auto & node : ordered_nodes)
1451 0 : for (auto comp : make_range(node->n_comp(sys_num,var)))
1452 : {
1453 0 : libmesh_assert_not_equal_to (node->dof_number(sys_num, var, comp),
1454 : DofObject::invalid_id);
1455 :
1456 0 : io_buffer.push_back((*vec)(node->dof_number(sys_num, var, comp)));
1457 : }
1458 :
1459 : // Then write the element DOF values
1460 0 : for (const auto & elem : ordered_elements)
1461 0 : for (auto comp : make_range(elem->n_comp(sys_num,var)))
1462 : {
1463 0 : libmesh_assert_not_equal_to (elem->dof_number(sys_num, var, comp),
1464 : DofObject::invalid_id);
1465 :
1466 0 : io_buffer.push_back((*vec)(elem->dof_number(sys_num, var, comp)));
1467 : }
1468 : }
1469 :
1470 : // Finally, write the SCALAR data on the last processor
1471 0 : for (auto var : make_range(this->n_vars()))
1472 0 : if (this->variable(var).type().family == SCALAR)
1473 : {
1474 0 : if (this->processor_id() == (this->n_processors()-1))
1475 : {
1476 0 : const DofMap & dof_map = this->get_dof_map();
1477 0 : std::vector<dof_id_type> SCALAR_dofs;
1478 0 : dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
1479 :
1480 0 : for (auto dof : SCALAR_dofs)
1481 0 : io_buffer.push_back((*vec)(dof));
1482 : }
1483 : }
1484 :
1485 : // 10.)
1486 : //
1487 : // Actually write the reordered additional vector
1488 : // for this system to disk
1489 :
1490 : // set up the comment
1491 : {
1492 0 : comment = "# System \"";
1493 0 : comment += this->name();
1494 0 : comment += "\" Additional Vector \"";
1495 0 : comment += vec_name;
1496 0 : comment += "\"";
1497 : }
1498 :
1499 0 : io.data (io_buffer, comment);
1500 :
1501 : // total_written_size += io_buffer.size();
1502 : }
1503 : }
1504 :
1505 : // const Real
1506 : // dt = pl.get_elapsed_time(),
1507 : // rate = total_written_size*sizeof(Number)/dt;
1508 :
1509 : // libMesh::err << "Write " << total_written_size << " \"Number\" values\n"
1510 : // << " Elapsed time = " << dt << '\n'
1511 : // << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
1512 :
1513 : // pl.pop("write_parallel_data");
1514 0 : }
1515 :
1516 :
1517 :
1518 8485 : void System::write_serialized_data (Xdr & io,
1519 : const bool write_additional_data) const
1520 : {
1521 : /**
1522 : * This method implements the output of the vectors
1523 : * contained in this System object, embedded in the
1524 : * output of an EquationSystems<T_sys>.
1525 : *
1526 : * 9.) The global solution vector, re-ordered to be node-major
1527 : * (More on this later.)
1528 : *
1529 : * for each additional vector in the object
1530 : *
1531 : * 10.) The global additional vector, re-ordered to be
1532 : * node-major (More on this later.)
1533 : */
1534 246 : parallel_object_only();
1535 492 : std::string comment;
1536 :
1537 : // PerfLog pl("IO Performance",false);
1538 : // pl.push("write_serialized_data");
1539 : // std::size_t total_written_size = 0;
1540 :
1541 : // total_written_size +=
1542 8485 : this->write_serialized_vector(io, *this->solution);
1543 :
1544 : // set up the comment
1545 8731 : if (this->processor_id() == 0)
1546 : {
1547 123 : comment = "# System \"";
1548 123 : comment += this->name();
1549 123 : comment += "\" Solution Vector";
1550 :
1551 1342 : io.comment (comment);
1552 : }
1553 :
1554 : // Only write additional vectors if wanted
1555 8485 : if (write_additional_data)
1556 : {
1557 62160 : for (auto & pair : this->_vectors)
1558 : {
1559 : // total_written_size +=
1560 54040 : this->write_serialized_vector(io, *pair.second);
1561 :
1562 : // set up the comment
1563 55584 : if (this->processor_id() == 0)
1564 : {
1565 772 : comment = "# System \"";
1566 772 : comment += this->name();
1567 772 : comment += "\" Additional Vector \"";
1568 772 : comment += pair.first;
1569 772 : comment += "\"";
1570 8492 : io.comment (comment);
1571 : }
1572 : }
1573 : }
1574 :
1575 : // const Real
1576 : // dt = pl.get_elapsed_time(),
1577 : // rate = total_written_size*sizeof(Number)/dt;
1578 :
1579 : // libMesh::out << "Write " << total_written_size << " \"Number\" values\n"
1580 : // << " Elapsed time = " << dt << '\n'
1581 : // << " Rate = " << rate/1.e6 << "(MB/sec)\n\n";
1582 :
1583 : // pl.pop("write_serialized_data");
1584 :
1585 :
1586 :
1587 :
1588 : // // test the new method
1589 : // {
1590 : // std::vector<std::string> names;
1591 : // std::vector<NumericVector<Number> *> vectors_to_write;
1592 :
1593 : // names.push_back("Solution Vector");
1594 : // vectors_to_write.push_back(this->solution.get());
1595 :
1596 : // // Only write additional vectors if wanted
1597 : // if (write_additional_data)
1598 : // {
1599 : // std::map<std::string, NumericVector<Number> *>::const_iterator
1600 : // pos = _vectors.begin();
1601 :
1602 : // for (; pos != this->_vectors.end(); ++pos)
1603 : // {
1604 : // names.push_back("Additional Vector " + pos->first);
1605 : // vectors_to_write.push_back(pos->second);
1606 : // }
1607 : // }
1608 :
1609 : // total_written_size =
1610 : // this->write_serialized_vectors (io, names, vectors_to_write);
1611 :
1612 : // const Real
1613 : // dt2 = pl.get_elapsed_time(),
1614 : // rate2 = total_written_size*sizeof(Number)/(dt2-dt);
1615 :
1616 : // libMesh::out << "Write (new) " << total_written_size << " \"Number\" values\n"
1617 : // << " Elapsed time = " << (dt2-dt) << '\n'
1618 : // << " Rate = " << rate2/1.e6 << "(MB/sec)\n\n";
1619 :
1620 : // }
1621 8485 : }
1622 :
1623 :
1624 :
1625 : template <typename iterator_type>
1626 125618 : std::size_t System::write_serialized_blocked_dof_objects (const std::vector<const NumericVector<Number> *> & vecs,
1627 : const dof_id_type n_objs,
1628 : const iterator_type begin,
1629 : const iterator_type end,
1630 : Xdr & io,
1631 : const unsigned int var_to_write) const
1632 : {
1633 3596 : parallel_object_only();
1634 :
1635 : //-------------------------------------------------------
1636 : // General order: (IO format 0.7.4 & greater)
1637 : //
1638 : // for (objects ...)
1639 : // for (vecs ....)
1640 : // for (vars ....)
1641 : // for (comps ...)
1642 : //
1643 : // where objects are nodes or elements, sorted to be
1644 : // partition independent,
1645 : // vecs are one or more *identically distributed* solution
1646 : // coefficient vectors, vars are one or more variables
1647 : // to write, and comps are all the components for said
1648 : // vars on the object.
1649 :
1650 : // We will write all variables unless requested otherwise.
1651 129214 : std::vector<unsigned int> vars_to_write(1, var_to_write);
1652 :
1653 125618 : if (var_to_write == libMesh::invalid_uint)
1654 : {
1655 125618 : vars_to_write.clear(); /**/ vars_to_write.reserve(this->n_vars());
1656 252088 : for (auto var : make_range(this->n_vars()))
1657 126470 : vars_to_write.push_back(var);
1658 : }
1659 :
1660 : const dof_id_type io_blksize = cast_int<dof_id_type>
1661 247640 : (std::min(max_io_blksize, static_cast<std::size_t>(n_objs)));
1662 :
1663 : const unsigned int
1664 7192 : sys_num = this->number(),
1665 7192 : num_vecs = cast_int<unsigned int>(vecs.size()),
1666 125618 : num_blks = cast_int<unsigned int>(std::ceil(static_cast<double>(n_objs)/
1667 118426 : static_cast<double>(io_blksize)));
1668 :
1669 : // libMesh::out << "io_blksize = " << io_blksize
1670 : // << ", num_objects = " << n_objs
1671 : // << ", num_blks = " << num_blks
1672 : // << std::endl;
1673 :
1674 125618 : std::size_t written_length=0; // The numer of values written. This will be returned
1675 132810 : std::vector<std::vector<dof_id_type>> xfer_ids(num_blks); // The global IDs and # of components for the local objects in all blocks
1676 132810 : std::vector<std::vector<Number>> send_vals(num_blks); // The raw values for the local objects in all blocks
1677 : std::vector<Parallel::Request>
1678 136406 : id_requests(num_blks), val_requests(num_blks); // send request handle for each block
1679 : std::vector<Parallel::MessageTag>
1680 136406 : id_tags(num_blks), val_tags(num_blks); // tag number for each block
1681 :
1682 : // ------------------------------------------------------
1683 : // First pass - count the number of objects in each block
1684 : // traverse all the objects and figure out which block they
1685 : // will ultimately live in.
1686 : std::vector<unsigned int>
1687 129214 : xfer_ids_size (num_blks,0),
1688 125618 : send_vals_size (num_blks,0);
1689 :
1690 33656292 : for (iterator_type it=begin; it!=end; ++it)
1691 : {
1692 : const dof_id_type
1693 33530674 : id = (*it)->id(),
1694 33530674 : block = id/io_blksize;
1695 :
1696 3049452 : libmesh_assert_less (block, num_blks);
1697 :
1698 36580126 : xfer_ids_size[block] += 2; // for each object, we store its id, as well as the total number of components for all variables
1699 :
1700 3049452 : unsigned int n_comp_tot=0;
1701 :
1702 67228172 : for (const auto & var : vars_to_write)
1703 33697498 : n_comp_tot += (*it)->n_comp(sys_num, var); // for each variable, we will store the nonzero components
1704 :
1705 36580126 : send_vals_size[block] += n_comp_tot*num_vecs;
1706 : }
1707 :
1708 : //-----------------------------------------
1709 : // Collect the values for all local objects,
1710 : // binning them into 'blocks' that will be
1711 : // sent to processor 0
1712 251236 : for (unsigned int blk=0; blk<num_blks; blk++)
1713 : {
1714 : // libMesh::out << "Writing object block " << blk << std::endl;
1715 :
1716 : // Each processor should build up its transfer buffers for its
1717 : // local objects in [first_object,last_object).
1718 : const dof_id_type
1719 125618 : first_object = blk*io_blksize,
1720 125618 : last_object = std::min(cast_int<dof_id_type>((blk+1)*io_blksize), n_objs);
1721 :
1722 : // convenience
1723 7212 : std::vector<dof_id_type> & ids (xfer_ids[blk]);
1724 7192 : std::vector<Number> & vals (send_vals[blk]);
1725 :
1726 : // we now know the number of values we will store for each block,
1727 : // so we can do efficient preallocation
1728 129214 : ids.clear(); /**/ ids.reserve (xfer_ids_size[blk]);
1729 129214 : vals.clear(); /**/ vals.reserve (send_vals_size[blk]);
1730 :
1731 129214 : if (send_vals_size[blk] != 0) // only send if we have nonzero components to write
1732 26796153 : for (iterator_type it=begin; it!=end; ++it)
1733 28091956 : if (((*it)->id() >= first_object) && // object in [first_object,last_object)
1734 14714978 : ((*it)->id() < last_object))
1735 : {
1736 14714978 : ids.push_back((*it)->id());
1737 :
1738 : // count the total number of nonzeros transferred for this object
1739 : {
1740 1338000 : unsigned int n_comp_tot=0;
1741 :
1742 29582380 : for (const auto & var : vars_to_write)
1743 14867402 : n_comp_tot += (*it)->n_comp(sys_num, var);
1744 :
1745 14714978 : ids.push_back (n_comp_tot*num_vecs); // even if 0 - processor 0 has no way of knowing otherwise...
1746 : }
1747 :
1748 : // pack the values to send
1749 31224368 : for (const auto & vec : vecs)
1750 35284644 : for (const auto & var : vars_to_write)
1751 : {
1752 18775254 : const unsigned int n_comp = (*it)->n_comp(sys_num, var);
1753 :
1754 34875939 : for (unsigned int comp=0; comp<n_comp; comp++)
1755 : {
1756 1454980 : libmesh_assert_greater_equal ((*it)->dof_number(sys_num, var, comp), vec->first_local_index());
1757 1454980 : libmesh_assert_less ((*it)->dof_number(sys_num, var, comp), vec->last_local_index());
1758 16216204 : vals.push_back((*vec)((*it)->dof_number(sys_num, var, comp)));
1759 : }
1760 : }
1761 : }
1762 :
1763 : #ifdef LIBMESH_HAVE_MPI
1764 125618 : id_tags[blk] = this->comm().get_unique_tag(100*num_blks + blk);
1765 125618 : val_tags[blk] = this->comm().get_unique_tag(200*num_blks + blk);
1766 :
1767 : // nonblocking send the data for this block
1768 129214 : this->comm().send (0, ids, id_requests[blk], id_tags[blk]);
1769 129214 : this->comm().send (0, vals, val_requests[blk], val_tags[blk]);
1770 : #endif
1771 : }
1772 :
1773 :
1774 129214 : if (this->processor_id() == 0)
1775 : {
1776 23360 : std::vector<std::vector<dof_id_type>> recv_ids (this->n_processors());
1777 25158 : std::vector<std::vector<Number>> recv_vals (this->n_processors());
1778 3596 : std::vector<unsigned int> obj_val_offsets; // map to traverse entry-wise rather than processor-wise
1779 3596 : std::vector<Number> output_vals; // The output buffer for the current block
1780 :
1781 : // a ThreadedIO object to perform asynchronous file IO
1782 1798 : ThreadedIO<Number> threaded_io(io, output_vals);
1783 17966 : std::unique_ptr<Threads::Thread> async_io;
1784 :
1785 39528 : for (unsigned int blk=0; blk<num_blks; blk++)
1786 : {
1787 : // Each processor should build up its transfer buffers for its
1788 : // local objects in [first_object,last_object).
1789 : const dof_id_type
1790 19764 : first_object = cast_int<dof_id_type>(blk*io_blksize),
1791 19764 : last_object = std::min(cast_int<dof_id_type>((blk+1)*io_blksize), n_objs),
1792 19764 : n_objects_blk = last_object - first_object;
1793 :
1794 : // offset array. this will define where each object's values
1795 : // map into the actual output_vals buffer. this must get
1796 : // 0-initialized because 0-component objects are not actually sent
1797 19764 : obj_val_offsets.resize (n_objects_blk); /**/ std::fill (obj_val_offsets.begin(), obj_val_offsets.end(), 0);
1798 :
1799 1798 : std::size_t n_val_recvd_blk=0;
1800 :
1801 : // receive this block of data from all processors.
1802 145382 : for (processor_id_type comm_step=0, tnp=this->n_processors(); comm_step != tnp; ++comm_step)
1803 : {
1804 : #ifdef LIBMESH_HAVE_MPI
1805 : // blocking receive indices for this block, imposing no particular order on processor
1806 129214 : Parallel::Status id_status (this->comm().probe (Parallel::any_source, id_tags[blk]));
1807 125618 : std::vector<dof_id_type> & ids (recv_ids[id_status.source()]);
1808 129214 : this->comm().receive (id_status.source(), ids, id_tags[blk]);
1809 : #else
1810 : std::vector<dof_id_type> & ids (recv_ids[0]);
1811 : ids = xfer_ids[blk];
1812 : #endif
1813 :
1814 : // note its possible we didn't receive values for objects in
1815 : // this block if they have no components allocated.
1816 14840596 : for (std::size_t idx=0, sz=ids.size(); idx<sz; idx+=2)
1817 : {
1818 : const dof_id_type
1819 14714978 : local_idx = ids[idx+0]-first_object,
1820 14714978 : n_vals_tot_allvecs = ids[idx+1];
1821 :
1822 1338000 : libmesh_assert_less (local_idx, n_objects_blk);
1823 1338000 : libmesh_assert_less (local_idx, obj_val_offsets.size());
1824 :
1825 16052978 : obj_val_offsets[local_idx] = n_vals_tot_allvecs;
1826 : }
1827 :
1828 : #ifdef LIBMESH_HAVE_MPI
1829 : // blocking receive values for this block, imposing no particular order on processor
1830 129214 : Parallel::Status val_status (this->comm().probe (Parallel::any_source, val_tags[blk]));
1831 125618 : std::vector<Number> & vals (recv_vals[val_status.source()]);
1832 129214 : this->comm().receive (val_status.source(), vals, val_tags[blk]);
1833 : #else
1834 : // straight copy without MPI
1835 : std::vector<Number> & vals (recv_vals[0]);
1836 : vals = send_vals[blk];
1837 : #endif
1838 :
1839 129214 : n_val_recvd_blk += vals.size();
1840 : }
1841 :
1842 : // We need the offsets into the output_vals vector for each object.
1843 : // fortunately, this is simply the partial sum of the total number
1844 : // of components for each object
1845 17966 : std::partial_sum(obj_val_offsets.begin(), obj_val_offsets.end(),
1846 : obj_val_offsets.begin());
1847 :
1848 : // wait on any previous asynchronous IO - this *must* complete before
1849 : // we start messing with the output_vals buffer!
1850 19764 : if (async_io.get()) async_io->join();
1851 :
1852 : // this is the actual output buffer that will be written to disk.
1853 : // at ths point we finally know wha size it will be.
1854 19764 : output_vals.resize(n_val_recvd_blk);
1855 :
1856 : // pack data from all processors into output values
1857 145382 : for (auto proc : make_range(this->n_processors()))
1858 : {
1859 125618 : const std::vector<dof_id_type> & ids (recv_ids [proc]);
1860 7192 : const std::vector<Number> & vals(recv_vals[proc]);
1861 7192 : std::vector<Number>::const_iterator proc_vals(vals.begin());
1862 :
1863 14840596 : for (std::size_t idx=0, sz=ids.size(); idx<sz; idx+=2)
1864 : {
1865 : const dof_id_type
1866 14714978 : local_idx = ids[idx+0]-first_object,
1867 16052978 : n_vals_tot_allvecs = ids[idx+1];
1868 :
1869 : // put this object's data into the proper location
1870 : // in the output buffer
1871 13376978 : std::vector<Number>::iterator out_vals(output_vals.begin());
1872 14714978 : if (local_idx != 0)
1873 16042172 : std::advance (out_vals, obj_val_offsets[local_idx-1]);
1874 :
1875 30815663 : for (unsigned int val=0; val<n_vals_tot_allvecs; val++, ++out_vals, ++proc_vals)
1876 : {
1877 1454980 : libmesh_assert (out_vals != output_vals.end());
1878 1454980 : libmesh_assert (proc_vals != vals.end());
1879 16100685 : *out_vals = *proc_vals;
1880 : }
1881 : }
1882 : }
1883 :
1884 : // output_vals buffer is now filled for this block.
1885 : // write it to disk
1886 35932 : async_io = std::make_unique<Threads::Thread>(threaded_io);
1887 21562 : written_length += output_vals.size();
1888 : }
1889 :
1890 : // wait on any previous asynchronous IO - this *must* complete before
1891 : // our stuff goes out of scope
1892 19764 : async_io->join();
1893 32336 : }
1894 :
1895 125618 : Parallel::wait(id_requests);
1896 125618 : Parallel::wait(val_requests);
1897 :
1898 : // we need some synchronization here. Because this method
1899 : // can be called for a range of nodes, then a range of elements,
1900 : // we need some mechanism to prevent processors from racing past
1901 : // to the next range and overtaking ongoing communication. one
1902 : // approach would be to figure out unique tags for each range,
1903 : // but for now we just impose a barrier here. And might as
1904 : // well have it do some useful work.
1905 125618 : this->comm().broadcast(written_length);
1906 :
1907 251236 : return written_length;
1908 118426 : }
1909 :
1910 :
1911 :
1912 0 : unsigned int System::write_SCALAR_dofs (const NumericVector<Number> & vec,
1913 : const unsigned int var,
1914 : Xdr & io) const
1915 : {
1916 0 : unsigned int written_length=0;
1917 0 : std::vector<Number> vals; // The raw values for the local objects in the current block
1918 : // Collect the SCALARs for the current variable
1919 0 : if (this->processor_id() == (this->n_processors()-1))
1920 : {
1921 0 : const DofMap & dof_map = this->get_dof_map();
1922 0 : std::vector<dof_id_type> SCALAR_dofs;
1923 0 : dof_map.SCALAR_dof_indices(SCALAR_dofs, var);
1924 : const unsigned int n_scalar_dofs = cast_int<unsigned int>
1925 0 : (SCALAR_dofs.size());
1926 :
1927 0 : for (unsigned int i=0; i<n_scalar_dofs; i++)
1928 : {
1929 0 : vals.push_back( vec(SCALAR_dofs[i]) );
1930 : }
1931 : }
1932 :
1933 : #ifdef LIBMESH_HAVE_MPI
1934 0 : if (this->n_processors() > 1)
1935 : {
1936 : const Parallel::MessageTag val_tag =
1937 0 : this->comm().get_unique_tag(1);
1938 :
1939 : // Post the receive on processor 0
1940 0 : if (this->processor_id() == 0)
1941 : {
1942 0 : this->comm().receive(this->n_processors()-1, vals, val_tag);
1943 : }
1944 :
1945 : // Send the data to processor 0
1946 0 : if (this->processor_id() == (this->n_processors()-1))
1947 : {
1948 0 : this->comm().send(0, vals, val_tag);
1949 : }
1950 0 : }
1951 : #endif
1952 :
1953 : // -------------------------------------------------------
1954 : // Write the output on processor 0.
1955 0 : if (this->processor_id() == 0)
1956 : {
1957 : const unsigned int vals_size =
1958 0 : cast_int<unsigned int>(vals.size());
1959 0 : io.data_stream (vals.data(), vals_size);
1960 0 : written_length += vals_size;
1961 : }
1962 :
1963 0 : return written_length;
1964 : }
1965 :
1966 :
1967 :
1968 62525 : dof_id_type System::write_serialized_vector (Xdr & io,
1969 : const NumericVector<Number> & vec) const
1970 : {
1971 1790 : parallel_object_only();
1972 :
1973 1790 : libmesh_assert (io.writing());
1974 :
1975 62525 : dof_id_type vec_length = vec.size();
1976 65210 : if (this->processor_id() == 0) io.data (vec_length, "# vector length");
1977 :
1978 1790 : dof_id_type written_length = 0;
1979 :
1980 : //---------------------------------
1981 : // Collect the values for all nodes
1982 1790 : written_length += cast_int<dof_id_type>
1983 123260 : (this->write_serialized_blocked_dof_objects (std::vector<const NumericVector<Number> *>(1,&vec),
1984 62525 : this->get_mesh().n_nodes(),
1985 125050 : this->get_mesh().local_nodes_begin(),
1986 123260 : this->get_mesh().local_nodes_end(),
1987 : io));
1988 :
1989 : //------------------------------------
1990 : // Collect the values for all elements
1991 62525 : written_length += cast_int<dof_id_type>
1992 123260 : (this->write_serialized_blocked_dof_objects (std::vector<const NumericVector<Number> *>(1,&vec),
1993 62525 : this->get_mesh().n_elem(),
1994 125050 : this->get_mesh().local_elements_begin(),
1995 64315 : this->get_mesh().local_elements_end(),
1996 : io));
1997 :
1998 : //-------------------------------------------
1999 : // Finally loop over all the SCALAR variables
2000 125192 : for (auto var : make_range(this->n_vars()))
2001 62667 : if (this->variable(var).type().family == SCALAR)
2002 : {
2003 0 : written_length +=
2004 0 : this->write_SCALAR_dofs (vec, var, io);
2005 : }
2006 :
2007 1790 : if (this->processor_id() == 0)
2008 895 : libmesh_assert_equal_to (written_length, vec_length);
2009 :
2010 62525 : return written_length;
2011 : }
2012 :
2013 :
2014 : template <typename InValType>
2015 284 : std::size_t System::read_serialized_vectors (Xdr & io,
2016 : const std::vector<NumericVector<Number> *> & vectors) const
2017 : {
2018 8 : parallel_object_only();
2019 :
2020 : // Error checking
2021 : // #ifndef NDEBUG
2022 : // // In parallel we better be reading a parallel vector -- if not
2023 : // // we will not set all of its components below!!
2024 : // if (this->n_processors() > 1)
2025 : // {
2026 : // libmesh_assert (vec.type() == PARALLEL ||
2027 : // vec.type() == GHOSTED);
2028 : // }
2029 : // #endif
2030 :
2031 8 : libmesh_assert (io.reading());
2032 :
2033 292 : if (this->processor_id() == 0)
2034 : {
2035 : // sizes
2036 48 : unsigned int num_vecs=0;
2037 48 : dof_id_type vector_length=0;
2038 :
2039 : // Get the number of vectors
2040 48 : io.data(num_vecs);
2041 : // Get the buffer size
2042 48 : io.data(vector_length);
2043 :
2044 52 : libmesh_error_msg_if
2045 : (num_vecs != vectors.size(),
2046 : "Xdr file header declares " << num_vecs << " vectors, but we were asked to read " << vectors.size());
2047 :
2048 48 : if (num_vecs != 0)
2049 : {
2050 48 : libmesh_error_msg_if (vectors[0] == nullptr, "vectors[0] should not be null");
2051 48 : libmesh_error_msg_if (vectors[0]->size() != vector_length, "Inconsistent vector sizes");
2052 : }
2053 : }
2054 :
2055 : // no need to actually communicate these.
2056 : // this->comm().broadcast(num_vecs);
2057 : // this->comm().broadcast(vector_length);
2058 :
2059 : // Cache these - they are not free!
2060 : const dof_id_type
2061 284 : n_nodes = this->get_mesh().n_nodes(),
2062 284 : n_elem = this->get_mesh().n_elem();
2063 :
2064 8 : std::size_t read_length = 0;
2065 :
2066 : //---------------------------------
2067 : // Collect the values for all nodes
2068 24 : read_length +=
2069 536 : this->read_serialized_blocked_dof_objects (n_nodes,
2070 284 : this->get_mesh().local_nodes_begin(),
2071 284 : this->get_mesh().local_nodes_end(),
2072 : InValType(),
2073 : io,
2074 : vectors);
2075 :
2076 : //------------------------------------
2077 : // Collect the values for all elements
2078 284 : read_length +=
2079 536 : this->read_serialized_blocked_dof_objects (n_elem,
2080 284 : this->get_mesh().local_elements_begin(),
2081 284 : this->get_mesh().local_elements_end(),
2082 : InValType(),
2083 : io,
2084 : vectors);
2085 :
2086 : //-------------------------------------------
2087 : // Finally loop over all the SCALAR variables
2088 5240 : for (NumericVector<Number> * vec : vectors)
2089 14172 : for (auto var : make_range(this->n_vars()))
2090 9216 : if (this->variable(var).type().family == SCALAR)
2091 : {
2092 0 : libmesh_assert_not_equal_to (vec, 0);
2093 :
2094 0 : read_length +=
2095 0 : this->read_SCALAR_dofs (var, io, vec);
2096 : }
2097 :
2098 : //---------------------------------------
2099 : // last step - must close all the vectors
2100 5240 : for (NumericVector<Number> * vec : vectors)
2101 : {
2102 140 : libmesh_assert_not_equal_to (vec, 0);
2103 4956 : vec->close();
2104 : }
2105 :
2106 284 : return read_length;
2107 : }
2108 :
2109 :
2110 :
2111 284 : std::size_t System::write_serialized_vectors (Xdr & io,
2112 : const std::vector<const NumericVector<Number> *> & vectors) const
2113 : {
2114 8 : parallel_object_only();
2115 :
2116 8 : libmesh_assert (io.writing());
2117 :
2118 : // Cache these - they are not free!
2119 : const dof_id_type
2120 284 : n_nodes = this->get_mesh().n_nodes(),
2121 284 : n_elem = this->get_mesh().n_elem();
2122 :
2123 8 : std::size_t written_length = 0;
2124 :
2125 292 : if (this->processor_id() == 0)
2126 : {
2127 : unsigned int
2128 48 : n_vec = cast_int<unsigned int>(vectors.size());
2129 : dof_id_type
2130 48 : vec_size = vectors.empty() ? 0 : vectors[0]->size();
2131 : // Set the number of vectors
2132 48 : io.data(n_vec, "# number of vectors");
2133 : // Set the buffer size
2134 48 : io.data(vec_size, "# vector length");
2135 : }
2136 :
2137 : //---------------------------------
2138 : // Collect the values for all nodes
2139 8 : written_length +=
2140 284 : this->write_serialized_blocked_dof_objects (vectors,
2141 : n_nodes,
2142 568 : this->get_mesh().local_nodes_begin(),
2143 560 : this->get_mesh().local_nodes_end(),
2144 : io);
2145 :
2146 : //------------------------------------
2147 : // Collect the values for all elements
2148 284 : written_length +=
2149 284 : this->write_serialized_blocked_dof_objects (vectors,
2150 : n_elem,
2151 568 : this->get_mesh().local_elements_begin(),
2152 560 : this->get_mesh().local_elements_end(),
2153 : io);
2154 :
2155 : //-------------------------------------------
2156 : // Finally loop over all the SCALAR variables
2157 5240 : for (const NumericVector<Number> * vec : vectors)
2158 14172 : for (auto var : make_range(this->n_vars()))
2159 9216 : if (this->variable(var).type().family == SCALAR)
2160 : {
2161 0 : libmesh_assert_not_equal_to (vec, 0);
2162 :
2163 0 : written_length +=
2164 0 : this->write_SCALAR_dofs (*vec, var, io);
2165 : }
2166 :
2167 284 : return written_length;
2168 : }
2169 :
2170 :
2171 :
2172 :
2173 : template LIBMESH_EXPORT void System::read_parallel_data<Number> (Xdr & io, const bool read_additional_data);
2174 : template LIBMESH_EXPORT void System::read_serialized_data<Number> (Xdr & io, const bool read_additional_data);
2175 : template LIBMESH_EXPORT numeric_index_type System::read_serialized_vector<Number> (Xdr & io, NumericVector<Number> * vec);
2176 : template LIBMESH_EXPORT std::size_t System::read_serialized_vectors<Number> (Xdr & io, const std::vector<NumericVector<Number> *> & vectors) const;
2177 : #ifdef LIBMESH_USE_COMPLEX_NUMBERS
2178 : template LIBMESH_EXPORT void System::read_parallel_data<Real> (Xdr & io, const bool read_additional_data);
2179 : template LIBMESH_EXPORT void System::read_serialized_data<Real> (Xdr & io, const bool read_additional_data);
2180 : template LIBMESH_EXPORT numeric_index_type System::read_serialized_vector<Real> (Xdr & io, NumericVector<Number> * vec);
2181 : template LIBMESH_EXPORT std::size_t System::read_serialized_vectors<Real> (Xdr & io, const std::vector<NumericVector<Number> *> & vectors) const;
2182 : #endif
2183 :
2184 : } // namespace libMesh
|
