refactor katz2003 code, modify docs

Author: danieldouglas92

doc/sphinx/parameters/Material_20model.md

@@ -3821,6 +3821,15 @@ Also note that the fluid reaction time scale has to be larger than or equal to t
 **Documentation:** Lower cutoff for the compaction viscosity. Units: $\text{Pa}\text{s}$.
 ::::
 
+::::{dropdown} __Parameter:__ {ref}`Minimum weakened viscosity<parameters:Material_20model/Reactive_20Fluid_20Transport_20Model/Minimum_20weakened_20viscosity>`
+:name: parameters:Material_20model/Reactive_20Fluid_20Transport_20Model/Minimum_20weakened_20viscosity
+**Default value:** 1e17
+
+**Pattern:** [Double 0...MAX_DOUBLE (inclusive)]
+
+**Documentation:** Lower cutoff for the viscosity of the solid phase in the presence of fluids. Units: $\text{Pa}\text{s}$.
+::::
+
 ::::{dropdown} __Parameter:__ {ref}`Reference fluid density<parameters:Material_20model/Reactive_20Fluid_20Transport_20Model/Reference_20fluid_20density>`
 :name: parameters:Material_20model/Reactive_20Fluid_20Transport_20Model/Reference_20fluid_20density
 **Default value:** 2500

source/material_model/reaction_model/katz2003_mantle_melting.cc

@@ -344,15 +344,6 @@ namespace aspect
                 melt_out->compaction_viscosities[i] *= visc_temperature_dependence;
               }
           }
-
-        if (this->include_melt_transport() && in.requests_property(MaterialProperties::viscosity))
-          {
-            for (unsigned int i=0; i<in.n_evaluation_points(); ++i)
-              {
-                const double porosity = std::min(1.0, std::max(in.composition[i][porosity_idx],0.0));
-                out.viscosities[i] *= std::exp(- alpha_phi * porosity);
-              }
-          }
       }
 
 

source/material_model/reactive_fluid_transport.cc

@@ -180,22 +180,22 @@ namespace aspect
     {
       base_model->evaluate(in,out);
 
-      if (fluid_solid_reaction_scheme != katz2003)
-        {
-          const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
+      const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
 
-          // Modify the viscosity from the base model based on the presence of fluid.
-          if (in.requests_property(MaterialProperties::viscosity))
+      // Modify the viscosity from the base model based on the presence of fluid.
+      if (in.requests_property(MaterialProperties::viscosity))
+        {
+          // Scale the base model viscosity value based on the porosity.
+          for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
             {
-              // Scale the base model viscosity value based on the porosity.
-              for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
-                {
-                  const double porosity = std::max(in.composition[q][porosity_idx],0.0);
-                  const double weakened_viscosity = out.viscosities[q] * (1.0 - porosity) * std::exp(- alpha_phi * porosity);
-                  out.viscosities[q] = std::max(weakened_viscosity, min_weakened_viscosity);
-                }
+              const double porosity = std::max(in.composition[q][porosity_idx],0.0);
+              const double weakened_viscosity = out.viscosities[q] * (1.0 - porosity) * std::exp(- alpha_phi * porosity);
+              out.viscosities[q] = std::max(weakened_viscosity, min_weakened_viscosity);
             }
+        }
 
+      if (fluid_solid_reaction_scheme != katz2003)
+        {
           // Fill the melt outputs if they exist. Note that the MeltOutputs class was originally
           // designed for two-phase flow material models in ASPECT that model the flow of melt,
           // but can be reused for a geofluid of arbitrary composition.