Add missing assertions in rheological models

chaos-polymtl · Mar 13, 2024 · 2c105ed · 2c105ed
1 parent 079b462
commit 2c105ed
Show file tree

Hide file tree

Showing 5 changed files with 53 additions and 18 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -7,7 +7,7 @@ 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 the field is not defined, the exception thrown is comprehensible. [#add_PR_number](add_PR_link)
+- 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, an exception thrown is comprehensible. [#add_PR_number](add_PR_link)
 
 ## [Master] - 2024-03-05
 

diff --git a/include/core/parameters.h b/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

diff --git a/include/core/physical_property_model.h b/include/core/physical_property_model.h
@@ -37,8 +37,8 @@ DeclException2(PhysicialPropertyModelFieldUndefined,
                std::string,
                std::string,
                << "Error in '" << arg1 << "' model. \n "
-               << "The " << arg2
-               << " field required by the model is not defined.");
+               << "The '" << arg2
+               << "' field required by the model is not defined.");
 
 /*
  * Fields on which physical property can depend. All fields are assumed

diff --git a/include/core/rheological_model.h b/include/core/rheological_model.h
@@ -405,8 +405,7 @@ class PowerLaw : public RheologicalModel
     const std::map<field, double> &field_values) const override
   {
     AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
-                PhysicialPropertyModelFieldUndefined("PowerLaw rheological",
-                                                     "shear_rate"));
+                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;
   }
@@ -424,8 +423,7 @@ class PowerLaw : public RheologicalModel
     std::vector<double>                        &property_vector) override
   {
     AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
-                PhysicialPropertyModelFieldUndefined("PowerLaw rheological",
-                                                     "shear_rate"));
+                PhysicialPropertyModelFieldUndefined("PowerLaw", "shear_rate"));
     const std::vector<double> &shear_rate_magnitude =
       field_vectors.at(field::shear_rate);
 
@@ -571,8 +569,7 @@ class Carreau : public RheologicalModel
     const std::map<field, double> &field_values) const override
   {
     AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
-                PhysicialPropertyModelFieldUndefined("Carreau rheological",
-                                                     "shear_rate"));
+                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;
   }
@@ -590,8 +587,7 @@ class Carreau : public RheologicalModel
     std::vector<double>                        &property_vector) override
   {
     AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
-                PhysicialPropertyModelFieldUndefined("Carreau rheological",
-                                                     "shear_rate"));
+                PhysicialPropertyModelFieldUndefined("Carreau", "shear_rate"));
     const std::vector<double> &shear_rate_magnitude =
       field_vectors.at(field::shear_rate);
 

diff --git a/source/core/rheological_model.cc b/source/core/rheological_model.cc
@@ -39,6 +39,8 @@ Newtonian::value(const std::map<field, double> & /*field_values*/)
 double
 PowerLaw::value(const std::map<field, double> &field_values)
 {
+  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);
@@ -49,6 +51,8 @@ PowerLaw::vector_value(
   const std::map<field, std::vector<double>> &field_vectors,
   std::vector<double>                        &property_vector)
 {
+  AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+              PhysicialPropertyModelFieldUndefined("PowerLaw", "shear_rate"));
   const auto shear_rate_magnitude = field_vectors.at(field::shear_rate);
 
   for (unsigned int i = 0; i < shear_rate_magnitude.size(); ++i)
@@ -60,7 +64,12 @@ PowerLaw::jacobian(const std::map<field, double> &field_values, const field id)
 {
   const double shear_rate_magnitude = field_values.at(field::shear_rate);
   if (id == field::shear_rate)
-    return calculate_derivative(shear_rate_magnitude);
+    {
+      AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
+                  PhysicialPropertyModelFieldUndefined("PowerLaw",
+                                                       "shear_rate"));
+      return calculate_derivative(shear_rate_magnitude);
+    }
   else
     return 0;
 }
@@ -71,6 +80,8 @@ PowerLaw::vector_jacobian(
   const field                                 id,
   std::vector<double>                        &jacobian_vector)
 {
+  AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+              PhysicialPropertyModelFieldUndefined("PowerLaw", "shear_rate"));
   const auto shear_rate_magnitude = field_vectors.at(field::shear_rate);
 
   if (id == field::shear_rate)
@@ -83,6 +94,8 @@ PowerLaw::vector_jacobian(
 double
 Carreau::value(const std::map<field, double> &field_values)
 {
+  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);
@@ -92,6 +105,8 @@ void
 Carreau::vector_value(const std::map<field, std::vector<double>> &field_vectors,
                       std::vector<double> &property_vector)
 {
+  AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+              PhysicialPropertyModelFieldUndefined("Carreau", "shear_rate"));
   const auto shear_rate_magnitude = field_vectors.at(field::shear_rate);
 
   for (unsigned int i = 0; i < shear_rate_magnitude.size(); ++i)
@@ -105,7 +120,12 @@ double
 Carreau::jacobian(const std::map<field, double> &field_values, const field id)
 {
   if (id == field::shear_rate)
-    return this->numerical_jacobian(field_values, field::shear_rate);
+    {
+      AssertThrow(field_values.find(field::shear_rate) != field_values.end(),
+                  PhysicialPropertyModelFieldUndefined("Carreau",
+                                                       "shear_rate"));
+      return this->numerical_jacobian(field_values, field::shear_rate);
+    }
   else
     return 0;
 }
@@ -117,9 +137,14 @@ Carreau::vector_jacobian(
   std::vector<double>                        &jacobian_vector)
 {
   if (id == field::shear_rate)
-    this->vector_numerical_jacobian(field_vectors,
-                                    field::shear_rate,
-                                    jacobian_vector);
+    {
+      AssertThrow(field_vectors.find(field::shear_rate) != field_vectors.end(),
+                  PhysicialPropertyModelFieldUndefined("Carreau",
+                                                       "shear_rate"));
+      this->vector_numerical_jacobian(field_vectors,
+                                      field::shear_rate,
+                                      jacobian_vector);
+    }
   else
     std::fill(jacobian_vector.begin(), jacobian_vector.end(), 0);
 }
@@ -129,6 +154,9 @@ Carreau::vector_jacobian(
 double
 PhaseChangeRheology::value(const std::map<field, double> &field_values)
 {
+  AssertThrow(field_values.find(field::temperature) != field_values.end(),
+              PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                   "temperature"));
   const double temperature = field_values.at(field::temperature);
 
   return kinematic_viscosity(temperature);
@@ -139,6 +167,9 @@ PhaseChangeRheology::vector_value(
   const std::map<field, std::vector<double>> &field_vectors,
   std::vector<double>                        &property_vector)
 {
+  AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+              PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                   "temperature"));
   const std::vector<double> &temperature_vec =
     field_vectors.at(field::temperature);
 
@@ -153,7 +184,12 @@ PhaseChangeRheology::jacobian(const std::map<field, double> &field_values,
                               const field                    id)
 {
   if (id == field::temperature)
-    return this->numerical_jacobian(field_values, field::temperature);
+    {
+      AssertThrow(field_values.find(field::temperature) != field_values.end(),
+                  PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                       "temperature"));
+      return this->numerical_jacobian(field_values, field::temperature);
+    }
   else
     return 0;
 }
@@ -164,5 +200,8 @@ PhaseChangeRheology::vector_jacobian(
   const field                                 id,
   std::vector<double>                        &jacobian_vector)
 {
+  AssertThrow(field_vectors.find(field::temperature) != field_vectors.end(),
+              PhysicialPropertyModelFieldUndefined("PhaseChangeRheology",
+                                                   "temperature"));
   vector_numerical_jacobian(field_vectors, id, jacobian_vector);
 }