refactor katz2003 code, modify docs

Author: danieldouglas92

doc/sphinx/parameters/Boundary_20composition_20model.md

@@ -35,7 +35,7 @@ This parameter only describes which boundaries have a fixed composition, but not
 :name: parameters:Boundary_20composition_20model/List_20of_20model_20names
 **Default value:**
 
-**Pattern:** [MultipleSelection ascii data|box|box with lithosphere boundary indicators|function|initial composition|spherical constant ]
+**Pattern:** [MultipleSelection ascii data|box|box with lithosphere boundary indicators|function|initial composition|landlab|spherical constant ]
 
 **Documentation:** A comma-separated list of boundary composition models that will be used to initialize the composition. These plugins are loaded in the order given, and modify the existing composition field via the operators listed in ’List of model operators’.
 
@@ -57,6 +57,8 @@ The format of these functions follows the syntax understood by the muparser libr
 
 Because this class simply takes what the initial composition had described, this class can not know certain pieces of information such as the minimal and maximal composition on the boundary. For operations that require this, for example in post-processing, this boundary composition model must therefore be told what the minimal and maximal values on the boundary are. This is done using parameters set in section “Boundary composition model/Initial composition”.
 
+‘landlab’: A model in which the composition is chosen constant on the sides of a box which are selected by the parameters Left/Right/Top/Bottom/Front/Back composition
+
 ‘spherical constant’: A model in which the composition is chosen constant on the inner and outer boundaries of a sphere, spherical shell, chunk or ellipsoidal chunk. Parameters are read from subsection ’Spherical constant’.
 ::::
 
@@ -73,7 +75,7 @@ Because this class simply takes what the initial composition had described, this
 :name: parameters:Boundary_20composition_20model/Model_20name
 **Default value:** unspecified
 
-**Pattern:** [Selection ascii data|box|box with lithosphere boundary indicators|function|initial composition|spherical constant|unspecified ]
+**Pattern:** [Selection ascii data|box|box with lithosphere boundary indicators|function|initial composition|landlab|spherical constant|unspecified ]
 
 **Documentation:** Select one of the following models:
 
@@ -93,6 +95,8 @@ The format of these functions follows the syntax understood by the muparser libr
 
 Because this class simply takes what the initial composition had described, this class can not know certain pieces of information such as the minimal and maximal composition on the boundary. For operations that require this, for example in post-processing, this boundary composition model must therefore be told what the minimal and maximal values on the boundary are. This is done using parameters set in section “Boundary composition model/Initial composition”.
 
+‘landlab’: A model in which the composition is chosen constant on the sides of a box which are selected by the parameters Left/Right/Top/Bottom/Front/Back composition
+
 ‘spherical constant’: A model in which the composition is chosen constant on the inner and outer boundaries of a sphere, spherical shell, chunk or ellipsoidal chunk. Parameters are read from subsection ’Spherical constant’.
 
 **Warning**: This parameter provides an old and deprecated way of specifying boundary composition models and shouldn’t be used. Please use ’List of model names’ instead.
@@ -319,6 +323,17 @@ If the function you are describing represents a vector-valued function with mult
 **Documentation:** Minimal composition. Units: none.
 ::::
 
+(parameters:Boundary_20composition_20model/LandLab)=
+## **Subsection:** Boundary composition model / LandLab
+::::{dropdown} __Parameter:__ {ref}`Dummy variable<parameters:Boundary_20composition_20model/LandLab/Dummy_20variable>`
+:name: parameters:Boundary_20composition_20model/LandLab/Dummy_20variable
+**Default value:**
+
+**Pattern:** [List of <[Double -MAX_DOUBLE...MAX_DOUBLE (inclusive)]> of length 0...4294967295 (inclusive)]
+
+**Documentation:** A comma separated list of composition boundary values at the left boundary (at minimal $x$-value). This list must have as many entries as there are compositional fields. Units: none.
+::::
+
 (parameters:Boundary_20composition_20model/Spherical_20constant)=
 ## **Subsection:** Boundary composition model / Spherical constant
 ::::{dropdown} __Parameter:__ {ref}`Inner composition<parameters:Boundary_20composition_20model/Spherical_20constant/Inner_20composition>`

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.