Improve anisotropic viscosity code

gassmoeller · gassmoeller · commit a266d9ca8a0f · 2025-11-11T15:31:57.000+01:00
diff --git a/cookbooks/anisotropic_viscosity/av_material.cc b/cookbooks/anisotropic_viscosity/av_material.cc
@@ -20,7 +20,7 @@
 
 #include <aspect/introspection.h>
 #include <aspect/material_model/interface.h>
-#include <aspect/plugins.h>
+#include <aspect/material_model/additional_outputs/anisotropic_viscosity.h>
 #include <aspect/simulator/assemblers/stokes_anisotropic_viscosity.h>
 #include <deal.II/base/parameter_handler.h>
 #include <deal.II/base/patterns.h>
diff --git a/include/aspect/heating_model/shear_heating_anisotropic_viscosity.h b/include/aspect/heating_model/shear_heating_anisotropic_viscosity.h
@@ -31,8 +31,8 @@ namespace aspect
   {
     /**
      * A class that implements a standard model for shear heating extended for an
-     * anisotropic viscosity tensor. If the material model provides a stress-
-     * strain director tensor, then the strain-rate is multiplied with this
+     * anisotropic viscosity tensor. If the material model provides a stress-strain
+     * director tensor, then the strain-rate is multiplied with this
      * tensor to compute the stress that is used when computing the shear heating.
      *
      * @ingroup HeatingModels
diff --git a/include/aspect/simulator/assemblers/stokes.h b/include/aspect/simulator/assemblers/stokes.h
@@ -78,7 +78,7 @@ namespace aspect
     };
 
     /**
-     * This class assembles the term that arises in the viscosity term of Stokes matrix for
+     * This class assembles the term that arises in the viscosity term of the Stokes matrix for
      * compressible models, because the divergence of the velocity is not longer zero.
      */
     template <int dim>
diff --git a/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h b/include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h
@@ -27,80 +27,11 @@
 
 namespace aspect
 {
-  namespace MaterialModel
-  {
-    /**
-      * Additional output fields for anisotropic viscosities to be added to
-      * the MaterialModel::MaterialModelOutputs structure and filled in the
-      * MaterialModel::Interface::evaluate() function.
-      */
-    template <int dim>
-    class AnisotropicViscosity : public NamedAdditionalMaterialOutputs<dim>
-    {
-      public:
-        AnisotropicViscosity(const unsigned int n_points);
-
-        std::vector<double>
-        get_nth_output(const unsigned int idx) const override;
-
-        /**
-         * Stress-strain "director" tensors at the given positions. This
-         * variable is used to implement anisotropic viscosity.
-         *
-         * @note The strain rate term in equation (1) of the manual will be
-         * multiplied by this tensor *and* the viscosity scalar ($\eta$), as
-         * described in the manual section titled "Constitutive laws". This
-         * variable is assigned the rank-four identity tensor by default.
-         * This leaves the isotropic constitutive law unchanged if the material
-         * model does not explicitly assign a value.
-         */
-        std::vector<SymmetricTensor<4,dim>> stress_strain_directors;
-    };
-
-    namespace
-    {
-      template <int dim>
-      std::vector<std::string> make_AnisotropicViscosity_additional_outputs_names()
-      {
-        std::vector<std::string> names;
-
-        for (unsigned int i = 0; i < Tensor<4,dim>::n_independent_components ; ++i)
-          {
-            TableIndices<4> indices(Tensor<4,dim>::unrolled_to_component_indices(i));
-            names.push_back("anisotropic_viscosity"+std::to_string(indices[0])+std::to_string(indices[1])+std::to_string(indices[2])+std::to_string(indices[3]));
-          }
-        return names;
-      }
-    }
-
-
-
-    template <int dim>
-    AnisotropicViscosity<dim>::AnisotropicViscosity (const unsigned int n_points)
-      :
-      NamedAdditionalMaterialOutputs<dim>(make_AnisotropicViscosity_additional_outputs_names<dim>()),
-      stress_strain_directors(n_points, dealii::identity_tensor<dim> ())
-    {}
-
-
-
-    template <int dim>
-    std::vector<double>
-    AnisotropicViscosity<dim>::get_nth_output(const unsigned int idx) const
-    {
-      std::vector<double> output(stress_strain_directors.size());
-      for (unsigned int i = 0; i < stress_strain_directors.size() ; ++i)
-        {
-          output[i]= stress_strain_directors[i][Tensor<4,dim>::unrolled_to_component_indices(idx)];
-        }
-      return output;
-    }
-  }
-
   namespace Assemblers
   {
     /**
-    * A class containing the functions to assemble the Stokes preconditioner.
+    * A class containing the functions to assemble the Stokes preconditioner for the
+    * case of anisotropic viscosities.
     */
     template <int dim>
     class StokesPreconditionerAnisotropicViscosity : public Assemblers::Interface<dim>,
@@ -120,7 +51,7 @@ namespace aspect
 
     /**
      * A class containing the functions to assemble the compressible adjustment
-     * to the Stokes preconditioner.
+     * to the Stokes preconditioner for the case of anisotropic viscosities.
      */
     template <int dim>
     class StokesCompressiblePreconditionerAnisotropicViscosity : public Assemblers::Interface<dim>,
@@ -134,7 +65,7 @@ namespace aspect
 
     /**
      * This class assembles the terms for the matrix and right-hand-side of the incompressible
-     * Stokes equation for the current cell.
+     * Stokes equation for the case of anisotropic viscosities.
      */
     template <int dim>
     class StokesIncompressibleTermsAnisotropicViscosity : public Assemblers::Interface<dim>,
@@ -153,8 +84,9 @@ namespace aspect
     };
 
     /**
-     * This class assembles the term that arises in the viscosity term of Stokes matrix for
-     * compressible models, because the divergence of the velocity is not longer zero.
+     * This class assembles the term that arises in the viscosity term of the Stokes matrix for
+     * compressible models for the case of anisotropic viscosities, because the divergence
+     * of the velocity is not longer zero.
      */
     template <int dim>
     class StokesCompressibleStrainRateViscosityTermAnisotropicViscosity : public Assemblers::Interface<dim>,
diff --git a/source/heating_model/shear_heating_anisotropic_viscosity.cc b/source/heating_model/shear_heating_anisotropic_viscosity.cc
@@ -21,8 +21,8 @@
 #include <aspect/heating_model/shear_heating_anisotropic_viscosity.h>
 
 #include <aspect/material_model/interface.h>
+#include <aspect/material_model/additional_outputs/anisotropic_viscosity.h>
 #include <aspect/heating_model/shear_heating.h>
-#include <aspect/simulator/assemblers/stokes_anisotropic_viscosity.h>
 
 #include <deal.II/base/symmetric_tensor.h>
 #include <deal.II/base/signaling_nan.h>
@@ -49,15 +49,15 @@ namespace aspect
       const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity =
         material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
 
+      Assert(anisotropic_viscosity != nullptr,
+             ExcMessage("This heating plugin should only be used with material models that provide "
+                        "an anisotropic viscosity tensor, but none was provided."));
+
       for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
         {
           // If there is an anisotropic viscosity, use it to compute the correct stress
-          const SymmetricTensor<2,dim> &directed_strain_rate = ((anisotropic_viscosity != nullptr)
-                                                                ?
-                                                                anisotropic_viscosity->stress_strain_directors[q]
-                                                                * material_model_inputs.strain_rate[q]
-                                                                :
-                                                                material_model_inputs.strain_rate[q]);
+          const SymmetricTensor<2,dim> &directed_strain_rate = anisotropic_viscosity->stress_strain_directors[q]
+                                                               * material_model_inputs.strain_rate[q];
 
           const SymmetricTensor<2,dim> stress =
             2 * material_model_outputs.viscosities[q] *
diff --git a/source/simulator/assemblers/stokes_anisotropic_viscosity.cc b/source/simulator/assemblers/stokes_anisotropic_viscosity.cc
@@ -21,6 +21,7 @@
 #include <aspect/simulator/assemblers/stokes_anisotropic_viscosity.h>
 
 #include <aspect/gravity_model/interface.h>
+#include <aspect/material_model/additional_outputs/anisotropic_viscosity.h>
 
 #include <deal.II/base/signaling_nan.h>
 
@@ -37,9 +38,13 @@ namespace aspect
       internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesPreconditioner<dim>&> (scratch_base);
       internal::Assembly::CopyData::StokesPreconditioner<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesPreconditioner<dim>&> (data_base);
 
-      std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity =
+      const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity =
         scratch.material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
 
+      Assert(anisotropic_viscosity != nullptr,
+             ExcMessage("This assembler should only be used with material models that provide "
+                        "an anisotropic viscosity tensor, but none was provided."));
+
       const Introspection<dim> &introspection = this->introspection();
       const FiniteElement<dim> &fe = this->get_fe();
       const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
@@ -78,30 +83,16 @@ namespace aspect
 
           const double eta = scratch.material_model_outputs.viscosities[q];
           const double one_over_eta = 1. / eta;
-
-          const bool use_tensor = (anisotropic_viscosity != nullptr);
-
-          const SymmetricTensor<4, dim> &stress_strain_director = (use_tensor)
-                                                                  ?
-                                                                  anisotropic_viscosity->stress_strain_directors[q]
-                                                                  :
-                                                                  dealii::identity_tensor<dim>();
-
-
-
+          const SymmetricTensor<4, dim> &stress_strain_director = anisotropic_viscosity->stress_strain_directors[q];
           const double JxW = scratch.finite_element_values.JxW(q);
 
           for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i)
             for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j)
               if (scratch.dof_component_indices[i] ==
                   scratch.dof_component_indices[j])
-                data.local_matrix(i, j) += ((
-                                              use_tensor ?
-                                              2.0 * eta * (scratch.grads_phi_u[i]
-                                                           * stress_strain_director
-                                                           * scratch.grads_phi_u[j]) :
-                                              2.0 * eta * (scratch.grads_phi_u[i]
-                                                           * scratch.grads_phi_u[j]))
+                data.local_matrix(i, j) += (2.0 * eta * (scratch.grads_phi_u[i]
+                                                         * stress_strain_director
+                                                         * scratch.grads_phi_u[j])
                                             + one_over_eta * pressure_scaling
                                             * pressure_scaling
                                             * (scratch.phi_p[i]
@@ -140,6 +131,10 @@ namespace aspect
       const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity =
         scratch.material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
 
+      Assert(anisotropic_viscosity != nullptr,
+             ExcMessage("This assembler should only be used with material models that provide "
+                        "an anisotropic viscosity tensor, but none was provided."));
+
       const Introspection<dim> &introspection = this->introspection();
       const FiniteElement<dim> &fe = this->get_fe();
       const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
@@ -174,26 +169,15 @@ namespace aspect
             }
 
           const double eta_two_thirds = scratch.material_model_outputs.viscosities[q] * 2.0 / 3.0;
-
-          const bool use_tensor = (anisotropic_viscosity != nullptr);
-
-          const SymmetricTensor<4, dim> &stress_strain_director = (use_tensor)
-                                                                  ?
-                                                                  anisotropic_viscosity->stress_strain_directors[q]
-                                                                  :
-                                                                  dealii::identity_tensor<dim>();
-
+          const SymmetricTensor<4, dim> &stress_strain_director = anisotropic_viscosity->stress_strain_directors[q];
           const double JxW = scratch.finite_element_values.JxW(q);
 
           for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i)
             for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j)
               if (scratch.dof_component_indices[i] ==
                   scratch.dof_component_indices[j])
-                data.local_matrix(i, j) += (- (use_tensor ?
-                                               eta_two_thirds * (scratch.div_phi_u[i] * trace(stress_strain_director * scratch.grads_phi_u[j]))
-                                               :
-                                               eta_two_thirds * (scratch.div_phi_u[i] * scratch.div_phi_u[j])
-                                              ))
+                data.local_matrix(i, j) += (- eta_two_thirds *
+                                            (scratch.div_phi_u[i] * trace(stress_strain_director * scratch.grads_phi_u[j])))
                                            * JxW;
         }
     }
@@ -212,6 +196,10 @@ namespace aspect
       const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity =
         scratch.material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
 
+      Assert(anisotropic_viscosity != nullptr,
+             ExcMessage("This assembler should only be used with material models that provide "
+                        "an anisotropic viscosity tensor, but none was provided."));
+
       const Introspection<dim> &introspection = this->introspection();
       const FiniteElement<dim> &fe = this->get_fe();
       const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
@@ -247,13 +235,7 @@ namespace aspect
                               :
                               numbers::signaling_nan<double>());
 
-          const bool use_tensor = (anisotropic_viscosity != nullptr);
-
-          const SymmetricTensor<4, dim> &stress_strain_director = (use_tensor)
-                                                                  ?
-                                                                  anisotropic_viscosity->stress_strain_directors[q]
-                                                                  :
-                                                                  dealii::identity_tensor<dim>();
+          const SymmetricTensor<4, dim> &stress_strain_director = anisotropic_viscosity->stress_strain_directors[q];
 
           const Tensor<1,dim>
           gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q));
@@ -274,18 +256,15 @@ namespace aspect
               if (scratch.rebuild_stokes_matrix)
                 for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
                   {
-                    data.local_matrix(i,j) += ( (use_tensor ?
-                                                 eta * 2.0 * (scratch.grads_phi_u[i] * stress_strain_director * scratch.grads_phi_u[j])
-                                                 :
-                                                 eta * 2.0 * (scratch.grads_phi_u[i] * scratch.grads_phi_u[j]))
-                                                // assemble \nabla p as -(p, div v):
-                                                - (pressure_scaling *
-                                                   scratch.div_phi_u[i] * scratch.phi_p[j])
-                                                // assemble the term -div(u) as -(div u, q).
-                                                // Note the negative sign to make this
-                                                // operator adjoint to the grad p term:
-                                                - (pressure_scaling *
-                                                   scratch.phi_p[i] * scratch.div_phi_u[j]))
+                    data.local_matrix(i,j) += (eta * 2.0 * (scratch.grads_phi_u[i] * stress_strain_director * scratch.grads_phi_u[j])
+                                               // assemble \nabla p as -(p, div v):
+                                               - (pressure_scaling *
+                                                  scratch.div_phi_u[i] * scratch.phi_p[j])
+                                               // assemble the term -div(u) as -(div u, q).
+                                               // Note the negative sign to make this
+                                               // operator adjoint to the grad p term:
+                                               - (pressure_scaling *
+                                                  scratch.phi_p[i] * scratch.div_phi_u[j]))
                                               * JxW;
                   }
             }
@@ -336,6 +315,10 @@ namespace aspect
       const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity =
         scratch.material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
 
+      Assert(anisotropic_viscosity != nullptr,
+             ExcMessage("This assembler should only be used with material models that provide "
+                        "an anisotropic viscosity tensor, but none was provided."));
+
       const Introspection<dim> &introspection = this->introspection();
       const FiniteElement<dim> &fe = this->get_fe();
       const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
@@ -357,25 +340,14 @@ namespace aspect
 
           // Viscosity scalar
           const double eta_two_thirds = scratch.material_model_outputs.viscosities[q] * 2.0 / 3.0;
-
-          const bool use_tensor = (anisotropic_viscosity != nullptr);
-
-          const SymmetricTensor<4, dim> &stress_strain_director = (use_tensor)
-                                                                  ?
-                                                                  anisotropic_viscosity->stress_strain_directors[q]
-                                                                  :
-                                                                  dealii::identity_tensor<dim>();
-
+          const SymmetricTensor<4, dim> &stress_strain_director = anisotropic_viscosity->stress_strain_directors[q];
           const double JxW = scratch.finite_element_values.JxW(q);
 
           for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
             for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
               {
-                data.local_matrix(i,j) += (- (use_tensor ?
-                                              eta_two_thirds * (scratch.div_phi_u[i] * trace(stress_strain_director * scratch.grads_phi_u[j]))
-                                              :
-                                              eta_two_thirds * (scratch.div_phi_u[i] * scratch.div_phi_u[j])
-                                             ))
+                data.local_matrix(i,j) += (-eta_two_thirds * (scratch.div_phi_u[i] * trace(stress_strain_director * scratch.grads_phi_u[j]))
+                                          )
                                           * JxW;
               }
         }
diff --git a/tests/anisotropic_viscosity.cc b/tests/anisotropic_viscosity.cc
@@ -21,6 +21,7 @@
 #include <aspect/simulator_access.h>
 
 #include <aspect/material_model/simple.h>
+#include <aspect/material_model/additional_outputs/anisotropic_viscosity.h>
 #include <aspect/heating_model/shear_heating.h>
 #include <aspect/heating_model/interface.h>
 #include <aspect/gravity_model/interface.h>

Original file line number	Diff line number	Diff line change
`@@ -31,8 +31,8 @@ namespace aspect`
`31`	`31`	`{`
`32`	`32`	`/**`
`33`	`33`	`* A class that implements a standard model for shear heating extended for an`
`34`		`- * anisotropic viscosity tensor. If the material model provides a stress-`
`35`		`- * strain director tensor, then the strain-rate is multiplied with this`
	`34`	`+ * anisotropic viscosity tensor. If the material model provides a stress-strain`
	`35`	`+ * director tensor, then the strain-rate is multiplied with this`
`36`	`36`	`* tensor to compute the stress that is used when computing the shear heating.`
`37`	`37`	`*`
`38`	`38`	`* @ingroup HeatingModels`