Tracer flow rate at boundaries (chaos-polymtl#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 <[email protected]> Co-authored-by: hepap <[email protected]> Former-commit-id: fd4481e
M-Badri · Aug 6, 2024 · e15c232 · e15c232
1 parent c87b389
commit e15c232
Show file tree

Hide file tree

Showing 14 changed files with 423 additions and 106 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -77,6 +77,8 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
           simulation_parameters.fem_parameters.tracer_order);
         mapping         = std::make_shared<MappingFE<dim>>(*fe);
         cell_quadrature = std::make_shared<QGaussSimplex<dim>>(fe->degree + 1);
+        face_quadrature =
+          std::make_shared<QGaussSimplex<dim - 1>>(fe->degree + 1);
       }
     else
       {
@@ -85,6 +87,7 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
           simulation_parameters.fem_parameters.tracer_order);
         mapping         = std::make_shared<MappingQ<dim>>(fe->degree);
         cell_quadrature = std::make_shared<QGauss<dim>>(fe->degree + 1);
+        face_quadrature = std::make_shared<QGauss<dim - 1>>(fe->degree + 1);
       }
 
     // Allocate solution transfer
@@ -366,6 +369,27 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
   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 <typename GlobalVectorType>
+  std::vector<double>
+  postprocess_tracer_flow_rate(const GlobalVectorType &current_solution_fd);
+
+  /**
+   * @brief Writes the tracer flow rates to an output file
+   */
+  void
+  write_tracer_flow_rates(const std::vector<double> tracer_flow_rate_vector);
+
   /**
    * @brief Writes the tracer statistics to an output file
    */
@@ -384,11 +408,13 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
   std::shared_ptr<SimulationControl> simulation_control;
   DoFHandler<dim>                    dof_handler;
 
-  // Finite element spce
+  // Finite element space
   std::shared_ptr<FiniteElement<dim>> fe;
+
   // Mapping and Quadrature
-  std::shared_ptr<Mapping<dim>>    mapping;
-  std::shared_ptr<Quadrature<dim>> cell_quadrature;
+  std::shared_ptr<Mapping<dim>>        mapping;
+  std::shared_ptr<Quadrature<dim>>     cell_quadrature;
+  std::shared_ptr<Quadrature<dim - 1>> face_quadrature;
 
 
   ConvergenceTable error_table;
@@ -422,8 +448,9 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
   // Assemblers for the matrix and rhs
   std::vector<std::shared_ptr<TracerAssemblerBase<dim>>> 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
@@ -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
@@ -1373,24 +1373,11 @@ CFDDEMSolver<dim>::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
@@ -1994,24 +1994,11 @@ FluidDynamicsVANS<dim>::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
@@ -4885,24 +4885,10 @@ GLSSharpNavierStokesSolver<dim>::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();