diff --git a/include/aspect/material_model/rheology/composite_visco_plastic.h b/include/aspect/material_model/rheology/composite_visco_plastic.h
index 4ffd7a5a97e..523142be801 100644
--- a/include/aspect/material_model/rheology/composite_visco_plastic.h
+++ b/include/aspect/material_model/rheology/composite_visco_plastic.h
@@ -38,6 +38,31 @@ namespace aspect
 
     namespace Rheology
     {
+      /**
+       * A namespace for selecting how to average the viscosity
+       * over multiple chemical compositions.
+       *
+       * Current options are:
+       *  'isostress': Reuss (harmonic) averaging of viscosities.
+       *  'isostrain': Voigt (arithmetic) averaging of viscosities.
+       */
+      namespace ViscosityAveraging
+      {
+        enum Kind
+        {
+          isostress,
+          isostrain
+        };
+
+        /**
+         * Read the lithosphere age model from the parameter file,
+         * using the parameter name given in @p parameter_name, and return the
+         * enum that corresponds to this operation.
+         */
+        ViscosityAveraging::Kind
+        parse (const std::string &parameter_name,
+               const ParameterHandler &prm);
+      }
 
       template <int dim>
       class CompositeViscoPlastic : public ::aspect::SimulatorAccess<dim>
@@ -83,6 +108,7 @@ namespace aspect
                              const std::vector<double> &phase_function_values = std::vector<double>(),
                              const std::vector<unsigned int> &n_phase_transitions_per_composition = std::vector<unsigned int>()) const;
 
+        private:
           /**
            * Compute the isostress viscosity over all compositional fields
            * based on the composite viscous creep law.
@@ -104,7 +130,7 @@ namespace aspect
           /**
            * Compute the total strain rate and the first derivative of log strain rate
            * with respect to log viscoplastic stress from individual log strain rate components
-           * over all compositional fields. Also updates the partial_strain_rates vector
+           * over all compositional fields. Also updates the partial_strain_rates vector.
            */
           std::pair<double, double>
           calculate_isostress_log_strain_rate_and_derivative(const std::vector<std::array<std::pair<double, double>, 4>> &logarithmic_strain_rates_and_stress_derivatives,
@@ -112,7 +138,26 @@ namespace aspect
                                                              std::vector<double> &partial_strain_rates) const;
 
           /**
-           * Compute the compositional field viscosity
+           * Compute the isostrain viscosity over all compositional fields
+           * based on the composite viscous creep law.
+           * If @p expected_n_phases_per_composition points to a vector of
+           * unsigned integers this is considered the number of phase transitions
+           * for each compositional field and viscosity will be first computed on
+           * each phase and then averaged for each compositional field.
+           */
+          double
+          compute_isostrain_viscosity (const double pressure,
+                                       const double temperature,
+                                       const double grain_size,
+                                       const std::vector<double> &volume_fractions,
+                                       const SymmetricTensor<2,dim> &strain_rate,
+                                       std::vector<double> &partial_strain_rates,
+                                       const std::vector<double> &phase_function_values = std::vector<double>(),
+                                       const std::vector<unsigned int> &n_phase_transitions_per_composition = std::vector<unsigned int>()) const;
+
+
+          /**
+           * Compute the viscosity for a single compositional field
            * based on the composite viscous creep law.
            * If @p n_phase_transitions_per_composition points to a vector of
            * unsigned integers this is considered the number of phase transitions
@@ -131,26 +176,22 @@ namespace aspect
 
           /**
            * Compute the total strain rate and the first derivative of log strain rate
-           * with respect to log viscoplastic stress from individual log strain rate components.
-           * Also updates the partial_strain_rates vector
+           * with respect to log viscoplastic stress from individual log strain rate components
+           * for a single compositional field.
+           * Also updates the partial_strain_rates vector.
            */
           std::pair<double, double>
           calculate_composition_log_strain_rate_and_derivative(const std::array<std::pair<double, double>, 4> &logarithmic_strain_rates_and_stress_derivatives,
                                                                const double viscoplastic_stress,
                                                                std::vector<double> &partial_strain_rates) const;
 
-        private:
           /**
            * Enumeration for selecting which type of viscosity averaging to use.
            */
-          enum ViscosityAveragingScheme
-          {
-            isostrain,
-            isostress
-          } viscosity_averaging_scheme;
+          ViscosityAveraging::Kind viscosity_averaging_scheme;
 
           /**
-           * Whether to use different deformation mechanisms
+           * Whether to use different deformation mechanisms.
            */
           bool use_diffusion_creep;
           bool use_dislocation_creep;
@@ -158,10 +199,18 @@ namespace aspect
           bool use_drucker_prager;
 
           /**
-           * Vector storing which flow mechanisms are active
+           * Vector storing which flow mechanisms are active.
            */
           std::vector<unsigned int> active_flow_mechanisms;
 
+          /**
+           * Total number of decomposed strain rates. This includes the maximum
+           * viscosity limiter but not the minimum viscosity limiter,
+           * which is arranged in parallel with the viscoplastic elements and
+           * therefore does not contribute to the total strain rate.
+           */
+          static constexpr unsigned int n_decomposed_strain_rates = 5;
+
           /**
            * Pointers to objects for computing deformation mechanism
            * strain rates and effective viscosities.
@@ -178,7 +227,7 @@ namespace aspect
 
           /**
            * The maximum viscosity, imposed via an isoviscous damper
-           * in series with the composite viscoplastic element
+           * in series with the composite viscoplastic element.
            */
           double maximum_viscosity;
 
@@ -221,7 +270,7 @@ namespace aspect
           unsigned int stress_max_iteration_number;
 
           /**
-           * Useful number to aid in adding together exponentials
+           * Useful number to aid in adding together exponentials.
            */
           const double logmin = std::log(std::numeric_limits<double>::min());
       };
diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/source/material_model/rheology/composite_visco_plastic.cc
index 4bffd3763c7..7ce830c0393 100644
--- a/source/material_model/rheology/composite_visco_plastic.cc
+++ b/source/material_model/rheology/composite_visco_plastic.cc
@@ -62,6 +62,29 @@ namespace aspect
       //}
 
 
+
+      namespace ViscosityAveraging
+      {
+        ViscosityAveraging::Kind
+        parse (const std::string &parameter_name,
+               const ParameterHandler &prm)
+        {
+          ViscosityAveraging::Kind scheme;
+          if (prm.get (parameter_name) == "isostress")
+            scheme = isostress;
+          else if (prm.get (parameter_name) == "isostrain")
+            scheme = isostrain;
+          else
+            {
+              AssertThrow(false, ExcMessage("Not a valid viscosity averaging scheme. Choose either isostress or isostrain."));
+            }
+
+          return scheme;
+        }
+      }
+
+
+
       template <int dim>
       CompositeViscoPlastic<dim>::CompositeViscoPlastic ()
         = default;
@@ -80,61 +103,48 @@ namespace aspect
                                                      const std::vector<unsigned int> &n_phase_transitions_per_composition) const
       {
         double viscosity = 0.;
-        partial_strain_rates.resize(5, 0.);
 
-        // Isostress averaging
-        if (viscosity_averaging_scheme == isostress)
-          {
-            viscosity += compute_isostress_viscosity (pressure,
-                                                      temperature,
-                                                      grain_size,
-                                                      volume_fractions,
-                                                      strain_rate,
-                                                      partial_strain_rates,
-                                                      phase_function_values,
-                                                      n_phase_transitions_per_composition);
-          }
+        // In case partial_strain_rates already contains non-zero values,
+        // fill with zeros and then resize to the right length
+        std::fill(partial_strain_rates.begin(), partial_strain_rates.end(), 0);
+        partial_strain_rates.resize(n_decomposed_strain_rates, 0.);
 
-        // Isostrain averaging
-        if (viscosity_averaging_scheme == isostrain)
+        switch (viscosity_averaging_scheme)
           {
-            double total_volume_fraction = 1.;
-            for (unsigned int composition=0; composition < number_of_chemical_compositions; ++composition)
-              {
-                // Only include the contribution to the viscosity
-                // from a given composition if the volume fraction exceeds
-                // a certain (small) fraction.
-                if (volume_fractions[composition] > 2.*std::numeric_limits<double>::epsilon())
-                  {
-                    std::vector<double> partial_strain_rates_composition(5, 0.);
-                    viscosity += (volume_fractions[composition]
-                                  * compute_composition_viscosity (pressure,
-                                                                   temperature,
-                                                                   grain_size,
-                                                                   composition,
-                                                                   strain_rate,
-                                                                   partial_strain_rates_composition,
-                                                                   phase_function_values,
-                                                                   n_phase_transitions_per_composition));
-                    for (unsigned int j=0; j < 5; ++j)
-                      partial_strain_rates[j] += volume_fractions[composition] * partial_strain_rates_composition[j];
-                  }
-                else
-                  {
-                    total_volume_fraction -= volume_fractions[composition];
-                  }
-
-                viscosity /= total_volume_fraction;
-                for (unsigned int j=0; j < 5; ++j)
-                  partial_strain_rates[j] /= total_volume_fraction;
-              }
+            case ViscosityAveraging::Kind::isostress:
+            {
+              viscosity = compute_isostress_viscosity (pressure,
+                                                       temperature,
+                                                       grain_size,
+                                                       volume_fractions,
+                                                       strain_rate,
+                                                       partial_strain_rates,
+                                                       phase_function_values,
+                                                       n_phase_transitions_per_composition);
+              break;
+            }
+            case ViscosityAveraging::Kind::isostrain:
+            {
+              viscosity = compute_isostrain_viscosity (pressure,
+                                                       temperature,
+                                                       grain_size,
+                                                       volume_fractions,
+                                                       strain_rate,
+                                                       partial_strain_rates,
+                                                       phase_function_values,
+                                                       n_phase_transitions_per_composition);
+              break;
+            }
+            default:
+            {
+              Assert (false, ExcNotImplemented());
+            }
           }
         return viscosity;
       }
 
 
 
-
       template <int dim>
       double
       CompositeViscoPlastic<dim>::compute_isostress_viscosity (const double pressure,
@@ -228,6 +238,19 @@ namespace aspect
 
         for (unsigned int i = 0; i < active_compositions.size(); ++i)
           {
+            // In the isostress viscosity averaging formalism,
+            // the total strain rate is the sum of the strain rates for each
+            // composition i multiplied by the volume fraction of that
+            // composition. The strain rate for each composition is the sum
+            // of strain rates for each deformation mechanism j:
+            // edot_total = sum_i (volume_fraction_i * (sum_j (edot_ij)))
+            // This can be refactored:
+            // edot_total = sum_i (sum_j (volume_fraction_i * (edot_ij)))
+            // Therefore the logarithm of the contribution of each
+            // component-deformation mechanism to the total strain rate
+            // is:
+            // ln(edot_total_contrib_ij) = ln(volume_fraction_i) + ln(edot_ij)
+
             const double log_volume_fraction = std::log(volume_fractions[active_compositions[i]]);
 
             if (use_diffusion_creep)
@@ -269,7 +292,9 @@ namespace aspect
         // between these components.
 
         // The following while loop contains a Newton iteration to obtain the
-        // viscoplastic stress that is consistent with the total strain rate.
+        // viscoplastic stress that is consistent with the total strain rate
+        // and results in the correct partitioning of strain rate amongst
+        // the different flow components.
         unsigned int stress_iteration = 0;
         while (std::abs(log_strain_rate_residual) > log_strain_rate_residual_threshold && stress_iteration < stress_max_iteration_number)
           {
@@ -382,6 +407,57 @@ namespace aspect
 
 
 
+      template <int dim>
+      double
+      CompositeViscoPlastic<dim>::compute_isostrain_viscosity (const double pressure,
+                                                               const double temperature,
+                                                               const double grain_size,
+                                                               const std::vector<double> &volume_fractions,
+                                                               const SymmetricTensor<2,dim> &strain_rate,
+                                                               std::vector<double> &partial_strain_rates,
+                                                               const std::vector<double> &phase_function_values,
+                                                               const std::vector<unsigned int> &n_phase_transitions_per_composition) const
+      {
+        // The isostrain viscosity is calculated by summing the viscosities
+        //of each composition weighted by their volume fractions.
+        double viscosity = 0.;
+        double total_volume_fraction = 1.;
+        for (unsigned int composition=0; composition < number_of_chemical_compositions; ++composition)
+          {
+            // Only include the contribution to the viscosity
+            // from a given composition if the volume fraction exceeds
+            // a certain (small) fraction.
+            if (volume_fractions[composition] > 2. * std::numeric_limits<double>::epsilon())
+              {
+                std::vector<double> partial_strain_rates_composition(n_decomposed_strain_rates, 0.);
+                viscosity += (volume_fractions[composition]
+                              * compute_composition_viscosity (pressure,
+                                                               temperature,
+                                                               grain_size,
+                                                               composition,
+                                                               strain_rate,
+                                                               partial_strain_rates_composition,
+                                                               phase_function_values,
+                                                               n_phase_transitions_per_composition));
+                for (unsigned int j=0; j < n_decomposed_strain_rates; ++j)
+                  partial_strain_rates[j] += volume_fractions[composition] * partial_strain_rates_composition[j];
+              }
+            else
+              {
+                total_volume_fraction -= volume_fractions[composition];
+              }
+          }
+
+        viscosity /= total_volume_fraction;
+
+        for (unsigned int j=0; j < n_decomposed_strain_rates; ++j)
+          partial_strain_rates[j] /= total_volume_fraction;
+
+        return viscosity;
+      }
+
+
+
       template <int dim>
       double
       CompositeViscoPlastic<dim>::compute_composition_viscosity (const double pressure,
@@ -393,6 +469,7 @@ namespace aspect
                                                                  const std::vector<double> &phase_function_values,
                                                                  const std::vector<unsigned int> &n_phase_transitions_per_composition) const
       {
+        // Calculate the viscosity for a single compositional field based on the local state.
         // If strain rate is zero (like during the first time step) set it to some very small number
         // to prevent a division-by-zero, and a floating point exception.
         // Otherwise, calculate the square-root of the norm of the second invariant of the deviatoric-
@@ -537,7 +614,8 @@ namespace aspect
           const double viscoplastic_stress,
           std::vector<double> &partial_strain_rates) const
       {
-        // The total strain rate
+        // Initialize the double containing the total strain rate
+        // for the single compositional field of interest.
         double viscoplastic_strain_rate_sum = 0.0;
 
         // The sum of the stress derivatives multiplied by the mechanism strain rates
@@ -593,7 +671,7 @@ namespace aspect
       CompositeViscoPlastic<dim>::declare_parameters (ParameterHandler &prm)
       {
         prm.declare_entry ("Viscosity averaging scheme", "isostress",
-                           Patterns::Selection("isostress|isostrain|unspecified"),
+                           Patterns::Selection("isostress|isostrain"),
                            "Determines the relationship between the conditions experienced by "
                            "each chemical compositional field in a composite material. "
                            "Select either isostress (the default option, Reuss averaging) "
@@ -661,14 +739,8 @@ namespace aspect
       CompositeViscoPlastic<dim>::parse_parameters (ParameterHandler &prm,
                                                     const std::unique_ptr<std::vector<unsigned int>> &expected_n_phases_per_composition)
       {
-        if (prm.get ("Viscosity averaging scheme") == "isostress")
-          viscosity_averaging_scheme = isostress;
-        else if (prm.get ("Viscosity averaging scheme") == "isostrain")
-          viscosity_averaging_scheme = isostrain;
-        else
-          {
-            AssertThrow(false, ExcMessage("Not a valid viscosity averaging scheme. Choose either isostress or isostrain."));
-          }
+        viscosity_averaging_scheme = ViscosityAveraging::parse("Viscosity averaging scheme",
+                                                               prm);
 
         // A background field is required by the subordinate material models
         number_of_chemical_compositions = this->introspection().n_chemical_composition_fields() + 1;
diff --git a/tests/composite_viscous_outputs.cc b/tests/composite_viscous_outputs.cc
index 0cec200924d..7b5efbcc0b7 100644
--- a/tests/composite_viscous_outputs.cc
+++ b/tests/composite_viscous_outputs.cc
@@ -20,6 +20,7 @@
 
 #include <aspect/simulator.h>
 #include <aspect/material_model/rheology/composite_visco_plastic.h>
+#include <aspect/simulator_signals.h>
 
 template <int dim>
 void f(const aspect::SimulatorAccess<dim> &simulator_access,
@@ -54,32 +55,32 @@ void f(const aspect::SimulatorAccess<dim> &simulator_access,
   prm.set("Include Drucker Prager plasticity in composite rheology", "true");
   prm.set("Peierls creep flow law", "viscosity approximation");
   prm.set("Maximum yield stress", "5e8");
-  composite_creep->parse_parameters(prm, n_phases);
+  composite_creep->parse_parameters(prm);
 
   std::unique_ptr<Rheology::DiffusionCreep<dim>> diffusion_creep;
   diffusion_creep = std::make_unique<Rheology::DiffusionCreep<dim>>();
   diffusion_creep->initialize_simulator (simulator_access.get_simulator());
   diffusion_creep->declare_parameters(prm);
-  diffusion_creep->parse_parameters(prm, n_phases);
+  diffusion_creep->parse_parameters(prm);
 
   std::unique_ptr<Rheology::DislocationCreep<dim>> dislocation_creep;
   dislocation_creep = std::make_unique<Rheology::DislocationCreep<dim>>();
   dislocation_creep->initialize_simulator (simulator_access.get_simulator());
   dislocation_creep->declare_parameters(prm);
-  dislocation_creep->parse_parameters(prm, n_phases);
+  dislocation_creep->parse_parameters(prm);
 
   std::unique_ptr<Rheology::PeierlsCreep<dim>> peierls_creep;
   peierls_creep = std::make_unique<Rheology::PeierlsCreep<dim>>();
   peierls_creep->initialize_simulator (simulator_access.get_simulator());
   peierls_creep->declare_parameters(prm);
-  peierls_creep->parse_parameters(prm, n_phases);
+  peierls_creep->parse_parameters(prm);
 
   std::unique_ptr<Rheology::DruckerPragerPower<dim>> drucker_prager_power;
   drucker_prager_power = std::make_unique<Rheology::DruckerPragerPower<dim>>();
   drucker_prager_power->initialize_simulator (simulator_access.get_simulator());
   drucker_prager_power->declare_parameters(prm);
   prm.set("Maximum yield stress", "5e8");
-  drucker_prager_power->parse_parameters(prm, n_phases);
+  drucker_prager_power->parse_parameters(prm);
   Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition);
 
   // The creep components are arranged in series with each other.
@@ -125,7 +126,7 @@ void f(const aspect::SimulatorAccess<dim> &simulator_access,
       temperature = 1000. + i*100.;
 
       // Compute the viscosity
-      viscosity = composite_creep->compute_composition_viscosity(pressure, temperature, grain_size, composition, strain_rate, partial_strain_rates);
+      viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates);
       total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.);
 
       // The creep strain rate is calculated by subtracting the strain rate
diff --git a/tests/composite_viscous_outputs_isostress.cc b/tests/composite_viscous_outputs_isostress.cc
index a1dd2ad1697..fcd433f9247 100644
--- a/tests/composite_viscous_outputs_isostress.cc
+++ b/tests/composite_viscous_outputs_isostress.cc
@@ -128,7 +128,7 @@ void f(const aspect::SimulatorAccess<dim> &simulator_access,
       temperature = 1000. + i*100.;
 
       // Compute the viscosity
-      viscosity = composite_creep->compute_isostress_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates);
+      viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates);
       total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.);
 
       // The creep strain rate is calculated by subtracting the strain rate
diff --git a/tests/composite_viscous_outputs_no_peierls.cc b/tests/composite_viscous_outputs_no_peierls.cc
index 6d3078eabb0..4b15fb65b54 100644
--- a/tests/composite_viscous_outputs_no_peierls.cc
+++ b/tests/composite_viscous_outputs_no_peierls.cc
@@ -34,6 +34,7 @@ void f(const aspect::SimulatorAccess<dim> &simulator_access,
 
   // First, we set up a few objects which are used by the rheology model.
   aspect::ParameterHandler prm;
+
   const std::vector<std::string> list_of_composition_names = simulator_access.introspection().get_composition_names();
   auto n_phases = std::make_unique<std::vector<unsigned int>>(1); // 1 phase per composition
   const unsigned int composition = 0;
@@ -54,26 +55,26 @@ void f(const aspect::SimulatorAccess<dim> &simulator_access,
   prm.set("Include Drucker Prager plasticity in composite rheology", "true");
   prm.set("Peierls creep flow law", "viscosity approximation");
   prm.set("Maximum yield stress", "5e8");
-  composite_creep->parse_parameters(prm, n_phases);
+  composite_creep->parse_parameters(prm);
 
   std::unique_ptr<Rheology::DiffusionCreep<dim>> diffusion_creep;
   diffusion_creep = std::make_unique<Rheology::DiffusionCreep<dim>>();
   diffusion_creep->initialize_simulator (simulator_access.get_simulator());
   diffusion_creep->declare_parameters(prm);
-  diffusion_creep->parse_parameters(prm, n_phases);
+  diffusion_creep->parse_parameters(prm);
 
   std::unique_ptr<Rheology::DislocationCreep<dim>> dislocation_creep;
   dislocation_creep = std::make_unique<Rheology::DislocationCreep<dim>>();
   dislocation_creep->initialize_simulator (simulator_access.get_simulator());
   dislocation_creep->declare_parameters(prm);
-  dislocation_creep->parse_parameters(prm, n_phases);
+  dislocation_creep->parse_parameters(prm);
 
   std::unique_ptr<Rheology::DruckerPragerPower<dim>> drucker_prager_power;
   drucker_prager_power = std::make_unique<Rheology::DruckerPragerPower<dim>>();
   drucker_prager_power->initialize_simulator (simulator_access.get_simulator());
   drucker_prager_power->declare_parameters(prm);
   prm.set("Maximum yield stress", "5e8");
-  drucker_prager_power->parse_parameters(prm, n_phases);
+  drucker_prager_power->parse_parameters(prm);
   Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition);
 
   // The creep components are arranged in series with each other.
@@ -118,7 +119,7 @@ void f(const aspect::SimulatorAccess<dim> &simulator_access,
       temperature = 1000. + i*100.;
 
       // Compute the viscosity
-      viscosity = composite_creep->compute_composition_viscosity(pressure, temperature, grain_size, composition, strain_rate, partial_strain_rates);
+      viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates);
       total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.);
 
       // The creep strain rate is calculated by subtracting the strain rate