From e15c232217a3c3aa124ac204b7849cdc712c5e16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olivier=20Gu=C3=A9vremont?= Date: Mon, 5 Aug 2024 21:58:13 -0400 Subject: [PATCH] Tracer flow rate at boundaries (#1197) Description An application of the tracer physics is to obtain the residence time distribution (RTD) of a flow through a domain. This PR adds the post-processing option that outputs the tracer flow rate at every boundary, providing the data for RTD analyses. In addition, a bug where the conditions to allow boundary conditions to update their values were too strict was rectified. Each physics is now responsible for deciding when to update their boundary conditions, and this condition was already based on whether boundary conditions were time-dependent. Co-authored-by: Amishga Alphonius <107414376+AmishgaAlphonius@users.noreply.github.com> Co-authored-by: hepap <47506601+hepap@users.noreply.github.com> Former-commit-id: fd4481eb5847ec29409f6750165dc97fa2682458 --- CHANGELOG.md | 8 + .../lethe-fluid/tracer_flow_rate.output | 35 +++ .../lethe-fluid/tracer_flow_rate.prm | 88 ++++++++ doc/source/parameters/cfd/post_processing.rst | 4 +- include/core/parameters.h | 6 + include/solvers/tracer.h | 35 ++- source/core/parameters.cc | 33 ++- source/fem-dem/cfd_dem_coupling.cc | 23 +- source/fem-dem/fluid_dynamics_vans.cc | 23 +- source/fem-dem/gls_sharp_navier_stokes.cc | 22 +- source/solvers/fluid_dynamics_matrix_free.cc | 19 +- source/solvers/gls_navier_stokes.cc | 22 +- source/solvers/navier_stokes_base.cc | 2 - source/solvers/tracer.cc | 209 ++++++++++++++++++ 14 files changed, 423 insertions(+), 106 deletions(-) create mode 100755 applications_tests/lethe-fluid/tracer_flow_rate.output create mode 100755 applications_tests/lethe-fluid/tracer_flow_rate.prm diff --git a/CHANGELOG.md b/CHANGELOG.md index ea318748c3..d2a3bb6f20 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,14 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/). - MINOR Renamed all the files related to the `lethe-fluid-matrix-free` application: `mf_navier_stokes` to `fluid_dynamics_matrix_free` and `mf_navier_stokes_operators` to `fluid_dynamics_matrix_free_operators`. The main class `MFNavierStokesSolver` was also renamed to `FluidDynamicsMatrixFree`. [#1222](https://github.com/chaos-polymtl/lethe/pull/1222) +### Added + +- MINOR A post-processing feature for the tracer flow rate through boundaries is added. This allows to produce data useful for residence time distribution analyses. [#1197](https://github.com/chaos-polymtl/lethe/pull/1197) + +### Fixed + +- MINOR Conditions for whether a boundary condition should be updated were made less restrictive. Now, each physics is responsible for updating or not their boundary conditions. [#1197](https://github.com/chaos-polymtl/lethe/pull/1197) + ### Removed - MINOR Removed the gear3 DEM integrator. It was never used, it did not possess any unit test and it was not even clear if it still worked. [#1221](https://github.com/chaos-polymtl/lethe/pull/1211) diff --git a/applications_tests/lethe-fluid/tracer_flow_rate.output b/applications_tests/lethe-fluid/tracer_flow_rate.output new file mode 100755 index 0000000000..9753acdea5 --- /dev/null +++ b/applications_tests/lethe-fluid/tracer_flow_rate.output @@ -0,0 +1,35 @@ +Running on 1 MPI rank(s)... + Number of active cells: 64 + Number of degrees of freedom: 243 + Volume of triangulation: 0.04 + Number of tracer degrees of freedom: 81 +Tracer statistics : + Min : 0 + Max : 0.788675 + Average : 0.0625 + Std-Dev : 0.19432 +------------------ +Tracer flow rates +------------------ +Tracer flow rate at the boundaries: + boundary 0: 8.2 + boundary 1: 0 + boundary 2: 0 + boundary 3: 0 + +***************************** +Steady iteration: 1/1 +***************************** +Tracer statistics : + Min : 1 + Max : 1 + Average : 1 + Std-Dev : 3.37941e-09 +------------------ +Tracer flow rates +------------------ +Tracer flow rate at the boundaries: + boundary 0: 0.2 + boundary 1: -0.200001 + boundary 2: 9.02388e-10 + boundary 3: 2.18592e-09 diff --git a/applications_tests/lethe-fluid/tracer_flow_rate.prm b/applications_tests/lethe-fluid/tracer_flow_rate.prm new file mode 100755 index 0000000000..29b75fb0ba --- /dev/null +++ b/applications_tests/lethe-fluid/tracer_flow_rate.prm @@ -0,0 +1,88 @@ +set dimension = 2 + +subsection simulation control + set output frequency = 1 +end + +subsection physical properties + subsection fluid 0 + set kinematic viscosity = 0.01 + set tracer diffusivity = 1 + end +end + +subsection multiphysics + set tracer = true +end + +subsection mesh + set type = dealii + set grid type = subdivided_hyper_rectangle + set grid arguments = 1, 1: -0.1, -0.1: 0.1, 0.1 : true + set initial refinement = 3 +end + +subsection boundary conditions + set number = 4 + set time dependent = true + subsection bc 0 + set id = 0 + set type = function + subsection u + set Function expression = 1 + end + end + subsection bc 1 + set id = 1 + set type = outlet + end + subsection bc 2 + set id = 2 + set type = noslip + end + subsection bc 3 + set id = 3 + set type = noslip + end +end + +subsection boundary conditions tracer + set number = 1 + subsection bc 0 + set id = 0 + set type = dirichlet + subsection dirichlet + set Function expression = 1 + end + end +end + +subsection post-processing + set verbosity = verbose + set calculate tracer statistics = true + set tracer statistics name = tracer_statistics + set calculate tracer flow rate = true + set tracer flow rate name = tracer_flow_rate_test +end + +subsection non-linear solver + subsection fluid dynamics + set verbosity = quiet + set tolerance = 1e-4 + set max iterations = 10 + end + subsection tracer + set verbosity = quiet + set tolerance = 1e-7 + set max iterations = 30 + end +end + +subsection linear solver + subsection fluid dynamics + set verbosity = quiet + end + subsection tracer + set verbosity = quiet + end +end diff --git a/doc/source/parameters/cfd/post_processing.rst b/doc/source/parameters/cfd/post_processing.rst index b829f1da24..34e69cd684 100644 --- a/doc/source/parameters/cfd/post_processing.rst +++ b/doc/source/parameters/cfd/post_processing.rst @@ -55,9 +55,11 @@ This subsection controls the post-processing other than the forces and torque on #--------------------------------------------------- # Multiphysics post-processing #--------------------------------------------------- - # Tracer statistics + # Tracer postprocessing set calculate tracer statistics = false set tracer statistics name = tracer_statistics + set calculate tracer flow rate = false + set tracer flow rate name = tracer_flow_rate # Thermal postprocessing set postprocessed fluid = both diff --git a/include/core/parameters.h b/include/core/parameters.h index d31ffd8ef3..a9423789c4 100644 --- a/include/core/parameters.h +++ b/include/core/parameters.h @@ -883,6 +883,9 @@ namespace Parameters /// Enable flow rate post-processing bool calculate_flow_rate; + /// Enable tracer flow rate post-processing + bool calculate_tracer_flow_rate; + /// Set initial time to start calculations for velocities double initial_time; @@ -901,6 +904,9 @@ namespace Parameters /// Prefix for flow rate output std::string flow_rate_output_name; + /// Prefix for tracer flow rate output + std::string tracer_flow_rate_output_name; + /// Prefix for the enstrophy output std::string enstrophy_output_name; diff --git a/include/solvers/tracer.h b/include/solvers/tracer.h index 75795644c9..f2a95e0008 100644 --- a/include/solvers/tracer.h +++ b/include/solvers/tracer.h @@ -77,6 +77,8 @@ class Tracer : public AuxiliaryPhysics simulation_parameters.fem_parameters.tracer_order); mapping = std::make_shared>(*fe); cell_quadrature = std::make_shared>(fe->degree + 1); + face_quadrature = + std::make_shared>(fe->degree + 1); } else { @@ -85,6 +87,7 @@ class Tracer : public AuxiliaryPhysics simulation_parameters.fem_parameters.tracer_order); mapping = std::make_shared>(fe->degree); cell_quadrature = std::make_shared>(fe->degree + 1); + face_quadrature = std::make_shared>(fe->degree + 1); } // Allocate solution transfer @@ -366,6 +369,27 @@ class Tracer : public AuxiliaryPhysics void calculate_tracer_statistics(); + + /** + * Post-processing. Calculate the tracer flow rate at every boundary. + * + * @param current_solution_fd current solution for the fluid dynamics, parsed + * by postprocess + * + * @return values of the tracer flow rate at each boundary + * + * @tparam GlobalVectorType The type of the global vector used for the tracer physic + */ + template + std::vector + postprocess_tracer_flow_rate(const GlobalVectorType ¤t_solution_fd); + + /** + * @brief Writes the tracer flow rates to an output file + */ + void + write_tracer_flow_rates(const std::vector tracer_flow_rate_vector); + /** * @brief Writes the tracer statistics to an output file */ @@ -384,11 +408,13 @@ class Tracer : public AuxiliaryPhysics std::shared_ptr simulation_control; DoFHandler dof_handler; - // Finite element spce + // Finite element space std::shared_ptr> fe; + // Mapping and Quadrature - std::shared_ptr> mapping; - std::shared_ptr> cell_quadrature; + std::shared_ptr> mapping; + std::shared_ptr> cell_quadrature; + std::shared_ptr> face_quadrature; ConvergenceTable error_table; @@ -422,8 +448,9 @@ class Tracer : public AuxiliaryPhysics // Assemblers for the matrix and rhs std::vector>> assemblers; - // Tracer statistics table + // Tracer post-processing tables TableHandler statistics_table; + TableHandler tracer_flow_rate_table; }; diff --git a/source/core/parameters.cc b/source/core/parameters.cc index 83b31d2914..8368f6fbc2 100644 --- a/source/core/parameters.cc +++ b/source/core/parameters.cc @@ -1806,6 +1806,12 @@ namespace Parameters Patterns::Bool(), "Enable calculation of flow rate at boundaries."); + prm.declare_entry( + "calculate tracer flow rate", + "false", + Patterns::Bool(), + "Enable calculation of tracer flow rate at boundaries."); + prm.declare_entry( "initial time", "0.0", @@ -1827,6 +1833,11 @@ namespace Parameters Patterns::FileName(), "File output volumetric flux"); + prm.declare_entry("tracer flow rate name", + "tracer_flow_rate", + Patterns::FileName(), + "Output file name for tracer flow rate"); + prm.declare_entry("enstrophy name", "enstrophy", Patterns::FileName(), @@ -1997,16 +2008,18 @@ namespace Parameters prm.get_bool("calculate apparent viscosity"); calculate_average_velocities = prm.get_bool("calculate average velocities"); - calculate_pressure_drop = prm.get_bool("calculate pressure drop"); - inlet_boundary_id = prm.get_integer("inlet boundary id"); - outlet_boundary_id = prm.get_integer("outlet boundary id"); - calculate_flow_rate = prm.get_bool("calculate flow rate"); - initial_time = prm.get_double("initial time"); - kinetic_energy_output_name = prm.get("kinetic energy name"); - pressure_drop_output_name = prm.get("pressure drop name"); - flow_rate_output_name = prm.get("flow rate name"); - enstrophy_output_name = prm.get("enstrophy name"); - pressure_power_output_name = prm.get("pressure power name"); + calculate_pressure_drop = prm.get_bool("calculate pressure drop"); + inlet_boundary_id = prm.get_integer("inlet boundary id"); + outlet_boundary_id = prm.get_integer("outlet boundary id"); + calculate_flow_rate = prm.get_bool("calculate flow rate"); + calculate_tracer_flow_rate = prm.get_bool("calculate tracer flow rate"); + initial_time = prm.get_double("initial time"); + kinetic_energy_output_name = prm.get("kinetic energy name"); + pressure_drop_output_name = prm.get("pressure drop name"); + flow_rate_output_name = prm.get("flow rate name"); + tracer_flow_rate_output_name = prm.get("tracer flow rate name"); + enstrophy_output_name = prm.get("enstrophy name"); + pressure_power_output_name = prm.get("pressure power name"); viscous_dissipation_output_name = prm.get("viscous dissipation name"); apparent_viscosity_output_name = prm.get("apparent viscosity name"); output_frequency = prm.get_integer("output frequency"); diff --git a/source/fem-dem/cfd_dem_coupling.cc b/source/fem-dem/cfd_dem_coupling.cc index c7c96c5d58..400defedd3 100644 --- a/source/fem-dem/cfd_dem_coupling.cc +++ b/source/fem-dem/cfd_dem_coupling.cc @@ -1373,24 +1373,11 @@ CFDDEMSolver::solve() while (this->simulation_control->integrate()) { this->simulation_control->print_progression(this->pcout); - bool refinement_step; - if (this->simulation_parameters.mesh_adaptation.refinement_at_frequency) - refinement_step = - this->simulation_control->get_step_number() % - this->simulation_parameters.mesh_adaptation.frequency != - 0; - else - refinement_step = this->simulation_control->get_step_number() == 0; - if (refinement_step || - this->simulation_parameters.mesh_adaptation.type == - Parameters::MeshAdaptation::Type::none || - this->simulation_control->is_at_start()) - { - // We allow the physics to update their boundary conditions - // according to their own parameters - this->update_boundary_conditions(); - this->multiphysics->update_boundary_conditions(); - } + + // We allow the physics to update their boundary conditions + // according to their own parameters + this->update_boundary_conditions(); + this->multiphysics->update_boundary_conditions(); this->dynamic_flow_control(); diff --git a/source/fem-dem/fluid_dynamics_vans.cc b/source/fem-dem/fluid_dynamics_vans.cc index 5f14dd6ca0..cf36f158d8 100644 --- a/source/fem-dem/fluid_dynamics_vans.cc +++ b/source/fem-dem/fluid_dynamics_vans.cc @@ -1994,24 +1994,11 @@ FluidDynamicsVANS::solve() while (this->simulation_control->integrate()) { this->simulation_control->print_progression(this->pcout); - bool refinement_step; - if (this->simulation_parameters.mesh_adaptation.refinement_at_frequency) - refinement_step = - this->simulation_control->get_step_number() % - this->simulation_parameters.mesh_adaptation.frequency != - 0; - else - refinement_step = this->simulation_control->get_step_number() == 0; - if ((refinement_step || - this->simulation_parameters.mesh_adaptation.type == - Parameters::MeshAdaptation::Type::none || - this->simulation_control->is_at_start())) - { - // We allow the physics to update their boundary conditions - // according to their own parameters - this->update_boundary_conditions(); - this->multiphysics->update_boundary_conditions(); - } + + // We allow the physics to update their boundary conditions + // according to their own parameters + this->update_boundary_conditions(); + this->multiphysics->update_boundary_conditions(); this->dynamic_flow_control(); diff --git a/source/fem-dem/gls_sharp_navier_stokes.cc b/source/fem-dem/gls_sharp_navier_stokes.cc index b69bb05db1..6875b33928 100644 --- a/source/fem-dem/gls_sharp_navier_stokes.cc +++ b/source/fem-dem/gls_sharp_navier_stokes.cc @@ -4885,24 +4885,10 @@ GLSSharpNavierStokesSolver::solve() this->forcing_function->set_time( this->simulation_control->get_current_time()); - bool refinement_step; - if (this->simulation_parameters.mesh_adaptation.refinement_at_frequency) - refinement_step = - this->simulation_control->get_step_number() % - this->simulation_parameters.mesh_adaptation.frequency != - 0; - else - refinement_step = this->simulation_control->get_step_number() == 0; - if (refinement_step || - this->simulation_parameters.mesh_adaptation.type == - Parameters::MeshAdaptation::Type::none || - this->simulation_control->is_at_start()) - { - // We allow the physics to update their boundary conditions - // according to their own parameters - this->update_boundary_conditions(); - this->multiphysics->update_boundary_conditions(); - } + // We allow the physics to update their boundary conditions + // according to their own parameters + this->update_boundary_conditions(); + this->multiphysics->update_boundary_conditions(); if (some_particles_are_coupled == false) integrate_particles(); diff --git a/source/solvers/fluid_dynamics_matrix_free.cc b/source/solvers/fluid_dynamics_matrix_free.cc index 3e25e5daaa..2ebccf696d 100644 --- a/source/solvers/fluid_dynamics_matrix_free.cc +++ b/source/solvers/fluid_dynamics_matrix_free.cc @@ -1915,23 +1915,8 @@ FluidDynamicsMatrixFree::solve() this->forcing_function->set_time( this->simulation_control->get_current_time()); - bool refinement_step; - if (this->simulation_parameters.mesh_adaptation.refinement_at_frequency) - refinement_step = - this->simulation_control->get_step_number() % - this->simulation_parameters.mesh_adaptation.frequency != - 0; - else - refinement_step = this->simulation_control->get_step_number() == 0; - if ((refinement_step || - this->simulation_parameters.mesh_adaptation.type == - Parameters::MeshAdaptation::Type::none || - this->simulation_control->is_at_start()) && - this->simulation_parameters.boundary_conditions.time_dependent) - { - update_boundary_conditions(); - this->multiphysics->update_boundary_conditions(); - } + update_boundary_conditions(); + this->multiphysics->update_boundary_conditions(); this->simulation_control->print_progression(this->pcout); this->dynamic_flow_control(); diff --git a/source/solvers/gls_navier_stokes.cc b/source/solvers/gls_navier_stokes.cc index e0509d8ef7..87b8843f81 100644 --- a/source/solvers/gls_navier_stokes.cc +++ b/source/solvers/gls_navier_stokes.cc @@ -1826,24 +1826,10 @@ GLSNavierStokesSolver::solve() this->forcing_function->set_time( this->simulation_control->get_current_time()); - bool refinement_step; - if (this->simulation_parameters.mesh_adaptation.refinement_at_frequency) - refinement_step = - this->simulation_control->get_step_number() % - this->simulation_parameters.mesh_adaptation.frequency != - 0; - else - refinement_step = this->simulation_control->get_step_number() == 0; - if (refinement_step || - this->simulation_parameters.mesh_adaptation.type == - Parameters::MeshAdaptation::Type::none || - this->simulation_control->is_at_start()) - { - // We allow the physics to update their boundary conditions - // according to their own parameters - this->update_boundary_conditions(); - this->multiphysics->update_boundary_conditions(); - } + // We allow the physics to update their boundary conditions + // according to their own parameters + this->update_boundary_conditions(); + this->multiphysics->update_boundary_conditions(); this->simulation_control->print_progression(this->pcout); this->dynamic_flow_control(); diff --git a/source/solvers/navier_stokes_base.cc b/source/solvers/navier_stokes_base.cc index fa85fb5037..860dffcd6f 100644 --- a/source/solvers/navier_stokes_base.cc +++ b/source/solvers/navier_stokes_base.cc @@ -1548,8 +1548,6 @@ NavierStokesBase::postprocess_fd(bool firstIter) } } - // this->computing_timer.leave_subsection("Postprocessing"); - if (!firstIter) { // Calculate forces on the boundary conditions diff --git a/source/solvers/tracer.cc b/source/solvers/tracer.cc index 8d437822d9..776ff9cb22 100644 --- a/source/solvers/tracer.cc +++ b/source/solvers/tracer.cc @@ -429,6 +429,33 @@ Tracer::postprocess(bool first_iteration) 0) this->write_tracer_statistics(); } + + // Calculate tracer flow rate at every boundary + if (this->simulation_parameters.post_processing.calculate_tracer_flow_rate) + { + TimerOutput::Scope t(this->computing_timer, "Calculate tracer flow rate"); + + if (this->simulation_parameters.post_processing.verbosity == + Parameters::Verbosity::verbose) + { + announce_string(this->pcout, "Tracer flow rates"); + } + + std::vector tracer_flow_rates( + this->simulation_parameters.boundary_conditions.size); + if (multiphysics->fluid_dynamics_is_block()) + { + tracer_flow_rates = postprocess_tracer_flow_rate( + *multiphysics->get_block_solution(PhysicsID::fluid_dynamics)); + } + else + { + tracer_flow_rates = postprocess_tracer_flow_rate( + *multiphysics->get_solution(PhysicsID::fluid_dynamics)); + } + this->write_tracer_flow_rates(tracer_flow_rates); + } + if (this->simulation_parameters.timer.type == Parameters::Timer::Type::iteration) { @@ -530,6 +557,188 @@ Tracer::calculate_tracer_statistics() statistics_table.set_precision("std-dev", 12); } + +template +template +std::vector +Tracer::postprocess_tracer_flow_rate(const VectorType ¤t_solution_fd) +{ + const unsigned int n_q_points_face = this->face_quadrature->size(); + const MPI_Comm mpi_communicator = this->dof_handler.get_communicator(); + + // Initialize tracer information + std::vector tracer_values(n_q_points_face); + std::vector> tracer_gradient(n_q_points_face); + FEFaceValues fe_face_values_tracer(*this->mapping, + this->dof_handler.get_fe(), + *this->face_quadrature, + update_values | update_gradients | + update_JxW_values | + update_normal_vectors); + + // Initialize fluid dynamics information + const DoFHandler *dof_handler_fd = + multiphysics->get_dof_handler(PhysicsID::fluid_dynamics); + const FEValuesExtractors::Vector velocities(0); + std::vector> velocity_values(n_q_points_face); + FEFaceValues fe_face_values_fd(*this->mapping, + dof_handler_fd->get_fe(), + *this->face_quadrature, + update_values); + + // Initialize fluid properties + auto &properties_manager = + this->simulation_parameters.physical_properties_manager; + std::map> fields; + + const auto diffusivity_model = properties_manager.get_tracer_diffusivity(); + + std::vector tracer_diffusivity(n_q_points_face); + std::vector> face_quadrature_points; + + std::vector tracer_flow_rate_vector( + this->simulation_parameters.boundary_conditions.size, 0); + + // Get vector of all boundary conditions + const auto boundary_conditions_ids = + this->simulation_parameters.boundary_conditions.id; + + for (const auto &cell : this->dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned() && cell->at_boundary()) + { + for (const auto face : cell->face_indices()) + { + if (cell->face(face)->at_boundary()) + { + const auto boundary_id = + std::find(begin(boundary_conditions_ids), + end(boundary_conditions_ids), + cell->face(face)->boundary_id()); + + + unsigned int vector_index = + boundary_id - boundary_conditions_ids.begin(); + + // Gather tracer information + fe_face_values_tracer.reinit(cell, face); + fe_face_values_tracer.get_function_values( + this->present_solution, tracer_values); + fe_face_values_tracer.get_function_gradients( + this->present_solution, tracer_gradient); + + // We update the fields required by the diffusivity + // model + fields.clear(); + if (diffusivity_model->depends_on(field::levelset)) + { + std::vector levelset_values(n_q_points_face); + fields.insert( + std::pair>(field::levelset, + n_q_points_face)); + face_quadrature_points = + fe_face_values_tracer.get_quadrature_points(); + immersed_solid_signed_distance_function = + this->multiphysics + ->get_immersed_solid_signed_distance_function(); + this->immersed_solid_signed_distance_function->value_list( + face_quadrature_points, levelset_values); + set_field_vector(field::levelset, + levelset_values, + fields); + } + + diffusivity_model->vector_value(fields, tracer_diffusivity); + + // Get fluid dynamics active cell iterator + typename DoFHandler::active_cell_iterator cell_fd( + &(*(this->triangulation)), + cell->level(), + cell->index(), + dof_handler_fd); + + // Gather fluid dynamics information + fe_face_values_fd.reinit(cell_fd, face); + fe_face_values_fd[velocities].get_function_values( + current_solution_fd, velocity_values); + + // Loop on the quadrature points + for (unsigned int q = 0; q < n_q_points_face; q++) + { + Tensor<1, dim> normal_vector_tracer = + -fe_face_values_tracer.normal_vector(q); + + tracer_flow_rate_vector[vector_index] += + (-tracer_diffusivity[q] * tracer_gradient[q] * + normal_vector_tracer + + tracer_values[q] * velocity_values[q] * + normal_vector_tracer) * + fe_face_values_tracer.JxW(q); + } // end loop on quadrature points + } // end face is a boundary face + } // end loop on faces + } // end condition cell at boundary + } // end loop on cells + + + // Sum across all cores + for (unsigned int i_bc = 0; + i_bc < this->simulation_parameters.boundary_conditions.size; + ++i_bc) + tracer_flow_rate_vector[i_bc] = + Utilities::MPI::sum(tracer_flow_rate_vector[i_bc], mpi_communicator); + + return tracer_flow_rate_vector; +} + +template +void +Tracer::write_tracer_flow_rates( + const std::vector tracer_flow_rate_vector) +{ + // Fill table + this->tracer_flow_rate_table.add_value( + "time", this->simulation_control->get_current_time()); + this->tracer_flow_rate_table.set_precision("time", 12); + for (unsigned int i_bc = 0; + i_bc < this->simulation_parameters.boundary_conditions.size; + ++i_bc) + { + this->tracer_flow_rate_table.add_value("flow-rate-" + + Utilities::int_to_string(i_bc, + 2), + tracer_flow_rate_vector[i_bc]); + this->tracer_flow_rate_table.set_scientific( + "flow-rate-" + Utilities::int_to_string(i_bc, 2), true); + } + + // Console output + if (simulation_parameters.post_processing.verbosity == + Parameters::Verbosity::verbose) + { + this->pcout << "Tracer flow rate at the boundaries: " << std::endl; + for (unsigned int i_bc = 0; + i_bc < this->simulation_parameters.boundary_conditions.size; + ++i_bc) + this->pcout << "\t boundary " << i_bc << ": " + << tracer_flow_rate_vector[i_bc] << std::endl; + } + + auto mpi_communicator = triangulation->get_communicator(); + if ((simulation_control->get_step_number() % + this->simulation_parameters.post_processing.output_frequency == + 0) && + Utilities::MPI::this_mpi_process(mpi_communicator) == 0) + { + std::string filename = + simulation_parameters.simulation_control.output_folder + + simulation_parameters.post_processing.tracer_flow_rate_output_name + + ".dat"; + std::ofstream output(filename.c_str()); + this->tracer_flow_rate_table.write_text(output); + } +} + template void Tracer::write_tracer_statistics()