CPO-induced anisotropic viscosity cookbook

cookbooks/CPO_induced_anisotropic_viscosity/doc/CPO_induced_anisotropic_viscosity.md

@@ -0,0 +1,142 @@
+```{tags}
+category:cookbook
+feature:3d
+feature:cartesian
+feature:modular-equations
+feature:compositional-fields
+feature:particles
+```
+
+(sec:cookbooks:CPO_induced_anisotropic_viscosity)=
+# CPO induced anisotropic viscosity
+
+*This section was contributed by Yijun Wang and Ágnes Király.*
+
+This cookbook explains how to use the CPO-induced anisotropic viscosity material model to set up a shear box model.
+
+## 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.
+
+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:anisotropic_general_stress
+\dot{\varepsilon}_{ij} = \gamma J(\sigma_{ij})^{(n-1)} A_{ij} \sigma_{ij}
+```
+
+where $\gamma$ is the part of fluidity (the inverse of viscosity) which is temperature- and grain-size dependent:
+
+```{math}
+:label: eqn:fluidity
+\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 for olivine are taken from rock experiments performed by {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:
+
+```{math}
+:label: eqn:equivalent_yield_stress
+J(\sigma_{ij})=(F(\sigma_{11} - \sigma_{22})^2+G(\sigma_{22} - \sigma_{33})^2+H(\sigma_{33} - \sigma_{11})^2+2L\sigma_{23}^2+2M\sigma_{13}^2+2N\sigma_{12}^2)^{1/2}
+```
+
+and $A_{ij}$ is the anisotropic tensor of fluidity in Kelvin notation:
+
+```{math}
+:label: eqn:anisotropic_fluidity
+A_{ij}=\frac{2}{3} \left[
+\begin{matrix}
+F+H & -F & -H & 0 & 0 & 0 \\
+-F & G+F & -G & 0 & 0 & 0 \\
+-H & -G & H+G & 0 & 0 & 0 \\
+0 & 0 & 0 & L & 0 & 0 \\
+0 & 0 & 0 & 0 & M & 0 \\
+0 & 0 & 0 & 0 & 0 & N
+\end{matrix} \right]
+```
+
+$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 anisotropic viscosity, which is passed into the Stokes system to compute the stress. As a result, we adapt {math:numref}`eqn:anisotropic_general_stress` to be:
+
+```{math}
+:label: eqn:anisotropic_stress
+\sigma_{ij} = \frac{1}{\gamma J(\sigma_{ij})^{(n-1)}} * A_{ij}^{-1} * \dot\varepsilon_{ij}
+```
+
+Since the Hill coefficients are defined in the microscopic CPO reference frame, and parameters computed in ASPECT are in the macroscopic model reference frame, several reference frame conversions are needed. First, we need to rotate $\sigma_{ij}$ in {math:numref}`eqn:equivalent_yield_stress` from the model reference frame to the CPO reference frame so that $J(\sigma_{ij})$ is in the CPO reference frame. This is achieved by constructing a matrix from the eigenvectors corresponding with the largest eigenvalues of the covariance matrix for the a-, b-, and c-axis of olivine textures and then we assign the nearest orthogonal matrix to be the rotation matrix R:
+
+```{math}
+:label: eqn:rotation_matrix
+R = \left[
+\begin{matrix}
+\verb|max_eigenvector|\_{a1} & \verb|max_eigenvector|\_{b1} & \verb|max_eigenvector|\_{c1} \\
+\verb|max_eigenvector|\_{a2} & \verb|max_eigenvector|\_{b2} & \verb|max_eigenvector|\_{c2} \\
+\verb|max_eigenvector|\_{a3} & \verb|max_eigenvector|\_{b3} & \verb|max_eigenvector|\_{c3} \\
+\end{matrix} \right]
+```
+
+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 $A_{ij}^{-1}$ is the Kelvin notation of the rank-4 tensor, we apply the Kelvin notation representation of the R rotation matrix, $R_K$, on $A_{ij}^{-1}$:
+
+```{math}
+:label: eqn:rotation_matrix_kelvin
+R_K = \left[
+\begin{matrix}
+R_{11}^2 & R_{12}^2 & R_{13}^2 & \sqrt2* R_{12}* R_{13} & \sqrt2* R_{11}* R_{13} & \sqrt2* R_{11}* R_{12} \\
+R_{21}^2 & R_{22}^2 & R_{23}^2 & \sqrt2* R_{22}* R_{23} & \sqrt2* R_{21}* R_{23} & \sqrt2* R_{21}* R_{22} \\
+R_{31}^2 & R_{32}^2 & R_{33}^2 & \sqrt2* R_{32}* R_{33} & \sqrt2* R_{31}* R_{33} & \sqrt2* R_{31}* R_{32} \\
+\sqrt2* R_{21}* R_{31} & \sqrt2* R_{23}* R_{32} & \sqrt2* R_{23}* R_{33} & R_{22}* R_{33}+R_{23}* R_{32} & R_{21}* R_{33}+R_{23}* R_{31} & R_{21}* R_{32}+R_{22}* R_{31} \\
+\sqrt2* R_{11}* R_{31} & \sqrt2* R_{12}* R_{32} & \sqrt2* R_{13}* R_{33} & R_{12}* R_{33}+R_{13}* R_{32} & R_{11}* R_{33}+R_{13}* R_{31} & R_{11}* R_{32}+R_{12}* R_{31} \\
+\sqrt2* R_{11}* R_{21} & \sqrt2* R_{12}* R_{22} & \sqrt2* R_{13}* R_{23} & R_{12}* R_{23}+R_{13}* R_{22} & R_{11}* R_{23}+R_{13}* R_{32} & R_{11}* R_{22}+R_{12}* R_{21}
+\end{matrix} \right]
+```
+
+The final equation involving all reference frame conversions is:
+
+```{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
+```
+
+$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.
+
+
+
+## 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.
+
+## 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.
+
+Since the AV material model computes viscosity based on the evolving CPO stored on particles, several setup requirements must be met:
+
+- **Particles per cell**: Each computational cell must contain at least one particle, to allow interpolation of the CPO particle property. This is achieved by setting (in the Particles subsection):
+
+```{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:
+
+```{literalinclude} cpo_particle_property.part.prm
+```
+
+Note: These settings are similar to those used for simulations involving CPO alone. However, for the AV model, it is essential to set `Use rotation matrix = false` in the CPO Bingham Average subsection, so that the CPO is represented using Euler angles, as required.
+
+- **Compositional fields**: The eigenvalues and Euler angles of the CPO tensor are stored in compositional fields. This requires the following input file section:
+
+```{literalinclude} compositional_field.part.prm
+```
+
+In the `CPO induced Anisotropic Viscosity` material model subsection, all parameters have reasonable default values and do not need to be manually specified unless customization is needed.
+
+This shear box model uses an additional postprocessor, anisotropic stress, which is also implemented in this cookbook. It outputs a 3-by-3 matrix that can be visualized as a tensor, similar to the standard stress postprocessor. With the anisotropic viscosity material model, applying simple shear produces deformation in multiple directions. As a result, the anisotropic stress tensor appears as elongated and slightly tilted glyphs (indicating the principal stress directions), in contrast to the isotropic stress tensor (see figure below).
+
+```{figure-md} fig:anisotropic_stress_shearbox
+<img src="anisotropic_stress.png" style="width:100.0%" />
+
+Expected output of the shear box model using anisotropic viscosity material model, showing the anisotropic stress and stress postprocessor as tensor glyphs (blue disks) in Paraview. The arrows indicate the direction and magnitude of velocity.
+```

cookbooks/CPO_induced_anisotropic_viscosity/doc/compositional_field.part.prm

@@ -0,0 +1,32 @@
+subsection Initial composition model
+  set List of model names = function
+  subsection Function
+    set Function expression = 0;\
+                              0; 0; 0; 0;\
+                              0; 0; 0; 0;\
+                              0; 0; 0; 0;\
+                              0; 0
+  end
+end
+
+subsection Compositional fields
+  set Number of fields = 15
+  set Names of fields = scalar_vis, \
+                        phi1, eigvalue_a1, eigvalue_a2, eigvalue_a3,\
+                        theta, eigvalue_b1, eigvalue_b2, eigvalue_b3,\
+                        phi2, eigvalue_c1, eigvalue_c2, eigvalue_c3,\
+                        D, water
+  set Compositional field methods = prescribed field, \
+                                    particles, particles, particles, particles, \
+                                    particles, particles, particles, particles, \
+                                    particles, particles, particles, particles, \
+                                    static, static
+  set Mapped particle properties = phi1:cpo mineral 0 phi1[0], eigvalue_a1:cpo mineral 0 eigenvalues a axis[0], eigvalue_a2:cpo mineral 0 eigenvalues a axis[1], eigvalue_a3:cpo mineral 0 eigenvalues a axis[2], \
+                                   theta:cpo mineral 0 theta[0], eigvalue_b1:cpo mineral 0 eigenvalues b axis[0], eigvalue_b2:cpo mineral 0 eigenvalues b axis[1], eigvalue_b3:cpo mineral 0 eigenvalues b axis[2], \
+                                   phi2:cpo mineral 0 phi2[0], eigvalue_c1:cpo mineral 0 eigenvalues c axis[0], eigvalue_c2:cpo mineral 0 eigenvalues c axis[1], eigvalue_c3:cpo mineral 0 eigenvalues c axis[2]
+  set Types of fields = generic, \
+                        generic, generic, generic, generic, \
+                        generic, generic, generic, generic, \
+                        generic, generic, generic, generic, \
+                        generic, generic
+end

cookbooks/CPO_induced_anisotropic_viscosity/doc/cpo_particle_property.part.prm

@@ -0,0 +1,60 @@
+subsection Postprocess
+    set List of postprocessors = velocity statistics, composition statistics, memory statistics, visualization, particles, crystal preferred orientation
+
+    subsection Visualization
+        set Time between graphical output = 0.1
+        set List of output variables = material properties, strain rate, named additional outputs, shear stress, stress
+
+      subsection Material properties
+        set List of material properties = density, viscosity
+      end
+    end
+    subsection Particles
+        set Time between data output = 0.1
+        set Data output format       = gnuplot, vtu
+        set Exclude output properties = rotation_matrix, volume fraction
+    end
+    subsection Crystal Preferred Orientation
+        set Time between data output = 0.1
+        set Write in background thread = true
+        set Compress cpo data files = false
+        set Write out raw cpo data = mineral 0: volume fraction, mineral 0: Euler angles
+        set Write out draw volume weighted cpo data = mineral 0: volume fraction, mineral 0: Euler angles
+    end
+end
+
+subsection Particles
+    set List of particle properties =  integrated strain, integrated strain invariant, velocity, pT path, strain rate, crystal preferred orientation, cpo bingham average, cpo elastic tensor
+    set Allow cells without particles = false
+    set Interpolation scheme = nearest neighbor
+    set Minimum particles per cell = 1
+    set Maximum particles per cell = 5
+    set Load balancing strategy     = add particles
+    set Particle generator name = ascii file
+    subsection Generator
+        subsection Ascii file
+            set Data directory = ./
+            set Data file name = particle_one.dat
+        end
+    end
+    subsection Crystal Preferred Orientation
+        subsection Initial grains
+            set Minerals = Olivine: A-fabric
+            set Volume fractions minerals = 1.0
+        end
+        set Number of grains per particle = 1000
+        set Property advection method = Backward Euler
+        set Property advection tolerance = 1e-15
+        set CPO derivatives algorithm = D-Rex 2004
+        subsection D-Rex 2004
+            set Mobility = 10
+            set Stress exponents = 3.5
+            set Exponents p = 1.5
+            set Threshold GBS = 0.3
+        end
+    end
+    subsection CPO Bingham Average
+        set Random number seed = 200
+        set Use rotation matrix = false
+    end
+end

cookbooks/CPO_induced_anisotropic_viscosity/doc/min_particles_per_cell.part.prm

@@ -0,0 +1,2 @@
+set Minimum particles per cell  = 1
+set Load balancing strategy     = add particles

cookbooks/CPO_induced_anisotropic_viscosity/particle_one.dat

@@ -0,0 +1 @@
+0.5   0.5   0.5

cookbooks/CPO_induced_anisotropic_viscosity/plugin/CMakeLists.txt

@@ -0,0 +1,39 @@
+# Copyright (C) 2011 - 2024 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/>.
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.13.4)
+
+FIND_PACKAGE(Aspect 2.4.0 QUIET HINTS ${Aspect_DIR} ../ ../../ $ENV{ASPECT_DIR})
+
+IF (NOT Aspect_FOUND)
+  MESSAGE(FATAL_ERROR "\n"
+  "Could not find a valid ASPECT build/installation directory. "
+  "Please specify the directory where you are building ASPECT by passing\n"
+  "   -D Aspect_DIR=<path to ASPECT>\n"
+  "to cmake or by setting the environment variable ASPECT_DIR in your shell "
+  "before calling cmake. See the section 'How to write a plugin' in the "
+        "manual for more information.")
+ENDIF ()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+
+SET(TARGET "CPO_induced_anisotropic_viscosity")
+PROJECT(${TARGET})
+
+ADD_LIBRARY(${TARGET} SHARED  cpo_induced_anisotropic_viscosity.h cpo_induced_anisotropic_viscosity.cc anisotropic_stress.cc anisotropic_stress.h)
+ASPECT_SETUP_PLUGIN(${TARGET})