add some asserts

Author: danieldouglas92

include/aspect/material_model/rheology/compositional_plasticity_prefactors.h

@@ -34,8 +34,8 @@ namespace aspect
     namespace Rheology
     {
       /**
-       * A class that handles multiplication of viscosity for a given compositional
-       * field. The multiplication factors for each composition (viscosity
+       * A class that handles multiplication of the plastic parameters for a given compositional
+       * field. The multiplication factors for each composition (plasticity
        * prefactors) are also declared, parsed, and in some cases calculated in this class.
        */
       template <int dim>
@@ -61,7 +61,7 @@ namespace aspect
           parse_parameters (ParameterHandler &prm);
 
           /**
-           * Compute the viscosity.
+           * Compute the plasticity weakening factor.
            */
           std::pair<double, double>
           compute_weakening_factor (const MaterialModel::MaterialModelInputs<dim> &in,
@@ -73,18 +73,16 @@ namespace aspect
 
         private:
           /**
-           * The viscosity prefactors or terms used to calculate the viscosity
+           * The plasticity prefactors or terms used to calculate the plasticity
            * prefactors, which are read in from the input file by the
            * parse_parameters() function. Users can choose between different schemes.
-           * none: no viscosity change
-           * hk04_olivine_hydration: calculate the viscosity change due to hydrogen
-           * incorporation into olivine using Hirth & Kohlstaedt 2004 10.1029/138GM06.
-           * This method requires a composition called 'bound_fluid' which tracks the wt%
-           * water in the solid, which is used to compute an atomic ratio of H/Si ppm
-           * assuming 90 mol% forsterite and 10 mol% fayalite, and finally calculates
-           * a water fugacity.
+           * none: no change in the plasticity parameters
+           * porosity: calculate the change in the plasticity parameters due to the amount
+           * of porosity at each quadrature point.
+           * function: calculate the change in the plasticity parameters by specifying the
+           * prefactors via a function.
            * The prefactor for a given compositional field is multiplied with a
-           * base_viscosity value provided by the material model, which is then returned
+           * cohesions and the friction angles provided by the material model, which is then returned
            * to the material model.
            */
           enum WeakeningMechanism
@@ -100,13 +98,13 @@ namespace aspect
           std::vector<double> minimum_cohesions;
           std::vector<double> minimum_friction_angles;
           /**
-           * Parsed functions that specify the friction angle which must be
+           * Parsed functions that specify the plasticity prefactors which must be
            * given in the input file using the function method.
            */
           std::unique_ptr<Functions::ParsedFunction<dim>> prefactor_function;
 
           /**
-           * The coordinate representation to evaluate the function for the friction angle.
+           * The coordinate representation to evaluate the function for the plasticity prefactors.
            * Possible choices are depth, cartesian and spherical.
            */
           Utilities::Coordinates::CoordinateSystem coordinate_system_prefactor_function;

source/material_model/rheology/compositional_plasticity_prefactors.cc

@@ -61,7 +61,6 @@ namespace aspect
             {
               // Use the porosity to compute the cohesion and friction angle
               // The porosity is computed in the MaterialModel::MaterialModelInputs
-              // and is available as in.porosity[i]
               const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
               plasticity_parameters.first *= std::max((max_porosities_for_cohesion_prefactors[volume_fraction_index] - in.composition[i][porosity_idx]) / max_porosities_for_cohesion_prefactors[volume_fraction_index], 0.0);
               plasticity_parameters.second *= std::max((max_porosities_for_friction_angle_prefactors[volume_fraction_index] - in.composition[i][porosity_idx]) / max_porosities_for_friction_angle_prefactors[volume_fraction_index], 0.0);
@@ -72,7 +71,7 @@ namespace aspect
             }
             case function:
             {
-              // Use a given function input per composition to get the friction angle
+              // Use a given function input per composition to get the prefactors
               Utilities::NaturalCoordinate<dim> point =
                 this->get_geometry_model().cartesian_to_other_coordinates(position, coordinate_system_prefactor_function);
 
@@ -83,7 +82,7 @@ namespace aspect
               else
                 prefactor_function->set_time (this->get_time());
 
-              // determine the friction angle based on position and composition
+              // determine the prefactors based on position and composition
               // The volume_fraction_index is based on the number of chemical compositional fields.
               // However, this plugin reads a function for every compositional field, regardless of
               // its type. Therefore we have to get the correct index.
@@ -96,6 +95,7 @@ namespace aspect
               double prefactor_from_function =
                 prefactor_function->value(Utilities::convert_array_to_point<dim>(point.get_coordinates()),index);
 
+              // Multiply the cohesion and the friction angle by the prefactor.
               plasticity_parameters.first *= prefactor_from_function;
               plasticity_parameters.second *= prefactor_from_function * constants::degree_to_radians;
 
@@ -110,14 +110,14 @@ namespace aspect
       CompositionalPlasticityPrefactors<dim>::declare_parameters (ParameterHandler &prm)
       {
         prm.declare_entry ("Maximum porosity for cohesion prefactor", "1.0",
-                           Patterns::List(Patterns::Double(0., 1.0)),
+                           Patterns::List(Patterns::Double(std::numeric_limits<double>::epsilon(), 1.0)),
                            "List of the maximum porosity value for calculating the prefactor for the cohesion. "
                            "entered as a volume fraction. If chosen to be 0.1 (10 percent porosity), the "
                            "cohesion will be linearly reduced from its default value when the model porosity is 0"
                            "to a maximum reduction when the model porosity is 0.1. Units: none.");
 
         prm.declare_entry ("Maximum porosity for friction angle prefactor", "1.0",
-                           Patterns::List(Patterns::Double(0., 1.0)),
+                           Patterns::List(Patterns::Double(std::numeric_limits<double>::epsilon(), 1.0)),
                            "List of the maximum porosity value for calculating the prefactor for the friction "
                            "angle, entered as a volume fraction. If chosen to be 0.1 (10 percent porosity), the "
                            "friction angle will be linearly reduced from its default value when the model porosity is 0"
@@ -149,12 +149,14 @@ namespace aspect
       CompositionalPlasticityPrefactors<dim>::parse_parameters (ParameterHandler &prm)
       {
 
-        // if friction is specified as a function
         const unsigned int n_fields = this->n_compositional_fields() + 1;
         if (prm.get ("Plasticity prefactor scheme") == "none")
           prefactor_mechanism = none;
-        if (prm.get ("Plasticity prefactor scheme") == "porosity")
+        else if (prm.get ("Plasticity prefactor scheme") == "porosity")
           {
+            AssertThrow(this->introspection().compositional_name_exists("porosity"),
+                        ExcMessage("The 'porosity' Plasticity prefactor scheme work only if "
+                                   "is a compositional field called porosity."));
             prefactor_mechanism = porosity;
 
             // Retrieve the list of chemical names
@@ -178,7 +180,7 @@ namespace aspect
             minimum_friction_angles = Utilities::MapParsing::parse_map_to_double_array (prm.get("Minimum friction angles"),
                                                                                         options);
           }
-        if (prm.get ("Plasticity prefactor scheme") == "function")
+        else if (prm.get ("Plasticity prefactor scheme") == "function")
           {
             prefactor_mechanism = function;
             prm.enter_subsection("Prefactor function");
@@ -201,6 +203,8 @@ namespace aspect
             }
             prm.leave_subsection();
           }
+        else
+          AssertThrow(false, ExcMessage("Not a valid plasticity prefactor scheme"));
       }
     }
   }

tests/test_plastic_prefactor.prm

@@ -131,7 +131,7 @@ subsection Material model
     set Cohesions                                      = 100e6 
     set Angles of internal friction                    = 30
     set Plasticity prefactor scheme                    = porosity
-    set Maximum porosity for friction angle prefactor  = 0.2
+    set Maximum porosity for friction angle prefactor  = 0.
     set Maximum porosity for cohesion prefactor        = 0.1
 
     set Minimum cohesions                              = 10e6