Add viscosity additional outputs

doc/modules/changes/20260512_lhy11009

@@ -0,0 +1,3 @@
+Added: additional named outputs for diffusion and dislocation viscosities in the visco-plastic material model, including rheology-specific handling of inactive flow laws and visualization support through the additional material output system.
+<br>
+(Haoyuan Li, 2026/05/12)

include/aspect/material_model/rheology/visco_plastic.h

@@ -84,6 +84,35 @@ namespace aspect
 
     };
 
+    /**
+    * Additional output fields for diffusion and dislocation viscosities.
+    */
+    template <int dim>
+    class ViscosityAdditionalOutputs : public NamedAdditionalMaterialOutputs<dim>
+    {
+      public:
+        /**
+         * Constructor.
+         */
+        ViscosityAdditionalOutputs(const unsigned int n_points);
+
+        /**
+         * Return the nth output field.
+         */
+        std::vector<double>
+        get_nth_output(const unsigned int idx) const override;
+
+        /**
+         * Diffusion viscosities.
+         */
+        std::vector<double> diffusion_viscosities;
+
+        /**
+         * Dislocation viscosities.
+         */
+        std::vector<double> dislocation_viscosities;
+    };
+
     /**
      * A data structure with the output of calculate_isostrain_viscosities.
      */
@@ -117,6 +146,16 @@ namespace aspect
        * MaterialModel::PrescribedDilation::dilation_rhs_term.
        */
       std::vector<double> dilation_rhs_terms;
+
+      /**
+      * Diffusion viscosities for each composition.
+      */
+      std::vector<double> diffusion_viscosities;
+
+      /**
+       * Dislocation viscosities for each composition.
+       */
+      std::vector<double> dislocation_viscosities;
     };
 
     namespace Rheology
@@ -209,6 +248,20 @@ namespace aspect
                                     const MaterialModel::MaterialModelInputs<dim> &in,
                                     MaterialModel::MaterialModelOutputs<dim> &out,
                                     const IsostrainViscosities &isostrain_viscosities) const;
+          /**
+          * Create additional outputs for diffusion and dislocation viscosities.
+          */
+          void
+          create_viscosity_outputs(MaterialModel::MaterialModelOutputs<dim> &out) const;
+
+
+          /**
+           * Fill additional outputs for diffusion and dislocation viscosities.
+           */
+          void fill_viscosity_outputs(const unsigned int point_index,
+                                      const std::vector<double> &volume_fractions,
+                                      MaterialModel::MaterialModelOutputs<dim> &out,
+                                      const IsostrainViscosities &isostrain_viscosities) const;
 
           /**
            * Minimum strain rate used to stabilize the strain rate dependent rheology.

source/material_model/rheology/visco_plastic.cc

@@ -84,6 +84,42 @@ namespace aspect
 
 
 
+    namespace
+    {
+      std::vector<std::string> make_viscosity_additional_outputs_names()
+      {
+        std::vector<std::string> names;
+        names.emplace_back("diffusion_viscosity");
+        names.emplace_back("dislocation_viscosity");
+        return names;
+      }
+    }
+
+
+
+    template <int dim>
+    ViscosityAdditionalOutputs<dim>::ViscosityAdditionalOutputs(const unsigned int n_points)
+      : NamedAdditionalMaterialOutputs<dim>(make_viscosity_additional_outputs_names()),
+        diffusion_viscosities(n_points, numbers::signaling_nan<double>()),
+        dislocation_viscosities(n_points, numbers::signaling_nan<double>())
+    {}
+
+
+
+    template <int dim>
+    std::vector<double>
+    ViscosityAdditionalOutputs<dim>::get_nth_output(const unsigned int idx) const
+    {
+      AssertIndexRange(idx, 2);
+
+      if (idx == 0)
+        return diffusion_viscosities;
+      else
+        return dislocation_viscosities;
+    }
+
+
+
     namespace Rheology
     {
 
@@ -110,6 +146,8 @@ namespace aspect
         output_parameters.drucker_prager_parameters.resize(volume_fractions.size());
         output_parameters.dilation_lhs_terms.resize(volume_fractions.size(), numbers::signaling_nan<double>());
         output_parameters.dilation_rhs_terms.resize(volume_fractions.size(), numbers::signaling_nan<double>());
+        output_parameters.diffusion_viscosities.resize(volume_fractions.size(), numbers::signaling_nan<double>());
+        output_parameters.dislocation_viscosities.resize(volume_fractions.size(), numbers::signaling_nan<double>());
 
         // Assemble current and old stress tensor if elastic behavior is enabled
         SymmetricTensor<2, dim> stress_0_advected = numbers::signaling_nan<SymmetricTensor<2, dim>>();
@@ -220,12 +258,14 @@ namespace aspect
                   {
                     non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_diffusion, j, i,
                                                                                                   CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::diffusion);
+                    output_parameters.diffusion_viscosities[j] = non_yielding_viscosity;
                     break;
                   }
                   case dislocation:
                   {
                     non_yielding_viscosity = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i,
                                                                                                   CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
+                    output_parameters.dislocation_viscosities[j] = non_yielding_viscosity;
                     break;
                   }
                   case frank_kamenetskii:
@@ -244,6 +284,9 @@ namespace aspect
                                                                 CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
                     non_yielding_viscosity = (scaled_viscosity_diffusion * scaled_viscosity_dislocation)/
                                              (scaled_viscosity_diffusion + scaled_viscosity_dislocation);
+
+                    output_parameters.diffusion_viscosities[j] = scaled_viscosity_diffusion;
+                    output_parameters.dislocation_viscosities[j] = scaled_viscosity_dislocation;
                     break;
                   }
                   case minimum_diffusion_dislocation:
@@ -253,6 +296,9 @@ namespace aspect
                     const double scaled_viscosity_dislocation = compositional_viscosity_prefactors.compute_viscosity(in, viscosity_dislocation, j, i,
                                                                 CompositionalViscosityPrefactors<dim>::ModifiedFlowLaws::dislocation);
                     non_yielding_viscosity = std::min(scaled_viscosity_diffusion, scaled_viscosity_dislocation);
+
+                    output_parameters.diffusion_viscosities[j] = scaled_viscosity_diffusion;
+                    output_parameters.dislocation_viscosities[j] = scaled_viscosity_dislocation;
                     break;
                   }
                   default:
@@ -984,6 +1030,95 @@ namespace aspect
           }
       }
 
+      template <int dim>
+      void
+      Rheology::ViscoPlastic<dim>::create_viscosity_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
+      {
+        if (out.template has_additional_output_object<ViscosityAdditionalOutputs<dim>>() == false)
+          {
+            const unsigned int n_points = out.n_evaluation_points();
+
+            out.additional_outputs.push_back(
+              std::make_unique<ViscosityAdditionalOutputs<dim>>(n_points));
+          }
+      }
+
+      template <int dim>
+      void
+      Rheology::ViscoPlastic<dim>::fill_viscosity_outputs(
+        const unsigned int i,
+        const std::vector<double> &volume_fractions,
+        MaterialModel::MaterialModelOutputs<dim> &out,
+        const IsostrainViscosities &isostrain_viscosities) const
+      {
+        if (const std::shared_ptr<ViscosityAdditionalOutputs<dim>>
+            viscosity_out = out.template get_additional_output_object<ViscosityAdditionalOutputs<dim>>())
+          {
+            switch (viscous_flow_law)
+              {
+                case diffusion:
+                {
+                  viscosity_out->diffusion_viscosities[i]
+                    = MaterialUtilities::average_value(
+                        volume_fractions,
+                        isostrain_viscosities.diffusion_viscosities,
+                        viscosity_averaging);
+
+                  viscosity_out->dislocation_viscosities[i]
+                    = std::numeric_limits<double>::max();
+
+                  break;
+                }
+
+                case dislocation:
+                {
+                  viscosity_out->diffusion_viscosities[i]
+                    = std::numeric_limits<double>::max();
+
+                  viscosity_out->dislocation_viscosities[i]
+                    = MaterialUtilities::average_value(
+                        volume_fractions,
+                        isostrain_viscosities.dislocation_viscosities,
+                        viscosity_averaging);
+
+                  break;
+                }
+
+                case composite:
+                case minimum_diffusion_dislocation:
+                {
+                  viscosity_out->diffusion_viscosities[i]
+                    = MaterialUtilities::average_value(
+                        volume_fractions,
+                        isostrain_viscosities.diffusion_viscosities,
+                        viscosity_averaging);
+
+                  viscosity_out->dislocation_viscosities[i]
+                    = MaterialUtilities::average_value(
+                        volume_fractions,
+                        isostrain_viscosities.dislocation_viscosities,
+                        viscosity_averaging);
+
+                  break;
+                }
+
+                case frank_kamenetskii:
+                {
+                  viscosity_out->diffusion_viscosities[i]
+                    = std::numeric_limits<double>::max();
+
+                  viscosity_out->dislocation_viscosities[i]
+                    = std::numeric_limits<double>::max();
+
+                  break;
+                }
+
+                default:
+                  AssertThrow(false, ExcInternalError());
+              }
+          }
+      }
+
     }
   }
 }
@@ -995,6 +1130,7 @@ namespace aspect
   {
 #define INSTANTIATE(dim) \
   template class PlasticAdditionalOutputs<dim>; \
+  template class ViscosityAdditionalOutputs<dim>; \
   \
   namespace Rheology \
   { \

source/material_model/visco_plastic.cc

@@ -251,6 +251,8 @@ namespace aspect
               isostrain_viscosities.composition_yielding.clear();
               isostrain_viscosities.composition_viscosities.clear();
               isostrain_viscosities.drucker_prager_parameters.clear();
+              isostrain_viscosities.diffusion_viscosities.clear();
+              isostrain_viscosities.dislocation_viscosities.clear();
 
               out.viscosities[i] = numbers::signaling_nan<double>();
 
@@ -270,11 +272,14 @@ namespace aspect
           // TODO only when requests_property is set to reaction_terms
           rheology->strain_rheology.fill_reaction_outputs(in, i, rheology->min_strain_rate, plastic_yielding, out);
 
-          // Fill plastic outputs if they exist.
+          // Fill plastic outputs and additional viscosity outputs if they exist.
           // The values in isostrain_viscosities only make sense when the calculate_isostrain_viscosities function
           // has been called.
           if (in.requests_property(MaterialProperties::additional_outputs))
-            rheology->fill_plastic_outputs(i, volume_fractions, plastic_yielding, in, out, isostrain_viscosities);
+            {
+              rheology->fill_plastic_outputs(i, volume_fractions, plastic_yielding, in, out, isostrain_viscosities);
+              rheology->fill_viscosity_outputs(i, volume_fractions, out, isostrain_viscosities);
+            }
 
           if (this->get_parameters().enable_elasticity)
             {
@@ -476,6 +481,7 @@ namespace aspect
     ViscoPlastic<dim>::create_additional_named_outputs (MaterialModel::MaterialModelOutputs<dim> &out) const
     {
       rheology->create_plastic_outputs(out);
+      rheology->create_viscosity_outputs(out);
 
       if (this->get_parameters().enable_elasticity)
         rheology->elastic_rheology.create_elastic_additional_outputs(out);