CPO-induced anisotropic viscosity cookbook by Wang-yijun · Pull Request #6615 · geodynamics/aspect

Wang-yijun · 2025-07-17T15:58:26Z

This PR adds a cookbook demonstrating the use of CPO-induced anisotropic viscosity (AV), developed by me and @KiralyAgi. It includes:

The AV material model plugin and an anisotropic stress postprocessor plugin.
A simple shear box test model showcasing the AV material implementation.
A page in the documentation explaining the setup and requirements.

Let me know if any changes are needed!

gassmoeller

Very cool to have a cookbook for this feature. Sorry it took a while to get to this. I have left a first batch of comments, I havent yet looked at the code implementation in details, because I left some questions that would be useful to discuss beforehand.

Please in addition to my comments:

Add a test case in tests/ that runs the added cookbook reasonably fast. You can see how this is done in tests/magnetic_stripes.cc/prm.
Add a changelog entry in doc/modules

Let me know when someone can take another look.

gassmoeller · 2025-07-28T09:52:54Z

+
+## Introduction
+
+The directional-dependence (anisotropy) in the viscous properties of olivine crystals suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This CPO-induced anisotropic viscosity material model computes anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and modifies the subsequent deformation.


Please write out CPO at least once and explain what it means.

gassmoeller · 2025-07-28T09:53:23Z

+
+$$\dot \varepsilon_{ij} = \gamma J(\sigma_{ij})^{(n-1)} A_{ij} \sigma_{ij}$$
+
+where $\gamma$ is the part of fluidity which is temperature- and grain-size dependent:


also explain fluidity

gassmoeller · 2025-07-28T09:55:31Z

+
+$$\gamma=\gamma_0 exp \left(\frac{-Q}{RT} \right) /d^m$$
+
+$\gamma_0=1.1\times 10^{5}$ is the isotropic fluidity, $Q=530 kJ/mol$ is the activation energy, $R=8.314 m^3 \cdot Pa \cdot K^{−1} \cdot mol^{−1}$ is the gas constant, $d=0.001 m$ is the grain size, and $m=0.73$ is the grain size exponent. These values are taken from {cite:t}`hansen:etal:2016` and {cite:t}`HK04`. $J(\sigma_{ij})$ is the equivalent yield stress, where $\sigma_{ij}$ is the deviatoric (anisotropic) stress computed using the tensorial and scalar component of the anisotropic viscosity:


is suppose the values of gamma_0 and the other quantities is for some special material (olivine) this should be stated somewhere here. Also mention that these are based on rock experiments.

gassmoeller · 2025-07-28T09:57:22Z

+
+$$J(\sigma_{ij})$$ and $A_{ij}$ are computed using Hill coefficients $H, J, K, L, M,$ and $N$ {cite}`hill:1948`, which are obtained from regression analysis of a texture database constructed with olivine textures from laboratory experiments, shear box models, and subduction models (Kiraly et al., in rev.).
+
+In this material model plugin, strain rate, density, temperature, and other parameters are taken as input to compute the viscosity, which is passed into the Stokes system to compute the stress (Figure 1). As a result, we adapt eq. (1) to be:


If you need to reference equations you can do so using the commands explained here: https://aspect-documentation--6615.org.readthedocs.build/en/6615/user/extending/write-a-cookbook/quickref.html#id1. Eq(1) is not numbered, because you didnt use the {math} environment. And the same for the figure. Also check the other equation and figure references in this file.

Suggested change

In this material model plugin, strain rate, density, temperature, and other parameters are taken as input to compute the viscosity, which is passed into the Stokes system to compute the stress (Figure 1). As a result, we adapt eq. (1) to be:

In this material model plugin, strain rate, density, temperature, and other parameters are taken as input to compute the anisotropic viscosity, which is passed into the Stokes system to compute the stress (Figure 1). As a result, we adapt eq. (1) to be:

gassmoeller · 2025-07-28T10:04:34Z

+0 & 0 & 0 & L & 0 & 0 \\
+0 & 0 & 0 & 0 & M & 0 \\
+0 & 0 & 0 & 0 & 0 & N
+\end{matrix} \right]$$


This math environment doesnt properly close, see here: https://aspect-documentation--6615.org.readthedocs.build/en/6615/user/cookbooks/cookbooks/CPO_induced_anisotropic_viscosity/doc/CPO_induced_anisotropic_viscosity.html

You will likely need to use a full multiline math environment, see my comment below.

gassmoeller · 2025-07-28T11:59:44Z

+## References:
+- Fraters, M. R. T., & Billen, M. I. (2021).
+On the implementation and usability of crystal preferred orientation evolution in geodynamic modeling. Geochemistry, Geophysics, Geosystems, 22(10), e2021GC009846.
+
+- Hansen, L. N., Conrad, C. P., Boneh, Y., Skemer, P., Warren, J. M., & Kohlstedt, D. L. (2016).
+Viscous anisotropy of textured olivine aggregates: 2. Micromechanical model. Journal of Geophysical Research: Solid Earth, 121(10), 7137–7160.
+
+- Hirth, G., & Kohlstedt, D. (2004).
+Rheology of the upper mantle and the mantle wedge: A view from the experimentalists. In J. Eiler (Ed.), Geophysical Monograph Series (Vol. 138, pp. 83–105). American Geophysical Union.
+
+- Hill, R. (1948).
+A theory of the yielding and plastic flow of anisotropic metals. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 193(1033), 281–297.
+
+- Kaminski, E., Ribe, N. M., & Browaeys, J. T. (2004).
+D-Rex, a program for calculation of seismic anisotropy due to crystal lattice preferred orientation in the convective upper mantle. Geophysical Journal International, 158(2), 744–752.
+
+- Signorelli, J., Hassani, R., Tommasi, A., & Mameri, L. (2021).
+An effective parameterization of texture-induced viscous anisotropy in orthotropic materials with application for modeling geodynamical flows. Journal of Theoretical, Computational and Applied Mechanics, 6737.


I know you copied this from another cookbook, but these references should not be necessary. Make sure all the entries are part of references.bib and if they are correctly cited they will appear in our reference list: https://aspect-documentation--6615.org.readthedocs.build/en/6615/references.html. If you think it is useful that these entries appear here, you can instead add citations of these papers here (instead of manually writing the whole entries that can get out of date compared to the references list).

gassmoeller · 2025-07-28T12:04:48Z

why is the cookbook called CPO_induced..., but the plugin is called LPO_...? Please pick one, or explain why this is intentional.

Also regarding the name: Is the 3D part important? Wouldnt CPO_induced_anisotropic_viscosity.cc be a more explanatory file name?

gassmoeller · 2025-07-28T12:12:06Z

+    /*template <int matrix_size>
+    void check_eigenvalues_positive(const SymmetricTensor<2,matrix_size> &matrix);*/


remove these lines

gassmoeller · 2025-07-28T12:19:32Z

Please rename this file to anisotropic_stress.cc and the same for the header and class name. Let's avoid arbitrary abbreviations wherever possible, it makes the code easier to read.

gassmoeller · 2025-07-28T12:25:03Z

+    {
+
+      /**
+       * Computes the Bingham average of the CPO particle properties.


Is this comment still correct? How is this plugin different from the one included in the main code? (I mean what is the purpose of the Euler angles, why are they better than what is done in the other plugin?).

Wang-yijun · 2025-08-29T10:13:09Z

Very cool to have a cookbook for this feature. Sorry it took a while to get to this. I have left a first batch of comments, I havent yet looked at the code implementation in details, because I left some questions that would be useful to discuss beforehand.

Please in addition to my comments:

Add a test case in tests/ that runs the added cookbook reasonably fast. You can see how this is done in tests/magnetic_stripes.cc/prm.

Add a changelog entry in doc/modules

Let me know when someone can take another look.

Hi Rene, thanks for the feedback! I've addressed the comments above, added a test case in tests/, and included a changelog entry.
The automatic tests here failed, but I think it’s due to the setup not being linked to a dealii version with scalapack. On my local build, the cpo_induced_anisotropic_viscosity test passes.

gassmoeller

Another batch of comments. I am not quite done, but have to stop review for today, I hope I can get back to it tomorrow.

gassmoeller · 2025-09-02T11:22:32Z

Please also add the statistics file to the test output. This is useful because statistics reports the numbers with higher accuracy than the screen output (more digits).

gassmoeller · 2025-09-02T11:26:31Z

+
+## Introduction
+
+Individual crystals of the mineral olivine reorganized their orientations into the crystal-preferred orientations (CPO) under deformation. The directional-dependence (anisotropy) in the viscous properties of olivine crystals suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This CPO-induced anisotropic viscosity material model computes anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and modifies the subsequent deformation.


Suggested change

Individual crystals of the mineral olivine reorganized their orientations into the crystal-preferred orientations (CPO) under deformation. The directional-dependence (anisotropy) in the viscous properties of olivine crystals suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This CPO-induced anisotropic viscosity material model computes anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and modifies the subsequent deformation.

Individual crystals of the mineral olivine reorganize their orientations into the crystal-preferred orientations (CPO) under deformation. The viscous properties of olivine crystals are direction-dependent (anisotropic), which suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This cookbook model computes an anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and includes this information in the subsequent modeling process.

gassmoeller · 2025-09-02T11:32:55Z

+Our constitutive equation for the relationship between the strain rate and stress using the anisotropic viscosity tensor is adapted from {cite:t}`signorelli:etal:2021`:
+
+```{math}
+:label: eqn:one


Please use names specific to the equations, these labels can be referenced from anywhere in the documentation, and one is not self-explanatory when called from another place in the documentation. What about anisotropic_general_stress, and then modify the other labels below accordingly.

gassmoeller · 2025-09-02T11:40:44Z

+
+We compute the rotation matrix R on the particles and further convert it to Euler angles for computation and memory efficiency. These properties need to be interpolated from particles to fields to be used in the material model. As a result, the anisotropic viscosity material model requires at least one particle in each cell so that all cells can have the texture parameters (Euler angles and eigenvalues) for constructing the rotation matrix R and compute the Hill coefficients. In the material model, the interpolated Euler angles are converted to the rotation matrix again. We use the same notation R to describe the rotation matrix used in the material model in the following paragraphs.
+
+The inverse of the A tensor then needs to be rotated to the model reference frame. Since $inv(A)$ is the Kelvin notation of the rank-4 tensor, we apply the Kelvin notation representation of the R rotation matrix, $R_K$, on $inv(A_{ij})$:


Did you decide on purpose to keep using inv(A) instead of $A^{-1}$? Then please explain why. $A^{-1}$ is what is usually used to denote the inverse of matrix $A$, but of course if there is a reason to use something else that is fine (I just dont see the reason at the moment).

gassmoeller · 2025-09-02T11:45:28Z

+
+$$\sigma_{ij} = \frac{1}{\gamma J(R\prime*\sigma_{ij}* R)^{(n-1)}} * (R_K * inv(A_{ij}) * R_K\prime) * \dot\varepsilon_i$$
+
+We save $\frac{1}{\gamma J(R\prime*\sigma_{ij}* R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The tensorial part of anisotropic viscosity, $R_K * inv(A_{ij}) * R_K\prime$ is called the stress-strain director and is stored in the additional outputs. Due to the dependence of the scalar viscosity on stress and the stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is not stable and its value naturally fluctuates up and down. This fluctuation ultimately causes a numerical instability associated with the prediction of the anisotropic viscosity. Therefore, we damp the scalar viscosity using a non-linear Newtonian iteration, and in each iteration, only half of the change in the scalar viscosity is applied until the result is stable.


I think you have not addressed my earlier comment: https://github.com/geodynamics/aspect/pull/6615/files#r2236059538, and you still do not explain what the prime in $R_K^\prime$ stands for.

gassmoeller · 2025-09-02T12:50:00Z

+                            MaterialModel::MaterialModelOutputs<3> &out) const
+    {
+      const int dim=3;
+      const double sqrt2 = sqrt(2);


no need for this line, you can use numbers:SQRT2 instead (it is declared in #include <deal.II/base/numbers.h>.

gassmoeller · 2025-09-02T12:51:00Z

+
+
+


lets use 1 empty line for separating code inside functions. Also see here, for more coding conventions: https://dealii.org/developer/doxygen/deal.II/CodingConventions.html

gassmoeller · 2025-09-02T12:51:13Z

+
+      for (unsigned int q=0; q<in.n_evaluation_points(); ++q)
+        {
+          //change these according to diffusion dislocation material model I guess


Suggested change

//change these according to diffusion dislocation material model I guess

// Change these according to diffusion dislocation material model I guess

gassmoeller · 2025-09-02T12:51:20Z

+          out.entropy_derivative_pressure[q] = eos_outputs.entropy_derivative_pressure[0];
+          out.entropy_derivative_temperature[q] = eos_outputs.entropy_derivative_temperature[0];
+
+          //Calculate effective viscosity


Suggested change

//Calculate effective viscosity

// Calculate effective viscosity

gassmoeller · 2025-09-02T13:13:05Z

+              stress = 2*scalar_viscosity * V_r4 * deviatoric_strain_rate / 1e6; // Use stress in MPa //Added *2 2025.06.18.
+
+              while (std::abs(residual) > threshold && n_iterations < max_iteration)
+                {
+                  stress = (1./2.) * (stress + 2*scalar_viscosity * V_r4 * deviatoric_strain_rate / 1e6); //Added *2 2025.06.18.


Please remove the comments. Presumably you tested that the *2 is the correct thing to do? Then the comment is not really helpful.

Wang-yijun · 2025-09-29T19:26:25Z

Another batch of comments. I am not quite done, but have to stop review for today, I hope I can get back to it tomorrow.

Hi Rene, I've went through your comments. Looking forward for the next batch!

gassmoeller

Next batch of comments, mostly on the .prm file and on the anisotropic viscosity assemblers. I will try to get to the heating plugin and the rest of the .cc file next, then you can simplify this PR further.

gassmoeller · 2025-10-23T12:57:29Z

@@ -0,0 +1,133 @@
+(sec:cookbooks:CPO_induced_anisotropic_viscosity)=


We just introduced a new system to mark the features used in cookbooks and benchmarks in #6713, could you add tags here as well? I think the following should be appropriate:

```{tags} category:cookbook feature:3d feature:cartesian feature:modular-equations feature:compositional-fields feature:particles

Sorry for the extra effort.

gassmoeller · 2025-10-23T13:01:51Z

+                                                  "in the compressible case. If elasticity is used, the "
+                                                  "elastic contribution is being accounted for. The shear "
+                                                  "stress differs from the full stress tensor "
+                                                  "by the absence of the pressure. Note that the convention "
+                                                  "of positive compressive stress is followed. ")


This part needs an update. Instead of elasticity you know consider anisotropic viscosity. Also mention that this postprocessor only works with the specific material model of this cookbook.

gassmoeller · 2025-10-23T13:02:52Z

+                                                  "A visualization output object that generates output "
+                                                  "for the 6 (in 3d) components of the shear stress "
+                                                  "tensor, i.e., for the components of the tensor "
+                                                  "$-2\\eta\\varepsilon(\\mathbf u)$ "
+                                                  "in the incompressible case and "
+                                                  "$-2\\eta\\left[\\varepsilon(\\mathbf u)-"
+                                                  "\\tfrac 13(\\textrm{tr}\\;\\varepsilon(\\mathbf u))\\mathbf I\\right]$ "


Please update this to be consistent with the explanation you give in anisotropic_stress.h:38-47.

gassmoeller · 2025-10-23T13:03:13Z

+       * \frac 13 (\nabla \cdot \mathbf u) \mathbf I)$ * stress_strain_directors, and differs from the
+       * full stress by the absence of the pressure.  The second term in the
+       * difference is zero if the model is incompressible.
+       * If elasticity is used, the elastic contribution is being accounted for.


Is this true? Where is elasticity accounted for?

gassmoeller · 2025-10-23T13:10:49Z

+        void initialize() override;
+        void evaluate (const MaterialModel::MaterialModelInputs<dim> &in,
+                       MaterialModel::MaterialModelOutputs<dim> &out) const override;
+        static void declare_parameters (ParameterHandler &prm);
+        void parse_parameters (ParameterHandler &prm) override;
+        bool is_compressible () const override;
+        void create_additional_named_outputs(MaterialModel::MaterialModelOutputs<dim> &out) const override;


Please add an empty line between each function declaration. Since they are all default functions of the material model, you dont necessarily need documentation here, but empty lines make it easier to read. Also put them in the same order as in the Simple material model, except initialize at the beginning (we are not consistent about the position of create_additional_named_outputs, put it after initialize or after parse_parameters.

Then sort the functions in the .cc file in the same order. This doesnt matter for functionality, but matters a lot for simplifying the maintenance of ASPECT over time.

gassmoeller · 2025-10-23T13:26:17Z

+            set Mobility = 10 ## Annotation CPO Mobility
+            set Stress exponents = 3.5
+            set Exponents p = 1.5
+            set Threshold GBS = 0.3 ## Annotation TGBS


Suggested change

set Threshold GBS = 0.3 ## Annotation TGBS

set Threshold GBS = 0.3

gassmoeller · 2025-10-23T13:26:26Z

+            set Threshold GBS = 0.3 ## Annotation TGBS
+        end
+    end
+    subsection CPO Bingham Average


Suggested change

subsection CPO Bingham Average

subsection CPO Bingham Average

gassmoeller · 2025-10-23T13:27:02Z

+end
+
+subsection Solver parameters
+  set Temperature solver tolerance = 1e-10


Is the temperature solver tolerance change necessary? It is a setting that is rather specific. If so, document why you changed it.

gassmoeller · 2025-10-23T13:27:56Z

+# subsection Checkpointing
+#   set Steps between checkpoint = 50
+# end


gassmoeller · 2025-10-23T15:39:30Z

+  {
+    using namespace dealii;
+
+    /**


When #6728 is merged you can delete lines 83 - 330 and instead include the header aspect/simulator/assemblers/stokes_anisotropic_viscosity.h. That should simplify your plugin.

gassmoeller

A few more simplifications that are possible now that #6728 is merged.

And you need to indent and fix your tests. Let me know when you want me to take another look.

The test complains about ProcessGrid, are you missing a header (https://github.com/geodynamics/aspect/actions/runs/19302883915/job/55202043862?pr=6615#step:6:536)? Or is our tester not configured with Scalapack?

gassmoeller · 2025-11-13T15:45:24Z

+  namespace HeatingModel
+  {
+    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;
+    };
+  }


With the changes in #6728 you can now also delete this part. The shear heating plugin for anisotropic viscosity is now a normal heating plugin that you can activate in the Heating model subsection using set List of model names = anisotropic shear heating.
So in essence I dont think this header file is still necessary at all.

gassmoeller · 2025-11-13T15:46:54Z

+ */
+
+#include "cpo_induced_anisotropic_viscosity.h"
+#include "../../cookbooks/anisotropic_viscosity/av_material.h"


this line shouldnt be necessary any more

gassmoeller · 2025-11-13T15:47:17Z

+  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);
+    }
+  }
+


The whole namespace heating model can be deleted now as it is now part of the base ASPECT.

gassmoeller · 2025-11-13T15:48:16Z

+  namespace HeatingModel
+  {
+    ASPECT_REGISTER_HEATING_MODEL(ShearHeatingAnisotropicViscosity,
+                                  "anisotropic shear heating for cpo_induced_anisotropic_viscosity",
+                                  "Implementation of a standard model for shear heating. "
+                                  "Adds the term: "
+                                  "$  2 \\eta \\left( \\varepsilon - \\frac{1}{3} \\text{tr} "
+                                  "\\varepsilon \\mathbf 1 \\right) : \\left( \\varepsilon - \\frac{1}{3} "
+                                  "\\text{tr} \\varepsilon \\mathbf 1 \\right)$ to the "
+                                  "right-hand side of the temperature equation.")
+  }


delete this too

gassmoeller · 2025-11-13T15:48:34Z

+  namespace Assemblers
+  {
+#define INSTANTIATE(dim) \
+  template class StokesPreconditionerAnisotropicViscosity<dim>; \
+  template class StokesIncompressibleTermsAnisotropicViscosity<dim>; \
+
+    ASPECT_INSTANTIATE(INSTANTIATE)
+  }


everything here can be deleted as well.

Wang-yijun · 2025-11-17T11:07:12Z

Hi Rene! I finished responding to your last batch of comments.
For the test, I think the tester version of dealii doesn't include scalapack. I only see this error when I try to run the test without scalapack included.

gassmoeller

Here is the next batch of comments ;-). On the bright side now that the PR is significantly shorter I made it through the entire set of changes, so now I have seen all of the code. When you have addressed my comments the next reviews should be significantly easier and we are getting closer to a cookbook we can include in ASPECT.
Let me know when you want me to take another look.

gassmoeller · 2025-11-18T11:05:15Z

 end

+subsection Heating model
+  set List of model names = anisotropic shear heating


Good find, I must have missed that in my PR.

gassmoeller · 2025-11-19T13:23:38Z

 end

+subsection Heating model
+  set List of model names = anisotropic shear heating


This seem reasonable, but it was not active for the cookbook before my PR either. Did you check with @KiralyAgi if this cookbook was supposed to have shear heating? I dont remember the details, and it is not mentioned in the cookbook parameter file or documentation.

gassmoeller · 2025-11-19T13:26:56Z

+
+## Introduction
+
+Individual crystals of the mineral olivine reorganize their orientations into the crystal-preferred orientations (CPO) under deformation. The viscous properties of olivine crystals are direction-dependent (anisotropic), which suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This cookbook model computes an anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and includes this information in the subsequent modeling process.


Suggested change

Individual crystals of the mineral olivine reorganize their orientations into the crystal-preferred orientations (CPO) under deformation. The viscous properties of olivine crystals are direction-dependent (anisotropic), which suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This cookbook model computes an anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and includes this information in the subsequent modeling process.

Individual crystals of the mineral olivine reorganize their orientations into crystal-preferred orientations (CPO) under deformation. The viscous properties of olivine crystals are direction-dependent (anisotropic), which suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This cookbook model computes an anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and includes this information in the subsequent modeling process.

gassmoeller · 2025-11-19T13:39:47Z

+\sigma_{ij} = \frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}} * (R_K * A_{ij}^{-1} * R_K') * \dot\varepsilon_i
+```
+
+$R'$ and $R_K'$ is the transpose of matrix $R$ and $R_K$ respectively. We save $\frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The tensorial part of anisotropic viscosity, $R_K * A_{ij}^{-1} * R_K'$ is called the stress-strain director and is stored in the additional outputs. Due to the dependence of the scalar viscosity on stress and the stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is not stable and its value naturally oscillates. This fluctuation ultimately causes a numerical instability associated with the prediction of the anisotropic viscosity. Therefore, we damp the scalar viscosity using a non-linear Newton iteration, and in each iteration, only half of the change in the scalar viscosity is applied until the result is stable.


I would simplify a bit, and rephrase slightly:

Suggested change

$R'$ and $R_K'$ is the transpose of matrix $R$ and $R_K$ respectively. We save $\frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The tensorial part of anisotropic viscosity, $R_K * A_{ij}^{-1} * R_K'$ is called the stress-strain director and is stored in the additional outputs. Due to the dependence of the scalar viscosity on stress and the stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is not stable and its value naturally oscillates. This fluctuation ultimately causes a numerical instability associated with the prediction of the anisotropic viscosity. Therefore, we damp the scalar viscosity using a non-linear Newton iteration, and in each iteration, only half of the change in the scalar viscosity is applied until the result is stable.

$R'$ and $R_K'$ is the transpose of matrix $R$ and $R_K$ respectively. We save $\frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The tensorial part of anisotropic viscosity, $R_K * A_{ij}^{-1} * R_K'$ is called the stress-strain director and is stored in the additional outputs. Because the scalar viscosity depends on stress and stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is potentially unstable and its value can oscillate across nonlinear iterations. This oscillation can ultimately cause a numerical instability. Therefore, we damp the scalar viscosity computation using a non-linear Newton iteration, by applying only half of the change in each iteration until the result converges.

Actually after reading the code I dont think you perform a Newton iteration. Your algorithm is approximately:

while (residual > threshold) new_viscosity = 0.5 * (old_viscosity + new_estimated_viscosity);

This is a damped fixed point iteration, not a Newton iteration, which would use the derivative to determine the new viscosity. Please check and if necessary rephrase.

gassmoeller · 2025-11-19T13:48:35Z

+            SymmetricTensor<2,dim> aniso_stress;
+            aniso_stress= 2.*eta*deviatoric_strain_rate*anisotropic_viscosity->stress_strain_directors[q];


Just improving readability a bit and making the variable const:

Suggested change

SymmetricTensor<2,dim> aniso_stress;

aniso_stress= 2.*eta*deviatoric_strain_rate*anisotropic_viscosity->stress_strain_directors[q];

const SymmetricTensor<2,dim> anisotropic_stress =

2. * eta * deviatoric_strain_rate*anisotropic_viscosity->stress_strain_directors[q];

gassmoeller · 2025-11-19T16:00:43Z

+
+          prm.declare_entry ("Coefficients and intercept for F", "1.0390459583037057,  -0.767458622,  0.003066208,  0.19651133418307049,  0.413093763,  0.015463162,  -0.935925291,  -2.392877563,  0.051834768, 1.0799807050187482",
+                             Patterns::List(Patterns::Double()),
+                             "6 Coefficients and 1 intercept to compute the Hill Parameter F.");


This parameter and the other Hill coefficient input parameters need more explanation, as from the description of the cookbook and the class description I had expected these to be single value constants. In particular this description should contain the sentence from the cookbook description (that they are parameters of a regression analysis, and where one can find the source of the values (is that Agi's paper in review? is there a preprint available?). Also at least mention what the coefficients depend on, i.e. what are the coefficients multiplied with, or what do they represent?

gassmoeller · 2025-11-19T16:01:47Z

+          prm.declare_entry ("Grain size", "1000",
+                             Patterns::Double(),
+                             "Olivine anisotropic viscosity is dependent of grain size. Value is given in microns");


Please change to expect the value in meters. We are trying to avoid arbitrary units and let as many parameters as possible be provided in a consistent unit system.

gassmoeller · 2025-11-19T16:02:25Z

+
+

Remove two of the empty lines here to conform to our usual coding style.

gassmoeller · 2025-11-19T16:03:35Z

+  {
+    ASPECT_REGISTER_MATERIAL_MODEL(CPO_AV_3D,
+                                   "CPO-induced Anisotropic Viscosity",
+                                   "Olivine CPO related viscous anisotropy based on the Simple material model")


Also add "For more details see the documentation of this class in its header file and the cookbook documentation."

gassmoeller · 2025-11-19T16:04:03Z

+  namespace MaterialModel
+  {
+    ASPECT_REGISTER_MATERIAL_MODEL(CPO_AV_3D,
+                                   "CPO-induced Anisotropic Viscosity",


to be consistent with our other plugins, please use:

Suggested change

"CPO-induced Anisotropic Viscosity",

"CPO-induced anisotropic viscosity",

bangerth · 2025-12-11T15:37:21Z

+
+## Compiling requirement
+
+Since the determinant of $A_{ij}$ is 0, $A_{ij}$ is a singular, non-invertible matrix. We find the MoorePenrose pseudoinverse of the matrix $A_{ij}$ as an approximate of the inverse of $A_{ij}$, using the SCALAPACK package provided in deal.II. Thus it is required to link ASPECT with a deal.II version with the scalapack package included in order to run this cookbook. When using candi you can enable the scalapack package by including `scalapack` in the list of installed packages or uncommenting the line in `candi.cfg` that corresponds to the scalapack installation. Alternatively, you can install ScaLAPACK yourself and enable the configuration variable `DEAL_II_WITH_SCALAPACK` during the cmake configuration of deal.II.


How large is this matrix, and is it stored in parallel? I'm asking because SCALAPACK is such an old-fashioned package that I'd really like to avoid using it. If this is a matrix of modest size, I imagine we can find replacements for the SCALAPACK functionality that don't require us tying ASPECT to that package.

Separately, you have at least three different spellings of SCALAPACK in this paragraph -- would you mind picking one? :-)

6 by 6, see https://github.com/geodynamics/aspect/pull/6615/changes#diff-62582cf7a89944126803ad4a66c4db0cfff1bc34124706ca677cdaafd7326c13R252-R253

I see. That's small. I think SCALAPACK is overkill for that. Could we instead use something like https://stackoverflow.com/questions/69255272/what-are-the-best-driver-routines-in-lapack-for-implementing-the-moore-penrose-p or https://www.netlib.org/lapack/explore-html/da/d55/group__gelss_gac6159de3953ae0386c2799294745ac90.html for this? I don't recall whether we already require LAPACK, but I'd much rather require that than LAPACK.

We already make use of LAPACK for some of the particle interpolators, so it wouldnt be a limitation to use it here as well. (in particular the least squares interpolators)

I deleted the scalapack part and used compute_svd in lapack to create a helper function pseudoinverse. Now it doesn't require scalapack package anymore. Thanks for the link!

Wang-yijun · 2025-12-15T12:56:27Z

Here is the next batch of comments ;-). On the bright side now that the PR is significantly shorter I made it through the entire set of changes, so now I have seen all of the code. When you have addressed my comments the next reviews should be significantly easier and we are getting closer to a cookbook we can include in ASPECT. Let me know when you want me to take another look.

Hi Rene! I've responded to these comments and you can take a look any time now. I'll fix the second_invariant function and test later.

gassmoeller

Getting closer. My comments are mostly small code improvements.
Please:

Address the remaining comments
Update your test results with the reference tester results
Run the indent script
Rebase to the latest development branch and squash your commits into 1 (or a few logically separated) commits.

If all goes well this should then be ready.

gassmoeller · 2025-12-19T11:57:49Z

+
+```{math}
+:label: eqn:anisotropic_stress_final
+\sigma_{ij} = \frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}} * (R_K * A_{ij}^{-1} * R_K') * \dot\varepsilon_i


Are you missing an j at the final epsilon?

Suggested change

\sigma_{ij} = \frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}} * (R_K * A_{ij}^{-1} * R_K') * \dot\varepsilon_i

\sigma_{ij} = \frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}} * (R_K * A_{ij}^{-1} * R_K') * \dot\varepsilon_{ij}

gassmoeller · 2025-12-19T11:58:43Z

+\sigma_{ij} = \frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}} * (R_K * A_{ij}^{-1} * R_K') * \dot\varepsilon_i
+```
+
+$R'$ and $R_K'$ is the transpose of matrix $R$ and $R_K$ respectively. We save $\frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The scalar viscosity of the current step is also stored into the prescribed field, so that at the beginning of the next step, the anisotropic stress is computed with the scalar viscosity of the previous step using {math:numref}`eqn:anisotropic_stress`. The tensorial part of anisotropic viscosity, $R_K * A_{ij}^{-1} * R_K'$ is called the stress-strain director and is stored in the additional outputs. DBecause the scalar viscosity depends on stress and the stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is potentially unstable and its value can oscillate across nonlinear iterations. This oscillation can ultimately cause a numerical instability associated with the prediction of the anisotropic viscosity. Therefore, we damp the scalar viscosity computation using a fixed-point iteration, by applying only half of the change in each iteration until the result converges.


Suggested change

$R'$ and $R_K'$ is the transpose of matrix $R$ and $R_K$ respectively. We save $\frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The scalar viscosity of the current step is also stored into the prescribed field, so that at the beginning of the next step, the anisotropic stress is computed with the scalar viscosity of the previous step using {math:numref}`eqn:anisotropic_stress`. The tensorial part of anisotropic viscosity, $R_K * A_{ij}^{-1} * R_K'$ is called the stress-strain director and is stored in the additional outputs. DBecause the scalar viscosity depends on stress and the stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is potentially unstable and its value can oscillate across nonlinear iterations. This oscillation can ultimately cause a numerical instability associated with the prediction of the anisotropic viscosity. Therefore, we damp the scalar viscosity computation using a fixed-point iteration, by applying only half of the change in each iteration until the result converges.

$R'$ and $R_K'$ is the transpose of matrix $R$ and $R_K$ respectively. We save $\frac{1}{\gamma J(R'*\sigma_{ij}* R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The scalar viscosity of the current step is also stored into the prescribed field, so that at the beginning of the next step, the anisotropic stress is computed with the scalar viscosity of the previous step using {math:numref}`eqn:anisotropic_stress`. The tensorial part of anisotropic viscosity, $R_K * A_{ij}^{-1} * R_K'$ is called the stress-strain director and is stored in the additional outputs. Because the scalar viscosity depends on stress and the stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is potentially unstable and its value can oscillate across nonlinear iterations. This oscillation can ultimately cause a numerical instability associated with the prediction of the anisotropic viscosity. Therefore, we damp the scalar viscosity computation using a fixed-point iteration, by applying only half of the change in each iteration until the result converges.

gassmoeller · 2025-12-19T11:59:49Z

+## Model setup
+
+The usage of the AV material model is demonstrated with a 3d simple shear box model, where its dimension is $1 \times 1 \times 1 $ (non-dimensionalized). The shear strain rate is set to
+$0.5$. The origin is the center of the box, and one Olivine particle with 1000 grains sits at the origin to track CPO developments for computation of anisotropic viscosity parameters.


Suggested change

$0.5$. The origin is the center of the box, and one Olivine particle with 1000 grains sits at the origin to track CPO developments for computation of anisotropic viscosity parameters.

$0.5$. The origin is the center of the box, and one particle representing 1000 olivine grains sits at the origin to track CPO developments for computation of anisotropic viscosity parameters.

gassmoeller · 2025-12-19T12:00:30Z

+```{literalinclude} min_particles_per_cell.part.prm
+```
+
+- **CPO particle property**: The CPO particle property must be stored for use by the AV model. This requires enabling the particle and crystal preferred orientation postprocessors and the relevant subsuctions for them, including the CPO Bingham Average plugin, which calculates the Hill coefficients:


Suggested change

- **CPO particle property**: The CPO particle property must be stored for use by the AV model. This requires enabling the particle and crystal preferred orientation postprocessors and the relevant subsuctions for them, including the CPO Bingham Average plugin, which calculates the Hill coefficients:

- **CPO particle property**: The CPO particle property must be stored for use by the AV model. This requires enabling the particle and crystal preferred orientation postprocessors and the relevant subsections for them, including the CPO Bingham Average plugin, which calculates the Hill coefficients:

gassmoeller · 2025-12-19T12:02:36Z

+      {
+        const unsigned int n_quadrature_points = input_data.solution_values.size();
+        Assert (computed_quantities.size() == n_quadrature_points,    ExcInternalError());
+        Assert ((computed_quantities[0].size() == Tensor<2,dim>::n_independent_components),


Just to stay consistent in style to the other asserts:

Suggested change

Assert ((computed_quantities[0].size() == Tensor<2,dim>::n_independent_components),

Assert (computed_quantities[0].size() == Tensor<2,dim>::n_independent_components,

gassmoeller · 2025-12-19T12:19:30Z

+
+
+
+// std::explicit instantiations


this was probably unintended search and replace of exp:

Suggested change

// std::explicit instantiations

// explicit instantiations

gassmoeller · 2025-12-19T12:20:43Z

+        void pseudoinverse(LAPACKFullMatrix<double> &A,
+                           LAPACKFullMatrix<double> &A_pinv) const;


this function can move to the private section of this file, and also needs some documentation. What is the input? What is the output? Which algorithm did you use to implement it?

gassmoeller · 2025-12-19T12:21:43Z

@@ -0,0 +1,187 @@
+# This cookbook demonstrates how CPO induced anisotropic viscosity evolves during
+# simple shearing. The model consists of 1 cell and one particle in the middle.
+# The particle includes 500 grains with initially random orientations (isotropic


500 grains? below it says 1000

gassmoeller · 2025-12-19T12:22:17Z

+set Additional shared libraries = ./plugin/libCPO_induced_anisotropic_viscosity.debug.so
+set CFL number = 0.1
+set Timing output frequency  = 10000
+set Dimension                = 3
+set Pressure normalization   = surface
+set Surface pressure         = 0
+set Nonlinear solver scheme  = single Advection, single Stokes
+set End time                 = 0.1
+set Use years instead of seconds = false
+set Output directory  = output_Shearbox_CPO_AV


Align the equal signs to improve readability

gassmoeller · 2025-12-19T12:24:34Z

         */
        std::vector<SymmetricTensor<4,dim>> stress_strain_directors;
+
+        static Tensor<2,3> euler_angles_to_rotation_matrix(double phi1, double theta, double phi2);


I dont see this function implemented anywhere. Is this a leftover?

…cosity.

Wang-yijun · 2026-01-07T14:38:28Z

Getting closer. My comments are mostly small code improvements. Please:

Address the remaining comments

Update your test results with the reference tester results

Run the indent script

Rebase to the latest development branch and squash your commits into 1 (or a few logically separated) commits.

If all goes well this should then be ready.

Hello Rene, the PR is ready for review now.

gassmoeller

Hi @Wang-yijun I took the liberty of fixing some small remaining issues myself (you can see them in my last commits, mostly improving code style and addressing a few comments I had left that you missed). I also reworded the commit messages of the commits you shared with Agi, they were quite long, and not very expressive of what happened.
If the testers succeed this should be ready to merge. Thanks for finishing this PR, I think it is an interesting application case 👍

gassmoeller · 2026-01-19T15:49:33Z

I had to remove the memory statistics postprocessor from your test. Memory statistics are a bit unreliable one the test system, because the amount of allocated memory can vary based on harware/ operating system (I guess). Anyway this is now ready to merge. Thanks for all your efforts!

gassmoeller reviewed Jul 28, 2025

View reviewed changes

gassmoeller mentioned this pull request Sep 2, 2025

Move anisotropic viscosity assemblers into main code #6661

Closed

gassmoeller reviewed Sep 2, 2025

View reviewed changes

gassmoeller mentioned this pull request Oct 23, 2025

Unify anisotropic viscosity Stokes assemblers / Fix anisotropic viscosity cookbook #6728

Merged

gassmoeller reviewed Oct 23, 2025

View reviewed changes

Wang-yijun force-pushed the CPO_AV_Cookbook branch from a434d2a to eb17df8 Compare November 12, 2025 15:00

gassmoeller reviewed Nov 13, 2025

View reviewed changes

Wang-yijun force-pushed the CPO_AV_Cookbook branch 2 times, most recently from 274c5d6 to 41bad45 Compare November 16, 2025 01:03

gassmoeller reviewed Nov 19, 2025

View reviewed changes

This was referenced Dec 10, 2025

Update deal.II devel tester image #6808

Merged

Add scalapack package to docker images #6811

Closed

bangerth reviewed Dec 11, 2025

View reviewed changes

gassmoeller reviewed Dec 19, 2025

View reviewed changes

Created initial version of a cookbook for CPO induced anisotropic vis…

e56efac

…cosity.

Wang-yijun force-pushed the CPO_AV_Cookbook branch 3 times, most recently from b7a8d77 to 45845a2 Compare January 7, 2026 12:46

gassmoeller added the ready to test label Jan 7, 2026

KiralyAgi and others added 3 commits January 7, 2026 17:02

Add stress visualization. Rework cookbook.

95a723e

Address comments

1c5e27b

Minor cleanup

e8a446d

gassmoeller force-pushed the CPO_AV_Cookbook branch from 45845a2 to 8add42d Compare January 7, 2026 16:04

gassmoeller approved these changes Jan 7, 2026

View reviewed changes

gassmoeller force-pushed the CPO_AV_Cookbook branch from 8add42d to 0a73146 Compare January 9, 2026 11:19

Update test output

ad0633d

gassmoeller force-pushed the CPO_AV_Cookbook branch from 0a73146 to ad0633d Compare January 19, 2026 13:51

gassmoeller merged commit 29449ac into geodynamics:main Jan 19, 2026
8 of 9 checks passed


		## Introduction

		The directional-dependence (anisotropy) in the viscous properties of olivine crystals suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This CPO-induced anisotropic viscosity material model computes anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and modifies the subsequent deformation.


		$$\dot \varepsilon_{ij} = \gamma J(\sigma_{ij})^{(n-1)} A_{ij} \sigma_{ij}$$

		where $\gamma$ is the part of fluidity which is temperature- and grain-size dependent:


		$$\gamma=\gamma_0 exp \left(\frac{-Q}{RT} \right) /d^m$$

		$\gamma_0=1.1\times 10^{5}$ is the isotropic fluidity, $Q=530 kJ/mol$ is the activation energy, $R=8.314 m^3 \cdot Pa \cdot K^{−1} \cdot mol^{−1}$ is the gas constant, $d=0.001 m$ is the grain size, and $m=0.73$ is the grain size exponent. These values are taken from {cite:t}`hansen:etal:2016` and {cite:t}`HK04`. $J(\sigma_{ij})$ is the equivalent yield stress, where $\sigma_{ij}$ is the deviatoric (anisotropic) stress computed using the tensorial and scalar component of the anisotropic viscosity:


		$$J(\sigma_{ij})$$ and $A_{ij}$ are computed using Hill coefficients $H, J, K, L, M,$ and $N$ {cite}`hill:1948`, which are obtained from regression analysis of a texture database constructed with olivine textures from laboratory experiments, shear box models, and subduction models (Kiraly et al., in rev.).

		In this material model plugin, strain rate, density, temperature, and other parameters are taken as input to compute the viscosity, which is passed into the Stokes system to compute the stress (Figure 1). As a result, we adapt eq. (1) to be:

		/*template <int matrix_size>
		void check_eigenvalues_positive(const SymmetricTensor<2,matrix_size> &matrix);*/


		## Introduction

		Individual crystals of the mineral olivine reorganized their orientations into the crystal-preferred orientations (CPO) under deformation. The directional-dependence (anisotropy) in the viscous properties of olivine crystals suggests that the effective viscosity for olivine rocks/aggregates is different when deformations occur in different directions relative to the CPO. This CPO-induced anisotropic viscosity material model computes anisotropic viscosity based on the CPO evolution predicted by D-Rex ({cite}`fraters_billen_2021_cpo`; {cite}`kaminski2004`) and modifies the subsequent deformation.


		We compute the rotation matrix R on the particles and further convert it to Euler angles for computation and memory efficiency. These properties need to be interpolated from particles to fields to be used in the material model. As a result, the anisotropic viscosity material model requires at least one particle in each cell so that all cells can have the texture parameters (Euler angles and eigenvalues) for constructing the rotation matrix R and compute the Hill coefficients. In the material model, the interpolated Euler angles are converted to the rotation matrix again. We use the same notation R to describe the rotation matrix used in the material model in the following paragraphs.

		The inverse of the A tensor then needs to be rotated to the model reference frame. Since $inv(A)$ is the Kelvin notation of the rank-4 tensor, we apply the Kelvin notation representation of the R rotation matrix, $R_K$, on $inv(A_{ij})$:


		$$\sigma_{ij} = \frac{1}{\gamma J(R\prime\sigma_{ij} R)^{(n-1)}} * (R_K * inv(A_{ij}) * R_K\prime) * \dot\varepsilon_i$$

		We save $\frac{1}{\gamma J(R\prime\sigma_{ij} R)^{(n-1)}}$ as the material model viscosity output, which we call the scalar viscosity. The tensorial part of anisotropic viscosity, $R_K * inv(A_{ij}) * R_K\prime$ is called the stress-strain director and is stored in the additional outputs. Due to the dependence of the scalar viscosity on stress and the stress is determined by the scalar viscosity from the previous time step, the computation of the scalar viscosity is not stable and its value naturally fluctuates up and down. This fluctuation ultimately causes a numerical instability associated with the prediction of the anisotropic viscosity. Therefore, we damp the scalar viscosity using a non-linear Newtonian iteration, and in each iteration, only half of the change in the scalar viscosity is applied until the result is stable.

	//change these according to diffusion dislocation material model I guess
	// Change these according to diffusion dislocation material model I guess

	//Calculate effective viscosity
	// Calculate effective viscosity

		@@ -0,0 +1,133 @@
		(sec:cookbooks:CPO_induced_anisotropic_viscosity)=

	set Threshold GBS = 0.3 ## Annotation TGBS
	set Threshold GBS = 0.3

	subsection CPO Bingham Average

	subsection CPO Bingham Average

		SymmetricTensor<2,dim> aniso_stress;
		aniso_stress= 2.etadeviatoric_strain_rate*anisotropic_viscosity->stress_strain_directors[q];

	"CPO-induced Anisotropic Viscosity",
	"CPO-induced anisotropic viscosity",


		## Compiling requirement

		Since the determinant of $A_{ij}$ is 0, $A_{ij}$ is a singular, non-invertible matrix. We find the MoorePenrose pseudoinverse of the matrix $A_{ij}$ as an approximate of the inverse of $A_{ij}$, using the SCALAPACK package provided in deal.II. Thus it is required to link ASPECT with a deal.II version with the scalapack package included in order to run this cookbook. When using candi you can enable the scalapack package by including `scalapack` in the list of installed packages or uncommenting the line in `candi.cfg` that corresponds to the scalapack installation. Alternatively, you can install ScaLAPACK yourself and enable the configuration variable `DEAL_II_WITH_SCALAPACK` during the cmake configuration of deal.II.

	\sigma_{ij} = \frac{1}{\gamma J(R'\sigma_{ij} R)^{(n-1)}} * (R_K * A_{ij}^{-1} * R_K') * \dot\varepsilon_i
	\sigma_{ij} = \frac{1}{\gamma J(R'\sigma_{ij} R)^{(n-1)}} * (R_K * A_{ij}^{-1} * R_K') * \dot\varepsilon_{ij}

	$0.5$. The origin is the center of the box, and one Olivine particle with 1000 grains sits at the origin to track CPO developments for computation of anisotropic viscosity parameters.
	$0.5$. The origin is the center of the box, and one particle representing 1000 olivine grains sits at the origin to track CPO developments for computation of anisotropic viscosity parameters.

	- CPO particle property: The CPO particle property must be stored for use by the AV model. This requires enabling the particle and crystal preferred orientation postprocessors and the relevant subsuctions for them, including the CPO Bingham Average plugin, which calculates the Hill coefficients:
	- CPO particle property: The CPO particle property must be stored for use by the AV model. This requires enabling the particle and crystal preferred orientation postprocessors and the relevant subsections for them, including the CPO Bingham Average plugin, which calculates the Hill coefficients:

	Assert ((computed_quantities[0].size() == Tensor<2,dim>::n_independent_components),
	Assert (computed_quantities[0].size() == Tensor<2,dim>::n_independent_components,

Conversation

Wang-yijun commented Jul 17, 2025

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Wang-yijun commented Aug 29, 2025

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Wang-yijun commented Sep 29, 2025

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

		void pseudoinverse(LAPACKFullMatrix<double> &A,
		LAPACKFullMatrix<double> &A_pinv) const;