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>
52 const MPI_Comm &mpi_communicator,
55 const std::string &subsection )
56 : ParameterAcceptor(subsection)
57 , mpi_communicator_(mpi_communicator)
58 , hyperbolic_system_(&hyperbolic_system)
59 , offline_data_(&offline_data)
61 , time_series_cycle_(1)
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>
95 std::cout <<
"Quantities<dim, Number>::prepare()" << std::endl;
99 time_series_cycle_ = cycle;
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 unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
128 const auto support_points =
129 dof_handler.get_fe().get_unit_support_points();
131 std::vector<dealii::types::global_dof_index> local_dof_indices(
135 for (auto cell : dof_handler.active_cell_iterators()) {
138 if (!cell->is_locally_owned())
141 cell->get_active_or_mg_dof_indices(local_dof_indices);
143 for (unsigned int j = 0; j < dofs_per_cell; ++j) {
145 Point<dim> position =
146 discretization.mapping().transform_unit_to_real_cell(
147 cell, support_points[j]);
154 if (std::abs(level_set_function.value(position)) > 1.e-12)
157 const auto global_index = local_dof_indices[j];
159 offline_data_->scalar_partitioner()->global_to_local(
163 const unsigned int row_length = sparsity_simd.row_length(index);
167 if (index >= n_owned)
170 const Number interior_mass =
171 offline_data_->lumped_mass_matrix().local_element(index);
173 preliminary_map[index] = {index, interior_mass, position};
181 for (
const auto &[index, tuple] : preliminary_map) {
182 map.push_back(tuple);
185 return std::make_pair(name, map);
190 for (
const auto &[name, interior_map] : interior_maps_) {
192 const auto &options = get_options_from_name(interior_manifolds_, name);
193 if (options.find(
"instantaneous") == std::string::npos &&
194 options.find(
"time_averaged") == std::string::npos)
203 const auto received =
204 Utilities::MPI::gather(mpi_communicator_, interior_map);
206 if (Utilities::MPI::this_mpi_process(mpi_communicator_) == 0) {
208 std::ofstream output(base_name_ +
"-" + name +
"-R" +
209 Utilities::to_string(cycle, 4) +
"-points.dat");
211 output << std::scientific << std::setprecision(14);
213 output <<
"#\n# position\tinterior mass\n";
215 unsigned int rank = 0;
216 for (
const auto &entries : received) {
217 output <<
"# rank " << rank++ <<
"\n";
218 for (
const auto &entry : entries) {
219 const auto &[index, mass_i, x_i] = entry;
220 output << x_i <<
"\t" << mass_i <<
"\n";
224 output << std::flush;
237 boundary_maps_.clear();
239 boundary_manifolds_.begin(),
240 boundary_manifolds_.end(),
241 std::inserter(boundary_maps_, boundary_maps_.end()),
242 [
this, n_owned](
auto it) {
243 const auto &[name, expression, option] = it;
244 FunctionParser<dim> level_set_function(expression);
246 std::vector<boundary_point> map;
248 for (const auto &entry : offline_data_->boundary_map()) {
250 if (entry.first >= n_owned)
254 if (offline_data_->affine_constraints().is_constrained(
255 offline_data_->scalar_partitioner()->local_to_global(
259 const auto &[normal, normal_mass, boundary_mass, id, position] =
261 if (std::abs(level_set_function.value(position)) < 1.e-12)
262 map.push_back({entry.first,
269 return std::make_pair(name, map);
276 for (
const auto &[name, boundary_map] : boundary_maps_) {
278 const auto &options = get_options_from_name(boundary_manifolds_, name);
279 if (options.find(
"instantaneous") == std::string::npos &&
280 options.find(
"time_averaged") == std::string::npos)
289 const auto received =
290 Utilities::MPI::gather(mpi_communicator_, boundary_map);
292 if (Utilities::MPI::this_mpi_process(mpi_communicator_) == 0) {
294 std::ofstream output(base_name_ +
"-" + name +
"-R" +
295 Utilities::to_string(cycle, 4) +
"-points.dat");
297 output << std::scientific << std::setprecision(14);
299 output <<
"#\n# position\tnormal\tnormal mass\tboundary mass\n";
301 unsigned int rank = 0;
302 for (
const auto &entries : received) {
303 output <<
"# rank " << rank++ <<
"\n";
304 for (
const auto &entry : entries) {
305 const auto &[index, n_i, nm_i, bm_i, id, x_i] = entry;
306 output << x_i <<
"\t" << n_i <<
"\t" << nm_i <<
"\t" << bm_i
311 output << std::flush;
319 const auto &names = View::primitive_component_names;
320 header_ = std::accumulate(
324 [](
const std::string &description,
const std::string &name) {
325 return description.empty()
326 ? (std::string(
"primitive state (") + name)
327 : (description +
", " + name);
329 ")\t and 2nd moments\n";
333 template <
typename Description,
int dim,
typename Number>
334 void Quantities<Description, dim, Number>::clear_statistics()
336 const auto reset = [](
const auto &manifold_map,
auto &statistics_map) {
337 for (
const auto &[name, data_map] : manifold_map) {
338 const auto n_entries = data_map.size();
339 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
340 statistics_map[name];
341 val_old.resize(n_entries);
342 val_new.resize(n_entries);
343 val_sum.resize(n_entries);
344 t_old = t_new = t_sum = 0.;
350 interior_statistics_.clear();
351 reset(interior_maps_, interior_statistics_);
352 interior_time_series_.clear();
354 boundary_statistics_.clear();
355 reset(boundary_maps_, boundary_statistics_);
356 boundary_time_series_.clear();
360 template <
typename Description,
int dim,
typename Number>
361 template <
typename po
int_type,
typename value_type>
362 value_type Quantities<Description, dim, Number>::internal_accumulate(
364 const std::vector<point_type> &points_vector,
365 std::vector<value_type> &val_new)
367 const auto &U = std::get<0>(state_vector);
369 value_type spatial_average;
370 Number mass_sum = Number(0.);
373 points_vector.begin(),
376 [&](
auto point) -> value_type {
377 const auto i = std::get<0>(point);
382 constexpr auto index =
383 std::is_same<point_type, interior_point>::value ? 1 : 3;
384 const auto mass_i = std::get<index>(point);
386 const auto U_i = U.get_tensor(i);
387 const auto view = hyperbolic_system_->template view<dim, Number>();
388 const auto primitive_state = view.to_primitive_state(U_i);
391 std::get<0>(result) = primitive_state;
393 std::get<1>(result) = schur_product(primitive_state, primitive_state);
396 std::get<0>(spatial_average) += mass_i * std::get<0>(result);
397 std::get<1>(spatial_average) += mass_i * std::get<1>(result);
404 mass_sum = Utilities::MPI::sum(mass_sum, mpi_communicator_);
406 std::get<0>(spatial_average) =
407 Utilities::MPI::sum(std::get<0>(spatial_average), mpi_communicator_);
408 std::get<1>(spatial_average) =
409 Utilities::MPI::sum(std::get<1>(spatial_average), mpi_communicator_);
413 std::get<0>(spatial_average) /= mass_sum;
414 std::get<1>(spatial_average) /= mass_sum;
416 return spatial_average;
420 template <
typename Description,
int dim,
typename Number>
421 template <
typename value_type>
422 void Quantities<Description, dim, Number>::internal_write_out(
423 std::ostream &output,
424 const std::vector<value_type> &values,
433 const auto received = Utilities::MPI::gather(mpi_communicator_, values);
435 if (Utilities::MPI::this_mpi_process(mpi_communicator_) == 0) {
437 output <<
"# " << header_;
439 unsigned int rank = 0;
440 for (
const auto &entries : received) {
441 output <<
"# rank " << rank++ <<
"\n";
442 for (
const auto &entry : entries) {
443 const auto &[state, state_square] = entry;
444 output << scale * state <<
"\t" << scale * state_square <<
"\n";
448 output << std::flush;
453 template <
typename Description,
int dim,
typename Number>
454 template <
typename value_type>
455 void Quantities<Description, dim, Number>::internal_write_out_time_series(
456 std::ostream &output,
457 const std::vector<std::tuple<Number, value_type>> &values,
460 if (Utilities::MPI::this_mpi_process(mpi_communicator_) == 0) {
463 output <<
"# time t\t" << header_;
465 for (
const auto &entry : values) {
466 const auto t = std::get<0>(entry);
467 const auto &[state, state_square] = std::get<1>(entry);
469 output << t <<
"\t" << state <<
"\t" << state_square <<
"\n";
472 output << std::flush;
477 template <
typename Description,
int dim,
typename Number>
482 std::cout <<
"Quantities<dim, Number>::accumulate()" << std::endl;
485 const auto accumulate = [&](
const auto &point_maps,
486 const auto &manifolds,
489 for (
const auto &[name, point_map] : point_maps) {
492 const auto &options = get_options_from_name(manifolds, name);
495 if (options.find(
"time_averaged") == std::string::npos &&
496 options.find(
"space_averaged") == std::string::npos)
499 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
502 std::swap(t_old, t_new);
503 std::swap(val_old, val_new);
507 const auto spatial_average =
508 internal_accumulate(state_vector, point_map, val_new);
520 const Number tau = t_new - t_old;
522 for (std::size_t i = 0; i < val_sum.size(); ++i) {
523 std::get<0>(val_sum[i]) += 0.5 * tau * std::get<0>(val_old[i]);
524 std::get<0>(val_sum[i]) += 0.5 * tau * std::get<0>(val_new[i]);
525 std::get<1>(val_sum[i]) += 0.5 * tau * std::get<1>(val_old[i]);
526 std::get<1>(val_sum[i]) += 0.5 * tau * std::get<1>(val_new[i]);
532 time_series[name].push_back({t, spatial_average});
536 accumulate(interior_maps_,
538 interior_statistics_,
539 interior_time_series_);
541 accumulate(boundary_maps_,
543 boundary_statistics_,
544 boundary_time_series_);
548 template <
typename Description,
int dim,
typename Number>
550 const StateVector &state_vector,
const Number t,
unsigned int cycle)
553 std::cout <<
"Quantities<dim, Number>::write_out()" << std::endl;
556 const auto write_out = [&](
const auto &point_maps,
557 const auto &manifolds,
560 for (
const auto &[name, point_map] : point_maps) {
563 const auto &options = get_options_from_name(manifolds, name);
566 base_name_ +
"-" + name +
"-R" + Utilities::to_string(cycle, 4);
572 if (options.find(
"instantaneous") != std::string::npos) {
574 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
579 if (options.find(
"time_averaged") == std::string::npos &&
580 options.find(
"space_averaged") == std::string::npos)
581 internal_accumulate(state_vector, point_map, val_new);
583 AssertThrow(t_new == t, dealii::ExcInternalError());
585 std::ofstream output(prefix +
"-instantaneous.dat");
587 output << std::scientific << std::setprecision(14);
588 output <<
"# at t = " << t << std::endl;
590 internal_write_out(output, val_new, Number(1.));
597 if (options.find(
"time_averaged") != std::string::npos) {
599 std::ofstream output(prefix +
"-time_averaged.dat");
601 auto &[val_old, val_new, val_sum, t_old, t_new, t_sum] =
604 output << std::scientific << std::setprecision(14);
605 output <<
"# averaged from t = " << t_new - t_sum
606 <<
" to t = " << t_new << std::endl;
608 internal_write_out(output, val_sum, Number(1.) / t_sum);
615 if (options.find(
"space_averaged") != std::string::npos) {
617 auto &series = time_series[name];
619 std::ofstream output;
620 const auto file_name = base_name_ +
"-" + name +
"-R" +
621 Utilities::to_string(time_series_cycle_, 4) +
622 "-space_averaged_time_series.dat";
624 output << std::scientific << std::setprecision(14);
627 first_cycle_ =
false;
628 output.open(file_name, std::ofstream::out | std::ofstream::trunc);
629 internal_write_out_time_series(output, series,
false);
633 output.open(file_name, std::ofstream::out | std::ofstream::app);
634 internal_write_out_time_series(output, series,
true);
642 write_out(interior_maps_,
644 interior_statistics_,
645 interior_time_series_);
647 write_out(boundary_maps_,
649 boundary_statistics_,
650 boundary_time_series_);
652 if (clear_temporal_statistics_on_writeout_)
typename Description::HyperbolicSystem HyperbolicSystem
typename View::StateVector StateVector
Quantities(const MPI_Comm &mpi_communicator, const HyperbolicSystem &hyperbolic_system, const OfflineData< dim, Number > &offline_data, const std::string &subsection="/Quantities")
void prepare(const std::string &name, unsigned int cycle)
#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)