[WIP] Cookbook crystal preferred orientation calculation for olivine based on Fraters and Billen 2021

cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/cpo.part.prm

@@ -0,0 +1,87 @@
+#specify posstprocessor objects that will be run at the end of each time step
+#The statistics is important for cpo plotter for reading the data
+#The particles and crystal preferred orientation are required for CPO calculations
+subsection Postprocess
+  set List of postprocessors = velocity statistics, visualization, particles, crystal preferred orientation
+
+  #specify data output for visualizing in paraview or visit
+  subsection Visualization
+    set Time between graphical output = 1e5
+    set List of output variables = stress, viscosity, strain rate
+  end
+
+  #specify CPO data output
+  subsection Crystal Preferred Orientation
+    set Time between data output = 1e5
+    set Write in background thread = true
+    set Compress cpo data files = false # so results are human readable
+    set Write out raw cpo data = mineral 0: volume fraction, mineral 0: Euler angles #, mineral 1: volume fraction, mineral 1: Euler angles #specify what to output, here will only output volume fraction and Euler angles for mineral 0 (olivine)
+    set Write out draw volume weighted cpo data = mineral 0: volume fraction, mineral 0: Euler angles, mineral 1: volume fraction, mineral 1: Euler angles
+
+    #specify output parameters for draw volume weighted cpo, the pole figures are plotted with these data
+  end
+
+  # Here is for data results read by paraview and visit
+  # Exclude output properties so as to allow fast reading by paraview
+  # because if number of particles is large, paraview will crash
+  subsection Particles
+    set Time between data output = 1e4
+    set Data output format       = vtu
+    set Exclude output properties = volume fraction, rotation_matrix
+  end
+end
+
+# subsection Particles was part of subsection Postprocess but from V3.0 it standalone
+subsection Particles
+  set List of particle properties = integrated strain, crystal preferred orientation, cpo elastic tensor, cpo bingham average, integrated strain invariant, velocity, pT path #what particle properties are calculated
+
+  # When specifying locations for particles using subsection Generator below, below 3 lines are commented
+  #set Load balancing strategy = add particles
+  #set Minimum particles per cell = 128
+  #set Maximum particles per cell = 512
+  # use a ascii file particle_one.dat to specify initial locations of particles
+  set Particle generator name = ascii file
+
+  #
+  subsection Generator
+    subsection Ascii file
+      # you may need to change the below path to the absolute path in your local computer
+      # the particle_one.dat is located in the cookbook folder
+      set Data directory = /aspect/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/
+      set Data file name = particle_one.dat
+    end
+  end
+
+  # parameters for DRex
+  subsection Crystal Preferred Orientation
+    # a random seed to 'randomize' initial Euler angles for the grains
+    set Random number seed = 301
+    set Number of grains per particle = 500 ## Annotation grain count
+    set Property advection method = Backward Euler ## Annotation Property advection method
+    set Property advection tolerance = 1e-15
+    set CPO derivatives algorithm = D-Rex 2004
+
+    subsection Initial grains
+      # prescribe types of mineral
+      # The A,B,C,D,E types of olivine fabric can be set here by changing
+      # A-fabric to B-fabric for example
+      set Minerals = Olivine: A-fabric, Enstatite ## Annotation Minerals
+
+      # The two minerals are not interacting with each other.
+      set Volume fractions minerals = 0.7,0.3 ## Annotation Volume fractions minerals
+    end
+
+    # DRex parameters following Kaminski et al., 2004
+    subsection D-Rex 2004
+      set Mobility = 125 ## 125+-75 Annotation LPO Mobility
+      set Stress exponents = 3.5
+      set Exponents p = 1.5
+      set Nucleation efficiency = 5 ## Annotation Nucleation efficiency
+      set Threshold GBS = 0.3 ## Annotation TGBS
+    end
+  end
+
+  subsection CPO Bingham Average
+    set Random number seed = 200
+  end
+end

cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/cpo_pole_figure_config.toml

@@ -0,0 +1,15 @@
+base_dir = "/Users/tian_bc/Documents/postdoc/postdoc_meetings_workshops/2024-ASPECT-hack/aspect-output/drex_olivine_branch_CPO_based_on_yijun/main240613/"
+experiment_dirs = ["out_drex_olivineE/"]
+compressed = false
+
+[pole_figures]
+    elastisity_header = false #not yet implemented
+    small_figure = true
+    no_description_text = false #false
+    color_scale = "Gray"
+    times = [1.0,1e5,3e5]
+    particle_ids = [0]
+    axes = ["AAxis","BAxis","CAxis"]
+    minerals = ["Olivine","Enstatite"]
+    grain_data_file_prefix = "weighted_CPO"
+    particle_data_dir = "particles_cpo"

cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/crystal_preferred_orientation_olivine_fraters_billen_2021.md

@@ -0,0 +1,129 @@
+# [WIP] Olivine Fabric Developments Under Simple Shear
+
+*This section was contributed by Xiaochuan Tian, with helps from Yijun Wang and Menno Fraters. It is
+based on a section in {cite:t}`fraters_billen_2021_cpo`
+by Menno Fraters and Magali Billen published in 2021.*
+
+This cookbook explains how to set up a numerical experiment for fabric
+developments of a single olivine particle under simple shear macroscopic strain.
+It uses ASPECT's
+crystal preferred orientation (CPO) implementation, which is described in detail in
+{cite:t}`fraters_billen_2021_cpo` and the paper is open accessed at
+[here](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2021GC009846).
+The fabric calculation is based on DRex (Kaminski et al., 2004) at
+[here](https://academic.oup.com/gji/article/158/2/744/2013756).
+This notebook describe how to use ASPECT CPO implementation to calculate the major
+olivine fabrics observed in the lab experiments. It focuses on reproducing Fig. 3
+of {cite}`fraters_billen_2021_cpo`.
+
+## Motivation
+The Earth's plates can move relative to each other because the underlying mantle
+can flow/(be deformed) to accommodate such motions.
+This flowing motion reflects dynamics inside the planet but cannot be observed directly
+as the depth of interest is beyond human's drilling ability. However, indirect observation
+from seismic anisotropy provides information about directionality of the speed of seismic
+waves that propagate inside the Earth. If such directionality is controlled by
+the mantle's motion/deformation, we can then use seismic anisotropy observations
+to infer mantle's deformation.
+
+Olivine is a major mineral of the Earth's upper mantle where continuous deformation
+takes place. People care about olivine fabric developments under simple shear as it may
+provide a bridge that links seismic anisotropy observations to the
+flows in the upper mantle that accommodate plate motions. Indeed, high temperature and
+pressure lab experiments using Griggs apparatus to investigate olivine fabrics developments
+under simple shear do find systematic types of fabrics based on the stress and water content conditions
+(see Karato et al., 2008 for details). Here is a
+[link](https://doi.org/10.1146/annurev.earth.36.031207.124120) to the paper.
+
+## Model setup
+
+Following {cite:t}`fraters_billen_2021_cpo`, we prescribe simple shear in a 3d Cartesian box/cube
+with dimensions of $1 \times 1 \times 1 $ $[m^3]$. The shear strain rate
+$\dot{\epsilon}_{xz} = -5\times 10^{-6} [s^{-1}]$. The olivine particle is placed
+right at the center of the cubic box so it stays stationary. The DRex implementation
+keep tracks of rotations of crystal grains within the particle under macroscopic deformation.
+
+```{figure-md} fig:model_setup_3D_box
+<img src="model_setup_3D_box.png" style="width:80.0%" />
+
+3D shear box with velocity vectors. The grey ball is the central olivine particle.
+```
+
+The model shows how crystal grains rotate and align under simple shear with pole figures.
+
+## The input file
+
+:::{tip}
+Update the input file with update_prm_files.sh inside
+/aspect/contrib/utilities if the model doesn't run.
+See Pull Request #5873
+[conversations](https://github.com/geodynamics/aspect/pull/5873#issuecomment-2167145176) for detail.
+:::
+
+Instead of solving Stokes equation with initial and boundary conditions, we instead prescribe
+a constant simple shear strain rate field by setting Nonlinear solver scheme to
+"single Advecgion, no Stokes" and
+prescribing the Stokes solution with a function. In this case, what we only care is the
+prescribed strain rate that is exerted onto the olivine particle. Parameters useful for the Stokes
+equation is not important in this case.
+
+```{literalinclude} prescribe_stokes.part.prm
+```
+
+The second part of the input file specify what to output and what minerals to keep track of
+CPO development. It also specifies the parameters for the DRex algorithm.
+
+```{literalinclude} cpo.part.prm
+```
+
+The complete input file is located at
+[/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/olivineA.prm](https://www.github.com/geodynamics/aspect/blob/main/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/olivineA.prm).
+
+## Plotting pole figure
+To plot pole figures from
+
+## Model results
+```{figure-md} fig:initCPO
+<img src="initial_random_cpo.png" style="width:50.0%" />
+
+Initial randomized CPO at time zero.
+```
+
+```{figure-md} fig:olivine_fabrics
+<img src="Olivine_fabrics.svg" style="width:60.0%" />
+
+Olivine A-E type fabrics under simple shear with a shear strain of 1.5
+at 3e5 seconds of model time under constant shear strain rate:
+ $\dot{\epsilon}_{xz} = -5\times 10^{-6} [s^{-1}]$.
+Note that the color scales range from 0 to 5 but actual results may be larger.
+```
+
+When we look at the visualization output of this model (see also
+{numref}`fig:olivine_fabrics`), we can see
+
+
+
+## Extending the model
+
+This simple 1 $m^{3}$ cube simple shear model for tracking olivine fabric development using
+the crystal preferred orientation implementation can be extended for various settings.
+Some ideas for adapting the model setup:
+
+-   Setup high pressure and temperature lab experiments for olivine CPO.
+    After benchmarking against the lab experiment results, it could then be
+    extended beyond parameter space limited by lab experiment.
+
+-   Add CPO calculations in large scale geodynamic models like mantle plume,
+    lower crust delamination and mid-ocean ridges.
+
+-   Include anisotropic viscosity to feed this CPO back to influence the stokes solution.
+
+## 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.
+
+- 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.
+
+- Karato, S., Jung, H., Katayama, I., & Skemer, P. (2008).
+Geodynamic significance of seismic anisotropy of the upper mantle: New insights from laboratory studies. Annu. Rev. Earth Planet. Sci., 36, 59–95.

cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/prescribe_stokes.part.prm

@@ -0,0 +1,53 @@
+# we prescribe the simple shear strain rate by setting solver scheme
+# as "single Advection, no Stokes"
+# see below "subsection Prescribed Stokes solution "
+set Nonlinear solver scheme  = single Advection, no Stokes
+# The model run to 3e5 seconds, at which the total strain will
+# be 1.5 with strain rate of 5e-6 1/s.
+set End time                 = 3e5
+set Use years in output instead of seconds = false
+
+# Set zero gravity
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 0
+  end
+end
+
+# 1m by 1m by 1m cube
+# center is at (0,0,0)
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 1
+    set Y extent = 1
+    set Z extent = 1
+    set Box origin X coordinate = -0.5
+    set Box origin Y coordinate = -0.5
+    set Box origin Z coordinate = -0.5
+  end
+end
+
+# Here prescribing velocity field so that it results in
+# simple shear strain rate epsilon_xz = -5e6 1/s
+# as z ranges from -0.5 to 0.5 meter.
+subsection Prescribed Stokes solution
+  set Model name = function
+
+  subsection Velocity function
+    set Variable names = x,y,z,t
+    set Function expression = z/1e5;0;0 ## Annotation set velocity condition
+  end
+end
+
+# The default X,Y,Z repetitions is 1, so will be 2 by 2 by 2 grids without refinement
+# Here with setting "Initial global refinement" to be 1, it becomes
+4 by 4 by 4 grids.
+subsection Mesh refinement
+  set Initial global refinement                = 1
+  set Initial adaptive refinement              = 0
+  set Time steps between mesh refinement       = 10000
+end