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>
19template <
int rank,
int dim,
typename Number>
20bool operator<(
const Tensor<rank, dim, Number> &left,
21 const Tensor<rank, dim, Number> &right)
23 return std::lexicographical_compare(
24 left.begin_raw(), left.end_raw(), right.begin_raw(), right.end_raw());
26DEAL_II_NAMESPACE_CLOSE
30 using namespace dealii;
35 const std::string &get_options_from_name(
const T &manifolds,
36 const std::string &name)
39 std::find_if(manifolds.begin(),
41 [&, name = std::cref(name)](
const auto &element) {
42 return std::get<0>(element) == name.get();
44 Assert(it != manifolds.end(), dealii::ExcInternalError());
45 return std::get<2>(*it);
50 template <
typename Description,
int dim,
typename Number>
56 const std::string &subsection )
57 : ParameterAcceptor(subsection)
58 , mpi_ensemble_(mpi_ensemble)
59 , offline_data_(&offline_data)
60 , hyperbolic_system_(&hyperbolic_system)
61 , parabolic_system_(¶bolic_system)
63 , mesh_files_have_been_written_(false)
65 add_parameter(
"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)");
73 add_parameter(
"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");
90 template <
typename Description,
int dim,
typename Number>
94 std::cout <<
"Quantities<dim, Number>::prepare()" << std::endl;
100 time_series_cycle_.reset();
102 const unsigned int n_owned = offline_data_->n_locally_owned();
103 const auto &sparsity_simd = offline_data_->sparsity_pattern_simd();
111 interior_maps_.clear();
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);
120 std::vector<interior_point> map;
121 std::map<int, interior_point> preliminary_map;
123 const auto &discretization = offline_data_->discretization();
124 const auto &dof_handler = offline_data_->dof_handler();
126 const auto support_points =
127 dof_handler.get_fe().get_unit_support_points();
129 std::vector<dealii::types::global_dof_index> local_dof_indices;
132 for (auto cell : dof_handler.active_cell_iterators()) {
135 if (!cell->is_locally_owned())
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);
142 const auto &mapping =
143 discretization.mapping()[cell->active_fe_index()];
145 for (unsigned int j = 0; j < dofs_per_cell; ++j) {
147 Point<dim> position =
148 mapping.transform_unit_to_real_cell(cell, support_points[j]);
155 if (std::abs(level_set_function.value(position)) > 1.e-12)
158 const auto global_index = local_dof_indices[j];
160 offline_data_->scalar_partitioner()->global_to_local(
164 const unsigned int row_length = sparsity_simd.row_length(index);
168 if (index >= n_owned)
171 const Number interior_mass =
172 offline_data_->lumped_mass_matrix().local_element(index);
174 preliminary_map[index] = {index, interior_mass, position};
182 for (
const auto &[index, tuple] : preliminary_map) {
183 map.push_back(tuple);
186 return std::make_pair(name, map);
198 boundary_maps_.clear();
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);
207 std::vector<boundary_point> map;
209 for (const auto &entry : offline_data_->boundary_map()) {
211 const auto &i = std::get<0>(entry);
218 if (offline_data_->affine_constraints().is_constrained(
219 offline_data_->scalar_partitioner()->local_to_global(i)))
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);
226 return std::make_pair(name, map);
233 mesh_files_have_been_written_ =
false;
236 const auto &names = View::primitive_component_names;
237 header_ = std::accumulate(
241 [](
const std::string &description,
const std::string &name) {
242 return description.empty()
243 ? (std::string(
"primitive state (") + name)
244 : (description +
", " + name);
246 ")\t and 2nd moments\n";
250 template <
typename Description,
int dim,
typename Number>
252 Quantities<Description, dim, Number>::write_mesh_files(
unsigned int cycle)
258 for (
const auto &[name, interior_map] : interior_maps_) {
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)
271 const auto received = Utilities::MPI::gather(
272 mpi_ensemble_.ensemble_communicator(), interior_map);
274 if (Utilities::MPI::this_mpi_process(
275 mpi_ensemble_.ensemble_communicator()) == 0) {
277 std::ofstream output(base_name_ +
"-" + name +
"-R" +
278 Utilities::to_string(cycle, 4) +
"-points.dat");
280 output << std::scientific << std::setprecision(14);
282 output <<
"#\n# position\tinterior mass\n";
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";
293 output << std::flush;
301 for (
const auto &[name, boundary_map] : boundary_maps_) {
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)
314 const auto received = Utilities::MPI::gather(
315 mpi_ensemble_.ensemble_communicator(), boundary_map);
317 if (Utilities::MPI::this_mpi_process(
318 mpi_ensemble_.ensemble_communicator()) == 0) {
320 std::ofstream output(base_name_ +
"-" + name +
"-R" +
321 Utilities::to_string(cycle, 4) +
"-points.dat");
323 output << std::scientific << std::setprecision(14);
325 output <<
"#\n# position\tnormal\tnormal mass\tboundary mass\n";
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
337 output << std::flush;
343 template <
typename Description,
int dim,
typename Number>
344 void Quantities<Description, dim, Number>::clear_statistics()
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.;
360 interior_statistics_.clear();
361 reset(interior_maps_, interior_statistics_);
362 interior_time_series_.clear();
364 boundary_statistics_.clear();
365 reset(boundary_maps_, boundary_statistics_);
366 boundary_time_series_.clear();
370 template <
typename Description,
int dim,
typename Number>
371 template <
typename po
int_type,
typename value_type>
372 value_type Quantities<Description, dim, Number>::internal_accumulate(
374 const std::vector<point_type> &points_vector,
375 std::vector<value_type> &val_new)
377 const auto &U = std::get<0>(state_vector);
379 value_type spatial_average;
380 Number mass_sum = Number(0.);
383 points_vector.begin(),
386 [&](
auto point) -> value_type {
387 const auto i = std::get<0>(point);
392 constexpr auto index =
393 std::is_same_v<point_type, interior_point> ? 1 : 3;
394 const auto mass_i = std::get<index>(point);
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);
401 std::get<0>(result) = primitive_state;
403 std::get<1>(result) = schur_product(primitive_state, primitive_state);
406 std::get<0>(spatial_average) += mass_i * std::get<0>(result);
407 std::get<1>(spatial_average) += mass_i * std::get<1>(result);
415 Utilities::MPI::sum(mass_sum, mpi_ensemble_.ensemble_communicator());
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());
424 std::get<0>(spatial_average) /= mass_sum;
425 std::get<1>(spatial_average) /= mass_sum;
427 return spatial_average;
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,
445 const auto received =
446 Utilities::MPI::gather(mpi_ensemble_.ensemble_communicator(), values);
448 if (Utilities::MPI::this_mpi_process(
449 mpi_ensemble_.ensemble_communicator()) == 0) {
451 std::ofstream output(file_name);
452 output << std::scientific << std::setprecision(14);
453 output << time_stamp <<
"# " << header_;
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";
464 output << std::flush;
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,
476 if (Utilities::MPI::this_mpi_process(
477 mpi_ensemble_.ensemble_communicator()) == 0) {
478 std::ofstream output;
479 output << std::scientific << std::setprecision(14);
482 output.open(file_name, std::ofstream::out | std::ofstream::app);
484 output.open(file_name, std::ofstream::out | std::ofstream::trunc);
485 output <<
"# time t\t" << header_;
488 for (
const auto &entry : values) {
489 const auto t = std::get<0>(entry);
490 const auto &[state, state_square] = std::get<1>(entry);
492 output << t <<
"\t" << state <<
"\t" << state_square <<
"\n";
495 output << std::flush;
501 template <
typename Description,
int dim,
typename Number>
506 std::cout <<
"Quantities<dim, Number>::accumulate()" << std::endl;
509 const auto accumulate = [&](
const auto &point_maps,
510 const auto &manifolds,
513 for (
const auto &[name, point_map] : point_maps) {
516 const auto &options = get_options_from_name(manifolds, name);
519 if (options.find(
"time_averaged") == std::string::npos &&
520 options.find(
"space_averaged") == std::string::npos)
523 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
526 std::swap(t_old, t_new);
527 std::swap(val_old, val_new);
531 const auto spatial_average =
532 internal_accumulate(state_vector, point_map, val_new);
544 const Number tau = t_new - t_old;
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]);
556 time_series[name].push_back({t, spatial_average});
560 accumulate(interior_maps_,
562 interior_statistics_,
563 interior_time_series_);
565 accumulate(boundary_maps_,
567 boundary_statistics_,
568 boundary_time_series_);
572 template <
typename Description,
int dim,
typename Number>
574 const StateVector &state_vector,
const Number t,
unsigned int cycle)
577 std::cout <<
"Quantities<dim, Number>::write_out()" << std::endl;
583 if (!mesh_files_have_been_written_) {
584 write_mesh_files(cycle);
585 mesh_files_have_been_written_ =
true;
593 const auto write_out = [&](
const auto &point_maps,
594 const auto &manifolds,
597 for (
const auto &[name, point_map] : point_maps) {
600 const auto &options = get_options_from_name(manifolds, name);
603 base_name_ +
"-" + name +
"-R" + Utilities::to_string(cycle, 4);
609 if (options.find(
"instantaneous") != std::string::npos) {
611 const std::string file_name = prefix +
"-instantaneous.dat";
613 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
616 std::stringstream time_stamp;
617 time_stamp << std::scientific << std::setprecision(14);
618 time_stamp <<
"# at t = " << t << std::endl;
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);
626 AssertThrow(t_new == t, dealii::ExcInternalError());
628 internal_write_out(file_name, time_stamp.str(), val_new, Number(1.));
635 if (options.find(
"time_averaged") != std::string::npos) {
637 const std::string file_name = prefix +
"-time_averaged.dat";
639 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
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;
650 file_name, time_stamp.str(), val_sum, Number(1.) / t_sum);
658 if (options.find(
"space_averaged") != std::string::npos) {
660 if (!time_series_cycle_.has_value()) {
661 time_series_cycle_ = cycle;
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";
670 auto &series = time_series[name];
671 internal_write_out_time_series(file_name, series, append);
677 write_out(interior_maps_,
679 interior_statistics_,
680 interior_time_series_);
682 write_out(boundary_maps_,
684 boundary_statistics_,
685 boundary_time_series_);
687 if (clear_temporal_statistics_on_writeout_)
typename Description::HyperbolicSystem HyperbolicSystem
typename Description::ParabolicSystem ParabolicSystem
typename View::StateVector StateVector
Quantities(const MPIEnsemble &mpi_ensemble, const OfflineData< dim, Number > &offline_data, const HyperbolicSystem &hyperbolic_system, const ParabolicSystem ¶bolic_system, const std::string &subsection="/Quantities")
#define RYUJIN_UNLIKELY(x)
std::tuple< MultiComponentVector< Number, problem_dim >, MultiComponentVector< Number, prec_dim >, BlockVector< Number > > StateVector
DEAL_II_NAMESPACE_OPEN bool operator<(const Tensor< rank, dim, Number > &left, const Tensor< rank, dim, Number > &right)