up

docs/make.jl

@@ -79,7 +79,12 @@ makedocs(;
         "User guide"=> Any[
             "Installation" => "man/installation.md",
             "Backend"   => "man/backend.md",
-            "Equations" => "man/equations.md",
+            "Equations" => Any[
+                "Governing equations"    => "man/equations_basic.md",
+                "Constitutive equations" => "man/constitutive_equations.md",
+                "APT equations"          => "man/equations_APT.md",
+                "Discretization"         => "man/equations_discretization.md",
+            ],
             "Boundary conditions" => "man/boundary_conditions.md",
             "Advection" => "man/advection.md",
             ],

docs/src/man/constitutive_equations.md

@@ -17,12 +17,37 @@ $$
 \begin{align}
 \boldsymbol{\dot\varepsilon} = 
 \frac{1}{2\eta_{\text{eff}}}\boldsymbol{\tau} +
-\frac{1}{2G} \frac{D\boldsymbol{\tau}}{Dt}  + \dot\lambda\frac{\partial Q}{\partial \boldsymbol{\tau_{II}}}
+\frac{1}{2G} \frac{D\boldsymbol{\tau}}{Dt}  + \dot\lambda\frac{\partial Q}{\partial \boldsymbol{\tau}_{II}}
 \end{align}
 $$
 
-where $\eta_{\text{eff}}$ is the effective viscosity, $G$ is the elastic shear modulus, $\dot\lambda$ is the plastic multiplier, and $Q$ is the plastic flow potential.
+where $\eta_{\text{eff}}$ is the effective viscosity, $G$ is the elastic shear modulus, $\dot\lambda$ is the plastic multiplier, and $Q$ is the plastic flow potential is:
+$$
+\begin{align}
+Q = \boldsymbol{\tau}_{II} - p\sin{\psi}
+\end{align}
+$$
+where $\psi$ is the dilation angle.
+
+## Effective viscosity
+
+The effective material of a material defined as:
+$$
+\begin{align}
+\eta_{\text{eff}} = \frac{1}{\frac{1}{\eta^{\text{diff}}} + \frac{1}{\eta^{\text{disl}}}}
+\end{align}
+$$
 
+where $\eta^{\text{diff}}$ and $\eta^{\text{disl}}$ are the diffusion and dislocation creep viscosities, which are computed from their respective strain rate equations:
+
+$$
+\begin{align}
+\dot{ε}_{II}^{\text{diff}} = A^{\text{diff}} τ_{II}^{n^{\text{diff}}} d^{p} f_{H_2O}^{r^{\text{diff}}} \exp \left(- {{E^{\text{diff}} + PV^{\text{diff}}} \over RT} \right) \\
+\dot{ε}_{II}^{\text{disl}} = A^{\text{disl}} τ_{II}^{n^{\text{disl}}} f_{H_2O}^{r^{\text{disl}}} \exp \left(- {{E^{\text{disl}} + PV^{\text{disl}}} \over RT} \right)
+\end{align}
+$$
+
+where $A$ material specific parameter, $n$ is the stress powerlaw exponent, $p$ is the negative defined grain size exponent, $f$ is the water fugacity, $r$ is the water fugacity exponent, $E$ is the activation energy, $PV$ is the activation volume, and $R$ is the universal gas constant.
 ## Elastic stress
 
 ### Method (1): Jaumann derivative
@@ -46,24 +71,6 @@ $$
 
 ### Method (2): Euler-Rodrigues rotation
 
-```julia
-# from Anton's talk
-@inline Base.@propagate_inbounds function rotate_elastic_stress3D(ωi, τ, dt)
-    # vorticity
-    ω = √(sum(x^2 for x in ωi))
-    # unit rotation axis
-    n = SVector{3,  Float64}(inv(ω) * ωi[i] for i in 1:3)
-    # integrate rotation angle
-    θ = dt * 0.5 * ω
-    # Euler Rodrigues rotation matrix
-    R = rodrigues_euler(θ, n)
-    # rotate tensor
-    τij = voigt2tensor(τ)
-    τij_rot = R * (τij * R')
-    tensor2voigt(τij_rot)
-end
-```
-
 1. Compute unit rotation axis
 $$
 \begin{align}
@@ -87,11 +94,11 @@ $$
 \begin{align}
     \boldsymbol{c} = \boldsymbol{n}\sin{\theta} \\
     \boldsymbol{R_1} = 
-    \left[
-        c0   -c_3   c_2 \\
-        c_3    c0  -c_1 \\
-       -c_2   c_1    c0 \\
-    \right] \\
+    \begin{bmatrix}
+        c0  & -c_3 &  c_2 \\
+        c_3 &  c0  & -c_1 \\
+       -c_2 &  c_1 &  c0 \\
+    \end{bmatrix} \\
     \boldsymbol{R_2} = (1-\cos{\theta}) \boldsymbol{n} \boldsymbol{n}^T \\
     \boldsymbol{R}=\boldsymbol{R_1}+\boldsymbol{R_2}
 \end{align}
@@ -107,20 +114,89 @@ $$
 
 ## Plastic formulation
 
-JustRelax.jl implements the regularised plasticity model from (refs). In this formulation, the yield function is given by
+JustRelax.jl implements the regularised plasticity model from [Thibault et al 2021](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2021GC009675). In this formulation, the yield function is given by
 
 $$
 \begin{align}
-F = \tau_y - \left( P \sin{\phi} + C \cos{\phi} + \eta_{\text{reg}} + \dot\lambda \right) \leq 0
+F = \tau_y - \left( P \sin{\phi} + C \cos{\phi} + \dot\lambda \eta_{\text{reg}}\right) \leq 0
 \end{align}
 $$
 
 where $\eta_{\text{reg}}$ is a regularization term and 
 
 $$
 \begin{align}
-\dot\lambda = \frac{F}{x}
+\dot\lambda = \frac{F^{\text{trial}}}{\eta + \eta_{\text{reg}} + K \Delta t \sin{\psi}\sin{\phi}}
 \end{align}
 $$
 
-Note that since we are using an iterative method to solve the APT Stokes equation, the non-linearities are dealt by the iterative scheme. Othwersie, one would need to solve a non-linear problem to compute $\dot\lambda$, which typically requires to compute $\frac{\partial F}{\partial \dot\lambda}$
+Note that since we are using an iterative method to solve the APT Stokes equation, the non-linearities are dealt by the iterative scheme. Othwersie, one would need to solve a non-linear problem to compute $\dot\lambda$, which requires to compute $\frac{\partial F}{\partial \dot\lambda}$
+
+# Selecting the constitutive model in JustRelax.jl
+
+All the local calculations corresponding to the effective rheology are implemented in GeoParams.jl. The composite rheology is implemented using the `CompositeRheology` object. An example of how to set up the a visco-elasto-viscoplastic rheology is shown below:
+
+```julia
+# elasticity
+el    = ConstantElasticity(; 
+    G = 40e9, # shear modulus [Pa]
+    ν = 0.45, # Poisson coeficient 
+)
+# Olivine dislocation law from Hirth and Kohlstedt 2003
+disl_wet_olivine  = SetDislocationCreep(Dislocation.wet_olivine1_Hirth_2003)
+# Olivine diffusion law from Hirth and Kohlstedt 2003
+diff_wet_olivine  = SetDiffusionCreep(Diffusion.wet_olivine_Hirth_2003)
+# plasticity
+ϕ     = 30   # friction angle
+Ψ     = 5    # dilation angle
+C     = 1e6  # cohesion [Pa]
+η_reg = 1e16 # viscosity regularization term [Pa s]
+pl    = DruckerPrager_regularised(; C = C, ϕ = ϕ, η_vp=η_reg, Ψ=Ψ)
+# composite rheology 
+rheology = CompositeRheology(
+    (el, disl_wet_olivine, diff_wet_olivine, pl)
+)
+```
+
+`rheology` then needs to be passed into a `MaterialParams` object with the help of `SetMaterialParams`:
+
+```julia
+rheology = (
+    SetMaterialParams(;
+        Phase             = 1,
+        CompositeRheology = rheology,
+        # other material properties here
+    ),
+)
+```
+
+### Multiple rheology phases
+
+It is common in geodynamic models to have multiple rheology phases. In this case, we just need to build a tuple containing every single material phase properties:
+```julia
+rheology = (
+    SetMaterialParams(;
+        Phase             = 1,
+        CompositeRheology = rheology_one,
+        # other material properties here
+    ),
+        SetMaterialParams(;
+        Phase             = 2,
+        CompositeRheology = rheology_two,
+        # other material properties here
+    ),
+)
+```
+
+### Computing the effective viscosity
+
+The effective viscosity is computed internally during the Stokes solver. However, it can also be computed externally with the `compute_viscosity!` function as follows:
+
+```julia
+# Rheology
+args             = (T=T, P=P, dt = dt) # or (T=thermal.Tc, P=stokes.P, dt=dt)
+viscosity_cutoff = (1e18, 1e23)
+compute_viscosity!(stokes, phase_ratios, args, rheology, viscosity_cutoff)
+```
+
+where `T` and `P` are the temperature and pressure fields defined at the **cell centers**, `dt` is the time step, and `phase_ratios` is an object containing the phase ratios corresponding to each material phase, of each cell.

docs/src/man/equations_APT.md

@@ -1,15 +1,15 @@
-The pseudo-transient method consists in augmenting the right-hand-side of the target PDE with a pseudo-time derivative (where $\psi$ is the pseudo-time) of the primary variables. We then solve the resulting system of equations with an iterative method. The pseudo-time derivative is then gradually reduced, until the original PDE is solved and the changes in the primary variables are below a preset tolerance.
+The Accelerated Pseudo-Transient (APT) method consists in augmenting the right-hand-side of the target PDE with a pseudo-time derivative (where $\psi$ is the pseudo-time) of the primary variables. We then solve the resulting system of equations with an iterative method. The pseudo-time derivative is then gradually reduced, until the original PDE is solved and the changes in the primary variables are below a preset tolerance.
 
 ## Heat diffusion
-The pseudo-transient heat-diffusion equation is:
+The APT heat-diffusion equation is:
 
 $$
 \begin{align}
 \widetilde{\rho}\frac{\partial T}{\partial \psi} + \rho C_p \frac{\partial T}{\partial t} = \nabla \cdot (\kappa\nabla T) = -\nabla q
 \end{align}
 $$
 
-We use a second order pseudo-transient scheme were continuation is also done on the flux, so that:
+We use a second order APT scheme were continuation is also done on the flux, so that:
 
 $$
 \begin{align}
@@ -19,7 +19,7 @@ $$
 
 ## Stokes equations
 
-For example, the pseudo-transient formulation of the Stokes equations yields:
+For example, the APT formulation of the Stokes equations yields:
 
 $$
 \begin{align}
@@ -34,7 +34,7 @@ $$
 $$
 
 ## Constitutive equations
-A pseudo-transient continuation is also done on the constitutive law:
+A APT continuation is also done on the constitutive law:
 
 $$
 \begin{align}
@@ -60,7 +60,7 @@ $$
 
 where the P-wave $\widetilde{V}=V_p$ is the characteristic velocity scale for Stokes, and $Re$ is the Reynolds number.
 
-### Physical parameters
+<!-- ### Physical parameters
 
 | Symbol                           | Parameter              |
 | :------------------------------- | :--------------------: |
@@ -77,7 +77,7 @@ where the P-wave $\widetilde{V}=V_p$ is the characteristic velocity scale for St
 | $G$                              | Shear modulus          |
 | $\alpha$                         | Thermal expansivity    |
 | $C_p$                            | Heat capacity          |
-| $\kappa$                         | Heat conductivity      |
+| $\kappa$                         | Heat conductivity      | -->
 
 ### Pseudo-transient parameters
 

docs/src/man/equations_basic.md

@@ -1,3 +1,5 @@
+# Governing equations
+
 ## Stokes equations
 
 Deformation of compressible viscous flow is desribed by the equations of conservation of momentum and mass:
@@ -16,6 +18,8 @@ $$
 
 where $\boldsymbol{\tau}$ is the deviatoric stress tensor, $p$ is pressure, $f$ is the external forces vector, $\boldsymbol{v}$ is the velocity vector, $\beta$ is the compressibility coefficient, $\alpha$ is the thermal expansivity coefficient and $T$ is temperature.
 
+## Constitutive equation
+
 To close the previous systems of equations, we further need the constitutive relationship between stress and deformationl. In its simplest linear form this is:
 
 $$
@@ -25,7 +29,6 @@ $$
 $$
 
 where $\eta$ is the shear viscosity and  $\boldsymbol{\dot\varepsilon}$ is the deviatoric strain tensor.
-
 ## Heat diffusion
 The pseudo-transient heat-diffusion equation is:
 
@@ -36,9 +39,3 @@ $$
 $$
 
 where $\rho$ is density, $C_p$ is specific heat capacity, $\kappa$ is thermal conductivity, $T$ is temperature $\boldsymbol{\tau}:\boldsymbol{\dot\varepsilon}$ is the energy dissipated by viscous deformation (shear heating), $\alpha T(\boldsymbol{v} \cdot \nabla P)$ is adiabatic heating, and $H$ is the sum any other source term, such as radiogenic heat production.
-
-# Discretization of the governing equations
-
-We discretize both the Stokes and heat diffussion equations using a Finite Differences approach on a staggered grid (ref Taras book here), as sketched below:
-
-![Staggered Velocity Grid](../assets/StaggeredVelocity.png)

docs/src/man/equations_discretization.md

@@ -2,8 +2,6 @@
 
 We discretize both the Stokes and heat diffussion equations using a Finite Differences approach on a staggered grid (ref Taras book here).
 
-
-
 ## Heat diffusion
 The heat diffusion equation is discretized on a staggered grid as sketched below:
 
@@ -39,8 +37,7 @@ The equations of conservation of mass and momentum are discretized on a staggere
 ![](../assets/stokes_stag2D.png)
 
 
-where dotted lines represent the velocity ghost nodes. The APT Stokes equations are then discretized as follows:
-
+where dotted lines represent the velocity ghost nodes. 
 <!-- $$
 \begin{align}
 \widetilde{\rho}\frac{\boldsymbol{u}}{\Delta\psi} + \nabla\cdot\boldsymbol{\tau} - 
@@ -49,6 +46,8 @@ where dotted lines represent the velocity ghost nodes. The APT Stokes equations
 \end{align}
 $$ -->
 
+### Conservation of momentum
+
 $$
 \begin{align}
 \widetilde{\rho}\frac{u^{n+1}_x - u^n_x}{\Delta\psi} + \nabla\cdot\boldsymbol{\tau} - 
@@ -72,4 +71,9 @@ $$
 \frac{1}{\Delta t} \left( \beta (p^{t+\Delta t} - p^{t}) + \alpha (T^{t+\Delta t} - T^{t}) \right) \\
 
 \end{align}
-$$
+$$
+
+### Conservation of mass
+
+
+### Constitutive equation