ryujin 2.1.1 revision 350e54cc11f3d21282bddcf3e6517944dbc508bf
quantities.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include "openmp.h"
9#include "quantities.h"
10
11#include <deal.II/base/function_parser.h>
12#include <deal.II/base/mpi.templates.h>
13#include <deal.II/base/work_stream.h>
14#include <deal.II/dofs/dof_tools.h>
15
16#include <fstream>
17
18DEAL_II_NAMESPACE_OPEN
19template <int rank, int dim, typename Number>
20bool operator<(const Tensor<rank, dim, Number> &left,
21 const Tensor<rank, dim, Number> &right)
22{
23 return std::lexicographical_compare(
24 left.begin_raw(), left.end_raw(), right.begin_raw(), right.end_raw());
25}
26DEAL_II_NAMESPACE_CLOSE
27
28namespace ryujin
29{
30 using namespace dealii;
31
32 namespace
33 {
34 template <typename T>
35 const std::string &get_options_from_name(const T &manifolds,
36 const std::string &name)
37 {
38 const auto it =
39 std::find_if(manifolds.begin(),
40 manifolds.end(),
41 [&, name = std::cref(name)](const auto &element) {
42 return std::get<0>(element) == name.get();
43 });
44 Assert(it != manifolds.end(), dealii::ExcInternalError());
45 return std::get<2>(*it);
46 }
47 } // namespace
48
49
50 template <typename Description, int dim, typename Number>
52 const MPIEnsemble &mpi_ensemble,
53 const OfflineData<dim, Number> &offline_data,
54 const HyperbolicSystem &hyperbolic_system,
55 const ParabolicSystem &parabolic_system,
56 const std::string &subsection /*= "Quantities"*/)
57 : ParameterAcceptor(subsection)
58 , mpi_ensemble_(mpi_ensemble)
59 , offline_data_(&offline_data)
60 , hyperbolic_system_(&hyperbolic_system)
61 , parabolic_system_(&parabolic_system)
62 , base_name_("")
63 , mesh_files_have_been_written_(false)
64 {
65 add_parameter("interior manifolds",
66 interior_manifolds_,
67 "List of level set functions describing interior manifolds. "
68 "The description is used to only output point values for "
69 "vertices belonging to a certain level set. "
70 "Format: '<name> : <level set formula> : <options> , [...] "
71 "(options: time_averaged, space_averaged, instantaneous)");
72
73 add_parameter("boundary manifolds",
74 boundary_manifolds_,
75 "List of level set functions describing boundary. The "
76 "description is used to only output point values for "
77 "boundary vertices belonging to a certain level set. "
78 "Format: '<name> : <level set formula> : <options> , [...] "
79 "(options: time_averaged, space_averaged, instantaneous)");
81 clear_temporal_statistics_on_writeout_ = true;
82 add_parameter("clear statistics on writeout",
83 clear_temporal_statistics_on_writeout_,
84 "If set to true then all temporal statistics (for "
85 "\"time_averaged\" quantities) accumulated so far are reset "
86 "each time a writeout of quantities is performed");
87 }
88
89
90 template <typename Description, int dim, typename Number>
92 {
93#ifdef DEBUG_OUTPUT
94 std::cout << "Quantities<dim, Number>::prepare()" << std::endl;
95#endif
96
97 base_name_ = name;
98
99 /* Force to write to a new time series file: */
100 time_series_cycle_.reset();
101
102 const unsigned int n_owned = offline_data_->n_locally_owned();
103 const auto &sparsity_simd = offline_data_->sparsity_pattern_simd();
104
105 /*
106 * Create interior maps and allocate statistics.
107 *
108 * We have to loop over the cells and populate the std::map interior_maps_.
109 */
110
111 interior_maps_.clear();
112 std::transform(
113 interior_manifolds_.begin(),
114 interior_manifolds_.end(),
115 std::inserter(interior_maps_, interior_maps_.end()),
116 [this, n_owned, &sparsity_simd](auto it) {
117 const auto &[name, expression, option] = it;
118 FunctionParser<dim> level_set_function(expression);
119
120 std::vector<interior_point> map;
121 std::map<int, interior_point> preliminary_map;
122
123 const auto &discretization = offline_data_->discretization();
124 const auto &dof_handler = offline_data_->dof_handler();
125
126 const auto support_points =
127 dof_handler.get_fe().get_unit_support_points();
128
129 std::vector<dealii::types::global_dof_index> local_dof_indices;
130
131 /* Loop over cells */
132 for (auto cell : dof_handler.active_cell_iterators()) {
133
134 /* skip if not locally owned */
135 if (!cell->is_locally_owned())
136 continue;
137
138 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
139 local_dof_indices.resize(dofs_per_cell);
140 cell->get_active_or_mg_dof_indices(local_dof_indices);
141
142 const auto &mapping =
143 discretization.mapping()[cell->active_fe_index()];
144
145 for (unsigned int j = 0; j < dofs_per_cell; ++j) {
146
147 Point<dim> position =
148 mapping.transform_unit_to_real_cell(cell, support_points[j]);
149
150 /*
151 * Insert index, interior mass value and position into
152 * a preliminary map if we satisfy level set condition.
153 */
154
155 if (std::abs(level_set_function.value(position)) > 1.e-12)
156 continue;
157
158 const auto global_index = local_dof_indices[j];
159 const auto index =
160 offline_data_->scalar_partitioner()->global_to_local(
161 global_index);
162
163 /* Skip constrained degrees of freedom: */
164 const unsigned int row_length = sparsity_simd.row_length(index);
165 if (row_length == 1)
166 continue;
167
168 if (index >= n_owned)
169 continue;
170
171 const Number interior_mass =
172 offline_data_->lumped_mass_matrix().local_element(index);
173 // FIXME: change to std::set
174 preliminary_map[index] = {index, interior_mass, position};
175 }
176 }
177
178 /*
179 * Now we populate the std::vector(interior_point) object called map.
180 */
181 // FIXME: use std::copy
182 for (const auto &[index, tuple] : preliminary_map) {
183 map.push_back(tuple);
184 }
185
186 return std::make_pair(name, map);
187 });
188
189 /*
190 * Create boundary maps and allocate statistics vector:
191 *
192 * We want to loop over the boundary_map() once and populate the map
193 * object boundary_maps_. We have to create a vector of
194 * boundary_manifolds.size() that holds a std::vector<boundary_point>
195 * for each map entry.
196 */
197
198 boundary_maps_.clear();
199 std::transform(
200 boundary_manifolds_.begin(),
201 boundary_manifolds_.end(),
202 std::inserter(boundary_maps_, boundary_maps_.end()),
203 [this, n_owned](auto it) {
204 const auto &[name, expression, option] = it;
205 FunctionParser<dim> level_set_function(expression);
206
207 std::vector<boundary_point> map;
208
209 for (const auto &entry : offline_data_->boundary_map()) {
210 // [i, normal, normal_mass, boundary_mass, id, position] = entry
211 const auto &i = std::get<0>(entry);
212
213 /* skip nonlocal */
214 if (i >= n_owned)
215 continue;
216
217 /* skip constrained */
218 if (offline_data_->affine_constraints().is_constrained(
219 offline_data_->scalar_partitioner()->local_to_global(i)))
220 continue;
221
222 const auto &position = std::get<5>(entry);
223 if (std::abs(level_set_function.value(position)) < 1.e-12)
224 map.push_back(entry);
225 }
226 return std::make_pair(name, map);
227 });
228
229 /* Clear statistics: */
230 clear_statistics();
231
232 /* Make sure we output new mesh files: */
233 mesh_files_have_been_written_ = false;
234
235 /* Prepare header string: */
236 const auto &names = View::primitive_component_names;
237 header_ = std::accumulate(
238 std::begin(names),
239 std::end(names),
240 std::string(),
241 [](const std::string &description, const std::string &name) {
242 return description.empty()
243 ? (std::string("primitive state (") + name)
244 : (description + ", " + name);
245 }) +
246 ")\t and 2nd moments\n";
247 }
248
249
250 template <typename Description, int dim, typename Number>
251 void
252 Quantities<Description, dim, Number>::write_mesh_files(unsigned int cycle)
253 {
254 /*
255 * Output interior maps:
256 */
257
258 for (const auto &[name, interior_map] : interior_maps_) {
259 /* Skip outputting the boundary map for spatial averages. */
260 const auto &options = get_options_from_name(interior_manifolds_, name);
261 if (options.find("instantaneous") == std::string::npos &&
262 options.find("time_averaged") == std::string::npos)
263 continue;
264
265 /*
266 * FIXME: This currently distributes boundary maps to all MPI ranks.
267 * This is unnecessarily wasteful. Ideally, we should do MPI IO with
268 * only MPI ranks participating who actually have boundary values.
269 */
270
271 const auto received = Utilities::MPI::gather(
272 mpi_ensemble_.ensemble_communicator(), interior_map);
273
274 if (Utilities::MPI::this_mpi_process(
275 mpi_ensemble_.ensemble_communicator()) == 0) {
276
277 std::ofstream output(base_name_ + "-" + name + "-R" +
278 Utilities::to_string(cycle, 4) + "-points.dat");
279
280 output << std::scientific << std::setprecision(14);
281
282 output << "#\n# position\tinterior mass\n";
283
284 unsigned int rank = 0;
285 for (const auto &entries : received) {
286 output << "# rank " << rank++ << "\n";
287 for (const auto &entry : entries) {
288 const auto &[index, mass_i, x_i] = entry;
289 output << x_i << "\t" << mass_i << "\n";
290 } /*entry*/
291 } /*entries*/
292
293 output << std::flush;
294 }
295 }
296
297 /*
298 * Output boundary maps:
299 */
300
301 for (const auto &[name, boundary_map] : boundary_maps_) {
302 /* Skip outputting the boundary map for spatial averages. */
303 const auto &options = get_options_from_name(boundary_manifolds_, name);
304 if (options.find("instantaneous") == std::string::npos &&
305 options.find("time_averaged") == std::string::npos)
306 continue;
307
308 /*
309 * FIXME: This currently distributes boundary maps to all MPI ranks.
310 * This is unnecessarily wasteful. Ideally, we should do MPI IO with
311 * only MPI ranks participating who actually have boundary values.
312 */
313
314 const auto received = Utilities::MPI::gather(
315 mpi_ensemble_.ensemble_communicator(), boundary_map);
316
317 if (Utilities::MPI::this_mpi_process(
318 mpi_ensemble_.ensemble_communicator()) == 0) {
319
320 std::ofstream output(base_name_ + "-" + name + "-R" +
321 Utilities::to_string(cycle, 4) + "-points.dat");
322
323 output << std::scientific << std::setprecision(14);
324
325 output << "#\n# position\tnormal\tnormal mass\tboundary mass\n";
326
327 unsigned int rank = 0;
328 for (const auto &entries : received) {
329 output << "# rank " << rank++ << "\n";
330 for (const auto &entry : entries) {
331 const auto &[index, n_i, nm_i, bm_i, id, x_i] = entry;
332 output << x_i << "\t" << n_i << "\t" << nm_i << "\t" << bm_i
333 << "\n";
334 } /*entry*/
335 } /*entries*/
336
337 output << std::flush;
338 }
339 }
340 }
341
342
343 template <typename Description, int dim, typename Number>
344 void Quantities<Description, dim, Number>::clear_statistics()
345 {
346 const auto reset = [](const auto &manifold_map, auto &statistics_map) {
347 for (const auto &[name, data_map] : manifold_map) {
348 const auto n_entries = data_map.size();
349 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
350 statistics_map[name];
351 val_old.resize(n_entries);
352 val_new.resize(n_entries);
353 val_sum.resize(n_entries);
354 t_old = t_new = t_sum = 0.;
355 }
356 };
357
358 /* Clear statistics and time series: */
359
360 interior_statistics_.clear();
361 reset(interior_maps_, interior_statistics_);
362 interior_time_series_.clear();
363
364 boundary_statistics_.clear();
365 reset(boundary_maps_, boundary_statistics_);
366 boundary_time_series_.clear();
367 }
368
369
370 template <typename Description, int dim, typename Number>
371 template <typename point_type, typename value_type>
372 value_type Quantities<Description, dim, Number>::internal_accumulate(
373 const StateVector &state_vector,
374 const std::vector<point_type> &points_vector,
375 std::vector<value_type> &val_new)
376 {
377 const auto &U = std::get<0>(state_vector);
378
379 value_type spatial_average;
380 Number mass_sum = Number(0.);
381
382 std::transform(
383 points_vector.begin(),
384 points_vector.end(),
385 val_new.begin(),
386 [&](auto point) -> value_type {
387 const auto i = std::get<0>(point);
388 /*
389 * Small trick to get the correct index for retrieving the
390 * boundary mass.
391 */
392 constexpr auto index =
393 std::is_same_v<point_type, interior_point> ? 1 : 3;
394 const auto mass_i = std::get<index>(point);
395
396 const auto U_i = U.get_tensor(i);
397 const auto view = hyperbolic_system_->template view<dim, Number>();
398 const auto primitive_state = view.to_primitive_state(U_i);
399
400 value_type result;
401 std::get<0>(result) = primitive_state;
402 /* Compute second moments of the primitive state: */
403 std::get<1>(result) = schur_product(primitive_state, primitive_state);
404
405 mass_sum += mass_i;
406 std::get<0>(spatial_average) += mass_i * std::get<0>(result);
407 std::get<1>(spatial_average) += mass_i * std::get<1>(result);
408
409 return result;
410 });
411
412 /* synchronize MPI ranks (MPI Barrier): */
413
414 mass_sum =
415 Utilities::MPI::sum(mass_sum, mpi_ensemble_.ensemble_communicator());
416
417 std::get<0>(spatial_average) = Utilities::MPI::sum(
418 std::get<0>(spatial_average), mpi_ensemble_.ensemble_communicator());
419 std::get<1>(spatial_average) = Utilities::MPI::sum(
420 std::get<1>(spatial_average), mpi_ensemble_.ensemble_communicator());
421
422 /* take average: */
423
424 std::get<0>(spatial_average) /= mass_sum;
425 std::get<1>(spatial_average) /= mass_sum;
426
427 return spatial_average;
428 }
429
430
431 template <typename Description, int dim, typename Number>
432 template <typename value_type>
433 void Quantities<Description, dim, Number>::internal_write_out(
434 const std::string &file_name,
435 const std::string &time_stamp,
436 const std::vector<value_type> &values,
437 const Number scale)
438 {
439 /*
440 * FIXME: This currently distributes interior maps to all MPI ranks.
441 * This is unnecessarily wasteful. Ideally, we should do MPI IO with
442 * only MPI ranks participating who actually have interior values.
443 */
444
445 const auto received =
446 Utilities::MPI::gather(mpi_ensemble_.ensemble_communicator(), values);
447
448 if (Utilities::MPI::this_mpi_process(
449 mpi_ensemble_.ensemble_communicator()) == 0) {
450
451 std::ofstream output(file_name);
452 output << std::scientific << std::setprecision(14);
453 output << time_stamp << "# " << header_;
454
455 unsigned int rank = 0;
456 for (const auto &entries : received) {
457 output << "# rank " << rank++ << "\n";
458 for (const auto &entry : entries) {
459 const auto &[state, state_square] = entry;
460 output << scale * state << "\t" << scale * state_square << "\n";
461 } /*entry*/
462 } /*entries*/
463
464 output << std::flush;
465 }
466 }
467
468
469 template <typename Description, int dim, typename Number>
470 template <typename value_type>
471 void Quantities<Description, dim, Number>::internal_write_out_time_series(
472 const std::string &file_name,
473 const std::vector<std::tuple<Number, value_type>> &values,
474 bool append)
475 {
476 if (Utilities::MPI::this_mpi_process(
477 mpi_ensemble_.ensemble_communicator()) == 0) {
478 std::ofstream output;
479 output << std::scientific << std::setprecision(14);
480
481 if (append) {
482 output.open(file_name, std::ofstream::out | std::ofstream::app);
483 } else {
484 output.open(file_name, std::ofstream::out | std::ofstream::trunc);
485 output << "# time t\t" << header_;
486 }
487
488 for (const auto &entry : values) {
489 const auto t = std::get<0>(entry);
490 const auto &[state, state_square] = std::get<1>(entry);
491
492 output << t << "\t" << state << "\t" << state_square << "\n";
493 }
494
495 output << std::flush;
496 output.close();
497 }
498 }
499
500
501 template <typename Description, int dim, typename Number>
503 const StateVector &state_vector, const Number t)
504 {
505#ifdef DEBUG_OUTPUT
506 std::cout << "Quantities<dim, Number>::accumulate()" << std::endl;
507#endif
508
509 const auto accumulate = [&](const auto &point_maps,
510 const auto &manifolds,
511 auto &statistics,
512 auto &time_series) {
513 for (const auto &[name, point_map] : point_maps) {
514
515 /* Find the correct option string in manifolds */
516 const auto &options = get_options_from_name(manifolds, name);
517
518 /* skip if we don't average in space or time: */
519 if (options.find("time_averaged") == std::string::npos &&
520 options.find("space_averaged") == std::string::npos)
521 continue;
522
523 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
524 statistics[name];
525
526 std::swap(t_old, t_new);
527 std::swap(val_old, val_new);
528
529 /* accumulate new values */
530
531 const auto spatial_average =
532 internal_accumulate(state_vector, point_map, val_new);
533
534 /* Average in time with trapezoidal rule: */
535
536 if (RYUJIN_UNLIKELY(t_old == Number(0.) && t_new == Number(0.))) {
537 /* We have not accumulated any statistics yet: */
538 t_old = t - 1.;
539 t_new = t;
540
541 } else {
542
543 t_new = t;
544 const Number tau = t_new - t_old;
545
546 for (std::size_t i = 0; i < val_sum.size(); ++i) {
547 std::get<0>(val_sum[i]) += 0.5 * tau * std::get<0>(val_old[i]);
548 std::get<0>(val_sum[i]) += 0.5 * tau * std::get<0>(val_new[i]);
549 std::get<1>(val_sum[i]) += 0.5 * tau * std::get<1>(val_old[i]);
550 std::get<1>(val_sum[i]) += 0.5 * tau * std::get<1>(val_new[i]);
551 }
552 t_sum += tau;
553 }
554
555 /* Record average in space: */
556 time_series[name].push_back({t, spatial_average});
557 }
558 };
559
560 accumulate(interior_maps_,
561 interior_manifolds_,
562 interior_statistics_,
563 interior_time_series_);
564
565 accumulate(boundary_maps_,
566 boundary_manifolds_,
567 boundary_statistics_,
568 boundary_time_series_);
569 }
570
571
572 template <typename Description, int dim, typename Number>
574 const StateVector &state_vector, const Number t, unsigned int cycle)
575 {
576#ifdef DEBUG_OUTPUT
577 std::cout << "Quantities<dim, Number>::write_out()" << std::endl;
578#endif
579
580 /*
581 * First, write out mesh files if this hasn't happened yet.
582 */
583 if (!mesh_files_have_been_written_) {
584 write_mesh_files(cycle);
585 mesh_files_have_been_written_ = true;
586 }
587
588 /*
589 * Next write out instantaneous and time_averaged maps, and flush the
590 * space_averaged values to the corresponding log files:
591 */
592
593 const auto write_out = [&](const auto &point_maps,
594 const auto &manifolds,
595 auto &statistics,
596 auto &time_series) {
597 for (const auto &[name, point_map] : point_maps) {
598
599 /* Find the correct option string in manifolds */
600 const auto &options = get_options_from_name(manifolds, name);
601
602 const auto prefix =
603 base_name_ + "-" + name + "-R" + Utilities::to_string(cycle, 4);
604
605 /*
606 * Compute and output instantaneous field:
607 */
608
609 if (options.find("instantaneous") != std::string::npos) {
610
611 const std::string file_name = prefix + "-instantaneous.dat";
612
613 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
614 statistics[name];
615
616 std::stringstream time_stamp;
617 time_stamp << std::scientific << std::setprecision(14);
618 time_stamp << "# at t = " << t << std::endl;
619
620 /* We have not computed any updated statistics yet: */
621
622 if (options.find("time_averaged") == std::string::npos &&
623 options.find("space_averaged") == std::string::npos)
624 internal_accumulate(state_vector, point_map, val_new);
625 else
626 AssertThrow(t_new == t, dealii::ExcInternalError());
627
628 internal_write_out(file_name, time_stamp.str(), val_new, Number(1.));
629 }
630
631 /*
632 * Output time averaged field:
633 */
634
635 if (options.find("time_averaged") != std::string::npos) {
636
637 const std::string file_name = prefix + "-time_averaged.dat";
638
639 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
640 statistics[name];
641
642 /* Check whether we have accumulated any statistics yet: */
643 if (t_sum != Number(0.)) {
644 std::stringstream time_stamp;
645 time_stamp << std::scientific << std::setprecision(14);
646 time_stamp << "# averaged from t = " << t_new - t_sum
647 << " to t = " << t_new << std::endl;
648
649 internal_write_out(
650 file_name, time_stamp.str(), val_sum, Number(1.) / t_sum);
651 }
652 }
653
654 /*
655 * Output space averaged field:
656 */
657
658 if (options.find("space_averaged") != std::string::npos) {
659 bool append = true;
660 if (!time_series_cycle_.has_value()) {
661 time_series_cycle_ = cycle;
662 append = false;
663 }
664
665 const auto file_name =
666 base_name_ + "-" + name + "-R" +
667 Utilities::to_string(time_series_cycle_.value(), 4) +
668 "-space_averaged_time_series.dat";
669
670 auto &series = time_series[name];
671 internal_write_out_time_series(file_name, series, /*append*/ append);
672 series.clear();
673 }
674 }
675 };
676
677 write_out(interior_maps_,
678 interior_manifolds_,
679 interior_statistics_,
680 interior_time_series_);
681
682 write_out(boundary_maps_,
683 boundary_manifolds_,
684 boundary_statistics_,
685 boundary_time_series_);
686
687 if (clear_temporal_statistics_on_writeout_)
688 clear_statistics();
689 }
690
691} /* namespace ryujin */
typename Description::HyperbolicSystem HyperbolicSystem
Definition: quantities.h:38
typename Description::ParabolicSystem ParabolicSystem
Definition: quantities.h:39
typename View::StateVector StateVector
Definition: quantities.h:46
Quantities(const MPIEnsemble &mpi_ensemble, const OfflineData< dim, Number > &offline_data, const HyperbolicSystem &hyperbolic_system, const ParabolicSystem &parabolic_system, const std::string &subsection="/Quantities")
#define RYUJIN_UNLIKELY(x)
Definition: openmp.h:132
std::tuple< MultiComponentVector< Number, problem_dim >, MultiComponentVector< Number, prec_dim >, BlockVector< Number > > StateVector
Definition: state_vector.h:51
DEAL_II_NAMESPACE_OPEN bool operator<(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, Number > &right)