Add fields check to physical property models (chaos-polymtl#1065)

Description of the problem
Some physical and interface property models require certain fields to compute property values and partial derivatives. When these fields were not defined a not so debug-friendly message was printed out.
Description of the solution
To help the debugging process, if an error were to occur, assertions were added to the property models in functions where the fields were required.
How Has This Been Tested?
By commenting lines in unit tests and code, the throw of the exception was verified in DEBUG mode.
Comments
Doxygen descriptions were also added to specific_heat_model.h.

Former-commit-id: bc03310

CHANGELOG.md

@@ -13,6 +13,8 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
 ### Added
 
+- MINOR Assertion on fields required for physical property models were added to the code. This is done to ensure that if a physical property model necessitating a certain field is employed and that field is not defined, the exception message thrown is comprehensible. [#1065](https://github.com/lethe-cfd/lethe/pull/1065)
+
 - MINOR A new `make_table_scalars_vectors` function has been added in `utilities.h`(`.cc`). [#1062](https://github.com/lethe-cfd/lethe/pull/1062)
 
 ## [Master] - 2024-03-11

include/core/density_model.h

@@ -94,26 +94,26 @@ class DensityConstant : public DensityModel
   /**
    * @brief Compute density value.
    *
-   * @param[in] fields_value Value of the various fields on which the property
-   * may depend.
+   * @param[in] field_values Value of the various fields on which the property
+   * may depend. The map stores a single value per field.
    *
-   * @note Here, the @p fields_value parameter is ignored since the density
+   * @note Here, the @p field_values parameter is ignored since the density
    * remains constant.
    *
    * @return Density value.
    */
   double
-  value(const std::map<field, double> &fields_value) override
+  value(const std::map<field, double> &field_values) override
   {
-    (void)fields_value;
+    (void)field_values;
     return density;
   }
 
   /**
    * @brief Compute a vector of density values.
    *
    * @param[in] field_vectors Vectors of the fields on which the density may
-   * depend.
+   * depend. The map stores a vector of values per field.
    *
    * @param[out] property_vector Vector of computed density values.
    *
@@ -133,7 +133,7 @@ class DensityConstant : public DensityModel
    * respect to a specified field.
    *
    * @param[in] field_values Values of the various fields on which the density
-   * may depend.
+   * may depend. The map stores a single value per field.
    *
    * @param[in] id Indicator of the field with respect to which the jacobian
    * should be computed.
@@ -158,7 +158,7 @@ class DensityConstant : public DensityModel
    * fields.
    *
    * @param[in] field_vectors Vector of values of the fields used to evaluate
-   * the density.
+   * the density. The map stores a vector of values per field.
    *
    * @param[in] id Identifier of the field with respect to which a derivative
    * should be computed.
@@ -167,7 +167,7 @@ class DensityConstant : public DensityModel
    * with respect to the field of the specified @p id. In this case, it returns
    * a vector of zeros since the density remains constant.
    *
-   * @note Here, the @p field_values and @p id parameters are ignored since the
+   * @note Here, the @p field_vectors and @p id parameters are ignored since the
    * density remains constant.
    *
    */
@@ -256,30 +256,38 @@ class DensityIsothermalIdealGas : public DensityModel
   /**
    * @brief Compute the density of the isothermal ideal gas.
    *
-   * @param[in] fields_value Value of the various fields on which the property
-   * may depend. In this case, the density depends on pressure.
+   * @param[in] field_values Values of the various fields on which the property
+   * may depend. In this case, the density depends on pressure. The map stores a
+   * single value per field.
    *
-   * @return Value of the density computed with the @p fields_value.
+   * @return Value of the density computed with the @p field_values.
    */
   double
-  value(const std::map<field, double> &fields_value) override
+  value(const std::map<field, double> &field_values) override
   {
-    const double pressure = fields_value.at(field::pressure);
+    AssertThrow(field_values.find(field::pressure) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "DensityIsothermalIdealGas", "pressure"));
+    const double pressure = field_values.at(field::pressure);
     return density_ref + psi * pressure;
   }
 
   /**
    * @brief Compute a vector of density values for an isothermal ideal gas.
    *
    * @param[in] field_vectors Vectors of the fields on which the density may
-   * depend. In this case, the density depends on pressure.
+   * depend. In this case, the density depends on pressure. The map stores a
+   * vector of values per field.
    *
    * @param[out] property_vector Vectors of computed density values.
    */
   void
   vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::pressure) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "DensityIsothermalIdealGas", "pressure"));
     const std::vector<double> &pressure = field_vectors.at(field::pressure);
     for (unsigned int i = 0; i < property_vector.size(); ++i)
       property_vector[i] = density_ref + psi * pressure[i];
@@ -290,7 +298,7 @@ class DensityIsothermalIdealGas : public DensityModel
    * respect to a specified field.
    *
    * @param[in] field_values Values of the various fields on which the density
-   * may depend.
+   * may depend. The map stores a single value per field.
    *
    * @param[in] id Indicator of the field with respect to which the jacobian
    * should be computed.
@@ -313,7 +321,7 @@ class DensityIsothermalIdealGas : public DensityModel
    * an isothermal ideal gas.
    *
    * @param[in] field_vectors Vector of values of the fields used to evaluate
-   * the density.
+   * the density. The map stores a vector of values per field.
    *
    * @param[in] id Identifier of the field with respect to which a derivative
    * should be computed.

include/core/evaporation_model.h

@@ -215,7 +215,7 @@ class EvaporationModelConstant : public EvaporationModel
     , ambient_gas_density_inv(1.0 / p_evaporation.ambient_gas_density)
     , liquid_density_inv(1.0 / p_evaporation.liquid_density)
   {
-    model_depends_on[temperature] = false;
+    model_depends_on[field::temperature] = false;
   }
 
 
@@ -413,7 +413,7 @@ class EvaporationModelTemperature : public EvaporationModel
     , liquid_density_inv(1.0 / p_evaporation.liquid_density)
     , universal_gas_constant(p_evaporation.universal_gas_constant)
   {
-    model_depends_on[temperature] = true;
+    model_depends_on[field::temperature] = true;
   }
 
   /**
@@ -425,6 +425,9 @@ class EvaporationModelTemperature : public EvaporationModel
   double
   saturation_pressure(const std::map<field, double> &field_values)
   {
+    AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "EvaporationModelTemperature", "temperature"));
     const double temperature_inv =
       1.0 / (field_values.at(field::temperature) + 1e-16);
 
@@ -450,6 +453,9 @@ class EvaporationModelTemperature : public EvaporationModel
   saturation_pressure(const std::map<field, std::vector<double>> &field_vectors,
                       std::vector<double> &saturation_pressure_vector)
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "EvaporationModelTemperature", "temperature"));
     const std::vector<double> &temperature =
       field_vectors.at(field::temperature);
 
@@ -506,6 +512,9 @@ class EvaporationModelTemperature : public EvaporationModel
   mass_flux(const std::map<field, std::vector<double>> &field_vectors,
             std::vector<double> &mass_flux_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "EvaporationModelTemperature", "temperature"));
     const std::vector<double> &temperature =
       field_vectors.at(field::temperature);
 
@@ -578,6 +587,9 @@ class EvaporationModelTemperature : public EvaporationModel
   heat_flux_jacobian(const std::map<field, double> &field_values,
                      const field                    id) override
   {
+    AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined(
+                  "EvaporationModelTemperature", "temperature"));
     const double temperature_inv =
       1.0 / (field_values.at(field::temperature) + 1e-16);
 
@@ -620,6 +632,10 @@ class EvaporationModelTemperature : public EvaporationModel
 
     if (id == field::temperature)
       {
+        AssertThrow(
+          field_vectors.find(field::temperature) != field_vectors.end(),
+          PhysicialPropertyModelFieldUndefined("EvaporationModelTemperature",
+                                               "temperature"));
         const std::vector<double> &temperature =
           field_vectors.at(field::temperature);
 

include/core/mobility_cahn_hilliard_model.h

@@ -199,6 +199,10 @@ class MobilityCahnHilliardModelQuartic : public MobilityCahnHilliardModel
   double
   value(const std::map<field, double> &fields_value) override
   {
+    AssertThrow(
+      fields_value.find(field::phase_order_cahn_hilliard) != fields_value.end(),
+      PhysicialPropertyModelFieldUndefined("MobilityCahnHilliardModelQuartic",
+                                           "phase_order_cahn_hilliard"));
     const double &phase_order_cahn_hilliard =
       fields_value.at(field::phase_order_cahn_hilliard);
 
@@ -221,6 +225,11 @@ class MobilityCahnHilliardModelQuartic : public MobilityCahnHilliardModel
   vector_value(const std::map<field, std::vector<double>> &field_vectors,
                std::vector<double> &property_vector) override
   {
+    AssertThrow(
+      field_vectors.find(field::phase_order_cahn_hilliard) !=
+        field_vectors.end(),
+      PhysicialPropertyModelFieldUndefined("MobilityCahnHilliardModelQuartic",
+                                           "phase_order_cahn_hilliard"));
     const std::vector<double> &phase_order_cahn_hilliard =
       field_vectors.at(field::phase_order_cahn_hilliard);
     for (unsigned int i = 0; i < property_vector.size(); ++i)
@@ -248,6 +257,10 @@ class MobilityCahnHilliardModelQuartic : public MobilityCahnHilliardModel
   double
   jacobian(const std::map<field, double> &fields_value, field /*id*/) override
   {
+    AssertThrow(
+      fields_value.find(field::phase_order_cahn_hilliard) != fields_value.end(),
+      PhysicialPropertyModelFieldUndefined("MobilityCahnHilliardModelQuartic",
+                                           "phase_order_cahn_hilliard"));
     const double &phase_order_cahn_hilliard =
       fields_value.at(field::phase_order_cahn_hilliard);
 
@@ -275,6 +288,11 @@ class MobilityCahnHilliardModelQuartic : public MobilityCahnHilliardModel
                   const field /*id*/,
                   std::vector<double> &jacobian_vector) override
   {
+    AssertThrow(
+      field_vectors.find(field::phase_order_cahn_hilliard) !=
+        field_vectors.end(),
+      PhysicialPropertyModelFieldUndefined("MobilityCahnHilliardModelQuartic",
+                                           "phase_order_cahn_hilliard"));
     const std::vector<double> &phase_order_cahn_hilliard =
       field_vectors.at(field::phase_order_cahn_hilliard);
     for (unsigned int i = 0; i < jacobian_vector.size(); ++i)

include/core/parameters.h

@@ -254,7 +254,7 @@ namespace Parameters
   };
 
   /**
-   * @brief Carreau rheological model to solve for non Newtonian
+   * @brief Carreau model to solve for non Newtonian
    * flows.
    */
   struct CarreauParameters

include/core/physical_property_model.h

@@ -33,6 +33,13 @@ DeclException2(SizeOfFields,
                << " is not equal to the number of values for another field "
                << arg2);
 
+DeclException2(PhysicialPropertyModelFieldUndefined,
+               std::string,
+               std::string,
+               << "Error in '" << arg1 << "' model. \n "
+               << "The '" << arg2
+               << "' field required by the model is not defined.");
+
 /*
  * Fields on which physical property can depend. All fields are assumed
  * to be at time t+dt other than those for which a _p suffix is explicitly

include/core/rheological_model.h

@@ -300,8 +300,8 @@ class PowerLaw : public RheologicalModel
     , n(n)
     , shear_rate_min(shear_rate_min)
   {
-    this->model_depends_on[shear_rate]    = true;
-    this->non_newtonian_rheological_model = true;
+    this->model_depends_on[field::shear_rate] = true;
+    this->non_newtonian_rheological_model     = true;
   }
 
   /**
@@ -404,6 +404,8 @@ class PowerLaw : public RheologicalModel
     const double                  &p_density_ref,
     const std::map<field, double> &field_values) const override
   {
+    AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined("PowerLaw", "shear_rate"));
     const double shear_rate_magnitude = field_values.at(field::shear_rate);
     return calculate_kinematic_viscosity(shear_rate_magnitude) * p_density_ref;
   }
@@ -420,6 +422,8 @@ class PowerLaw : public RheologicalModel
     const std::map<field, std::vector<double>> &field_vectors,
     std::vector<double>                        &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined("PowerLaw", "shear_rate"));
     const std::vector<double> &shear_rate_magnitude =
       field_vectors.at(field::shear_rate);
 
@@ -469,8 +473,8 @@ class Carreau : public RheologicalModel
     , a(p_a)
     , n(p_n)
   {
-    this->model_depends_on[shear_rate]    = true;
-    this->non_newtonian_rheological_model = true;
+    this->model_depends_on[field::shear_rate] = true;
+    this->non_newtonian_rheological_model     = true;
   }
 
   /**
@@ -564,6 +568,8 @@ class Carreau : public RheologicalModel
     const double                  &p_density_ref,
     const std::map<field, double> &field_values) const override
   {
+    AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined("Carreau", "shear_rate"));
     const double shear_rate_magnitude = field_values.at(field::shear_rate);
     return calculate_kinematic_viscosity(shear_rate_magnitude) * p_density_ref;
   }
@@ -580,6 +586,8 @@ class Carreau : public RheologicalModel
     const std::map<field, std::vector<double>> &field_vectors,
     std::vector<double>                        &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined("Carreau", "shear_rate"));
     const std::vector<double> &shear_rate_magnitude =
       field_vectors.at(field::shear_rate);
 
@@ -639,7 +647,7 @@ class PhaseChangeRheology : public RheologicalModel
   PhaseChangeRheology(const Parameters::PhaseChange p_phase_change_parameters)
     : param(p_phase_change_parameters)
   {
-    this->model_depends_on[temperature] = true;
+    this->model_depends_on[field::temperature] = true;
   }
 
   /**
@@ -700,6 +708,9 @@ class PhaseChangeRheology : public RheologicalModel
     const double                  &p_density_ref,
     const std::map<field, double> &field_values) const override
   {
+    AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                     "temperature"));
     const double temperature = field_values.at(field::temperature);
     return kinematic_viscosity(temperature) * p_density_ref;
   }
@@ -716,6 +727,9 @@ class PhaseChangeRheology : public RheologicalModel
     const std::map<field, std::vector<double>> &field_vectors,
     std::vector<double>                        &property_vector) override
   {
+    AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+                PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                     "temperature"));
     const std::vector<double> &temperature =
       field_vectors.at(field::temperature);