Fix heat flux postprocessor material_id processing (#988)

Description of the problem
In the previous implementation of the heat flux postprocessor, there was a confusion related to the material_id of a cell. As a cell's material_id is defined in the mesh, fluids moving from a cell to another did not have the right material_id when outputted. Futhermore, some dealii meshs have an intrinc way of numbering material_ids (e.g. colorized subdivided_hyper_rectangle) that were leading to inadequate comparison of the material_ids.
Description of the solution
As of this PR, in the heat flux postprocessor, the material_id is only used to indentify cell material when solids are present. When simulating with fluids and solids, cells of the mesh in which the fluids lie should have a material_id of 0.
For fluids-only simulations, at the moment, their heat fluxes will be output on the entire domain, as the current DataPostprocessor can only hold 1 dof_handler. During postprocessing, the user could clip the domain of interest of each fluid using the phase or filtered_phase fields outputted.
How Has This Been Tested?
Visual inspection of outputted results for a:

 single fluid simulation,

 1 fluid + 1 solid simulation, and

 2 fluids simulation.

Future changes
For multiple-fluids simulations, a DataPostprocessor handling multiple dof_handlers to allow postprocessing of quantities of interest of fluids in their respective subdomains will have to be added to the code.

include/solvers/postprocessors.h

@@ -511,73 +511,115 @@ class DensityPostprocessor : public DataPostprocessorScalar<dim>
 
 
 /**
- * @class Calculates heat flux (q = -k * grad T) of a material at every
- * quadrature point.
+ * @class Postprocessor that calculates local heat fluxes \f$(q = -k * \nabla
+ * T)\f$ of a material at every evaluation point of the domain.
+ *
  * @param p_thermal_conductivity_model Thermal conductivity model of the
  * material
  * @param p_material_string Type of material: "f" (fluid) or "s" (solid)
  * @param p_id Either fluid ID or solid ID (as initialized in the physical
  * properties) used to name the postprocessed variable vector.
- * @param p_material_id ID corresponding to the material. Material IDs the
- * cumulative IDs of fluids and solids by convention, fluids are counted before
- * solids. This is used to identify which cell lies in which material.
+ * @param p_mesh_material_id ID corresponding to the material as identified in
+ * the mesh. By convention, fluids are given the material_id 0 and solids follow
+ * with ids 1, 2, etc. This is used to identify which cell lies in which
+ * material.
+ * @param p_solid_phase_present Boolean indicating if the simulation has solids.
+ * If there is no solid, the @p material_id of cells is not checked.
  */
 template <int dim>
 class HeatFluxPostprocessor : public DataPostprocessorVector<dim>
 {
 public:
   HeatFluxPostprocessor(const std::shared_ptr<ThermalConductivityModel>
                                            &p_thermal_conductivity_model,
-                        const std::string  &p_material_string = "f",
-                        const unsigned int &p_id              = 0,
-                        const unsigned int &p_material_id     = 0)
+                        const std::string  &p_material_string     = "f",
+                        const unsigned int &p_id                  = 0,
+                        const unsigned int &p_mesh_material_id    = 0,
+                        const bool         &p_solid_phase_present = false)
     : DataPostprocessorVector<dim>("heat_flux_" + p_material_string +
                                      Utilities::to_string(p_id, 2),
                                    update_values | update_gradients)
     , thermal_conductivity_model(p_thermal_conductivity_model)
-    , material_id(p_material_id)
+    , mesh_material_id(p_mesh_material_id)
+    , solid_phase_present(p_solid_phase_present)
   {}
 
+  /**
+   * @brief Computes local heat flux quantities from temperature field and
+   * material properties.
+   *
+   * @param inputs Structure that contains all temperature information needed to
+   * evaluate local heat fluxes.
+   * @param computed_quantities Vector field containing local heat flux vectors
+   * evaluated at evaluation points.
+   */
   virtual void
   evaluate_scalar_field(
     const DataPostprocessorInputs::Scalar<dim> &inputs,
     std::vector<Vector<double>> &computed_quantities) const override
   {
-    const unsigned int n_quadrature_points = computed_quantities.size();
-    AssertDimension(inputs.solution_gradients.size(), n_quadrature_points);
-
-    // Get Cell
-    const auto cell = inputs.template get_cell<dim>();
+    const unsigned int n_evaluation_points = computed_quantities.size();
+    AssertDimension(inputs.solution_gradients.size(), n_evaluation_points);
 
-    // Check if the cell lies in the material of interest
-    if (cell->material_id() == material_id)
+    //  If solid materials are present, identify each material's subdomain.
+    if (solid_phase_present)
       {
-        std::map<field, double> field_values;
-        for (unsigned int q = 0; q < n_quadrature_points; ++q)
-          {
-            AssertDimension(computed_quantities[q].size(), dim);
-            field_values[field::temperature] = inputs.solution_values[q];
+        // Get Cell
+        const auto cell = inputs.template get_cell<dim>();
 
-            for (unsigned int i = 0; i < dim; ++i)
-              {
-                computed_quantities[q][i] =
-                  -thermal_conductivity_model->value(field_values) *
-                  inputs.solution_gradients[q][i];
-              }
+        // Check if the cell lies in the material of interest
+        if (cell->material_id() == mesh_material_id)
+          {
+            compute_quantities(inputs, computed_quantities);
+          }
+        // Apply null values to irrelevant subdomain
+        else
+          {
+            Vector<double>              null_vector(dim);
+            std::vector<Vector<double>> null_computed_quantities(
+              n_evaluation_points, null_vector);
+            computed_quantities = null_computed_quantities;
           }
       }
-    // Apply null values to irrelevant subdomain
+    // Only fluids are present
     else
       {
-        Vector<double>              null_vector(dim);
-        std::vector<Vector<double>> null_computed_quantities(
-          n_quadrature_points, null_vector);
-        computed_quantities = null_computed_quantities;
+        compute_quantities(inputs, computed_quantities);
       }
   }
 
 private:
+  /**
+   * @brief Computes local heat flux quantities. This function is used to
+   * improve readability of the HeatFluxPostprocessor::evaluate_scalar_field
+   * function and avoid repeated lines.
+   *
+   * @param inputs Structure that contains all temperature information needed to
+   * evaluate local heat fluxes.
+   * @param computed_quantities Vector field containing local heat flux vectors
+   * evaluated at evaluation points.
+   */
+  void
+  compute_quantities(const DataPostprocessorInputs::Scalar<dim> &inputs,
+                     std::vector<Vector<double>> &computed_quantities) const
+  {
+    std::map<field, double> field_values;
+    for (unsigned int p = 0; p < computed_quantities.size(); ++p)
+      {
+        AssertDimension(computed_quantities[p].size(), dim);
+        field_values[field::temperature] = inputs.solution_values[p];
+
+        for (unsigned int i = 0; i < dim; ++i)
+          {
+            computed_quantities[p][i] =
+              -thermal_conductivity_model->value(field_values) *
+              inputs.solution_gradients[p][i];
+          }
+      }
+  }
+
   std::shared_ptr<ThermalConductivityModel> thermal_conductivity_model;
-  unsigned int                              material_id;
+  const unsigned int                        mesh_material_id;
+  const bool                                solid_phase_present;
 };
 #endif

source/solvers/heat_transfer.cc

@@ -696,24 +696,37 @@ HeatTransfer<dim>::attach_solution_to_output(DataOut<dim> &data_out)
   const unsigned int n_solids =
     this->simulation_parameters.physical_properties_manager
       .get_number_of_solids();
+  unsigned int mesh_m_id = 0;
+
+  // Check if there are solids
+  const bool solid_phase_present = (n_solids > 0);
 
   // Postprocess heat fluxes
   heat_flux_postprocessors.clear();
   heat_flux_postprocessors.reserve(n_fluids + n_solids);
   // Heat fluxes in fluids
   for (unsigned int f_id = 0; f_id < n_fluids; ++f_id)
     {
-      heat_flux_postprocessors.push_back(HeatFluxPostprocessor<dim>(
-        thermal_conductivity_models[f_id], "f", f_id, f_id));
+      heat_flux_postprocessors.push_back(
+        HeatFluxPostprocessor<dim>(thermal_conductivity_models[f_id],
+                                   "f",
+                                   f_id,
+                                   mesh_m_id,
+                                   solid_phase_present));
       data_out.add_data_vector(this->dof_handler,
                                this->present_solution,
                                heat_flux_postprocessors[f_id]);
     }
   // Heat fluxes in solids
   for (unsigned int m_id = n_fluids; m_id < n_fluids + n_solids; ++m_id)
     {
-      heat_flux_postprocessors.push_back(HeatFluxPostprocessor<dim>(
-        thermal_conductivity_models[m_id], "s", m_id - n_fluids, m_id));
+      mesh_m_id += 1;
+      heat_flux_postprocessors.push_back(
+        HeatFluxPostprocessor<dim>(thermal_conductivity_models[m_id],
+                                   "s",
+                                   m_id - n_fluids,
+                                   mesh_m_id,
+                                   solid_phase_present));
       data_out.add_data_vector(this->dof_handler,
                                this->present_solution,
                                heat_flux_postprocessors[m_id]);