Improve anisotropic viscosity code

Author: gassmoeller

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>

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

include/aspect/material_model/additional_outputs/anisotropic_viscosity.h

@@ -0,0 +1,61 @@
+/*
+  Copyright (C) 2025 - by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file doc/COPYING.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _aspect_material_model_additional_outputs_anisotropic_viscosity_h
+#define _aspect_material_model_additional_outputs_anisotropic_viscosity_h
+
+
+#include <aspect/material_model/interface.h>
+
+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;
+    };
+  }
+}
+
+#endif

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>

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>,

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] *

source/material_model/additional_outputs/anisotropic_viscosity.cc

@@ -0,0 +1,80 @@
+/*
+  Copyright (C) 2025 - by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file doc/COPYING.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#include <aspect/material_model/additional_outputs/anisotropic_viscosity.h>
+
+namespace aspect
+{
+  namespace MaterialModel
+  {
+    namespace
+    {
+      template <int dim>
+      std::vector<std::string> make_anisotropic_viscosity_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_anisotropic_viscosity_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;
+    }
+  }
+}
+
+// explicit instantiations
+namespace aspect
+{
+  namespace MaterialModel
+  {
+#define INSTANTIATE(dim) \
+  template class AnisotropicViscosity<dim>;
+
+    ASPECT_INSTANTIATE(INSTANTIATE)
+
+#undef INSTANTIATE
+  }
+}