Unify anisotropic viscosity Stokes assemblers / Fix anisotropic viscosity cookbook

cookbooks/anisotropic_viscosity/AV_Rayleigh_Taylor.prm

@@ -62,7 +62,7 @@ subsection Initial composition model
 
   subsection Function
     set Variable names      = x,z
-    set Function constants  = pi=3.1415926;
+    set Function constants  = pi=3.1415926
     set Function expression = 0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02)); \
                             (0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02))) > 0.8 ? sin(45*pi/180) : 0.0; \
                             (0.5*(1+tanh((z-0.85-0.02*cos(2*pi*x/2))/0.02))) > 0.8 ? cos(45*pi/180) : 0.0;

include/aspect/heating_model/shear_heating_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 LICENSE.  If not see
+ <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef _aspect_heating_model_shear_heating_anisotropic_viscosity_h
+#define _aspect_heating_model_shear_heating_anisotropic_viscosity_h
+
+#include <aspect/heating_model/interface.h>
+
+#include <aspect/simulator_access.h>
+
+namespace aspect
+{
+  namespace HeatingModel
+  {
+    /**
+     * 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
+     * tensor to compute the stress that is used when computing the shear heating.
+     *
+     * @ingroup HeatingModels
+     */
+    template <int dim>
+    class ShearHeatingAnisotropicViscosity : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
+    {
+      public:
+        /**
+         * Compute the heating model outputs for this class.
+         */
+        void
+        evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
+                  const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
+                  HeatingModel::HeatingModelOutputs &heating_model_outputs) const override;
+
+        /**
+         * Allow the heating model to attach additional material model outputs.
+         */
+        void
+        create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &material_model_outputs) const override;
+    };
+  }
+}
+
+#endif

include/aspect/simulator/assemblers/stokes_anisotropic_viscosity.h

@@ -0,0 +1,172 @@
+/*
+  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_simulator_assemblers_stokes_anisotropic_viscosity_h
+#define _aspect_simulator_assemblers_stokes_anisotropic_viscosity_h
+
+
+#include <aspect/simulator/assemblers/interface.h>
+#include <aspect/simulator_access.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;
+    };
+
+    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.
+    */
+    template <int dim>
+    class StokesPreconditionerAnisotropicViscosity : public Assemblers::Interface<dim>,
+      public SimulatorAccess<dim>
+    {
+      public:
+        void
+        execute(internal::Assembly::Scratch::ScratchBase<dim>   &scratch,
+                internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
+
+        /**
+         * Create AnisotropicViscosities.
+         */
+        void
+        create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const override;
+    };
+
+    /**
+     * A class containing the functions to assemble the compressible adjustment
+     * to the Stokes preconditioner.
+     */
+    template <int dim>
+    class StokesCompressiblePreconditionerAnisotropicViscosity : public Assemblers::Interface<dim>,
+      public SimulatorAccess<dim>
+    {
+      public:
+        void
+        execute(internal::Assembly::Scratch::ScratchBase<dim>   &scratch,
+                internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
+    };
+
+    /**
+     * This class assembles the terms for the matrix and right-hand-side of the incompressible
+     * Stokes equation for the current cell.
+     */
+    template <int dim>
+    class StokesIncompressibleTermsAnisotropicViscosity : public Assemblers::Interface<dim>,
+      public SimulatorAccess<dim>
+    {
+      public:
+        void
+        execute(internal::Assembly::Scratch::ScratchBase<dim>   &scratch,
+                internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
+
+        /**
+         * Create AdditionalMaterialOutputsStokesRHS if we need to do so.
+         */
+        void
+        create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const override;
+    };
+
+    /**
+     * 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.
+     */
+    template <int dim>
+    class StokesCompressibleStrainRateViscosityTermAnisotropicViscosity : public Assemblers::Interface<dim>,
+      public SimulatorAccess<dim>
+    {
+      public:
+        void
+        execute(internal::Assembly::Scratch::ScratchBase<dim>   &scratch,
+                internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
+    };
+  }
+}
+
+
+#endif

source/heating_model/shear_heating_anisotropic_viscosity.cc

@@ -0,0 +1,123 @@
+/*
+ Copyright (C) 2015 - 2023 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 LICENSE.  If not see
+ <http://www.gnu.org/licenses/>.
+ */
+
+#include <aspect/heating_model/shear_heating_anisotropic_viscosity.h>
+
+#include <aspect/material_model/interface.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>
+
+namespace aspect
+{
+  namespace HeatingModel
+  {
+    template <int dim>
+    void
+    ShearHeatingAnisotropicViscosity<dim>::
+    evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
+              const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
+              HeatingModel::HeatingModelOutputs &heating_model_outputs) const
+    {
+      Assert(heating_model_outputs.heating_source_terms.size() == material_model_inputs.position.size(),
+             ExcMessage ("Heating outputs need to have the same number of entries as the material model inputs."));
+
+      // Some material models provide dislocation viscosities and boundary area work fractions
+      // as additional material outputs. If they are attached, use them.
+      const std::shared_ptr<const ShearHeatingOutputs<dim>> shear_heating_out =
+        material_model_outputs.template get_additional_output_object<ShearHeatingOutputs<dim>>();
+
+      const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity =
+        material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
+
+      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> stress =
+            2 * material_model_outputs.viscosities[q] *
+            (this->get_material_model().is_compressible()
+             ?
+             directed_strain_rate - 1./3. * trace(directed_strain_rate) * unit_symmetric_tensor<dim>()
+             :
+             directed_strain_rate);
+
+          const SymmetricTensor<2,dim> deviatoric_strain_rate =
+            (this->get_material_model().is_compressible()
+             ?
+             material_model_inputs.strain_rate[q]
+             - 1./3. * trace(material_model_inputs.strain_rate[q]) * unit_symmetric_tensor<dim>()
+             :
+             material_model_inputs.strain_rate[q]);
+
+          heating_model_outputs.heating_source_terms[q] = stress * deviatoric_strain_rate;
+
+          // If shear heating work fractions are provided, reduce the
+          // overall heating by this amount (which is assumed to be converted into other forms of energy)
+          if (shear_heating_out != nullptr)
+            heating_model_outputs.heating_source_terms[q] *= shear_heating_out->shear_heating_work_fractions[q];
+
+          heating_model_outputs.lhs_latent_heat_terms[q] = 0.0;
+        }
+    }
+
+
+
+    template <int dim>
+    void
+    ShearHeatingAnisotropicViscosity<dim>::
+    create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &material_model_outputs) const
+    {
+      const unsigned int n_points = material_model_outputs.viscosities.size();
+
+      if (material_model_outputs.template has_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>() == false)
+        {
+          material_model_outputs.additional_outputs.push_back(
+            std::make_unique<MaterialModel::AnisotropicViscosity<dim>> (n_points));
+        }
+
+      this->get_material_model().create_additional_named_outputs(material_model_outputs);
+    }
+  }
+}
+
+
+
+// explicit instantiations
+namespace aspect
+{
+  namespace HeatingModel
+  {
+    ASPECT_REGISTER_HEATING_MODEL(ShearHeatingAnisotropicViscosity,
+                                  "anisotropic shear heating",
+                                  "Implementation of 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 "
+                                  "tensor to compute the stress that is used when computing the shear heating.")
+  }
+}