diff --git a/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/3D_shear_box_vel_particle.pvsm b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/3D_shear_box_vel_particle.pvsm
new file mode 100644
index 00000000000..16fee894f0f
--- /dev/null
+++ b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/3D_shear_box_vel_particle.pvsm
@@ -0,0 +1,15593 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 
+
+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
+
+
+Initial randomized CPO at time zero.
+```
+
+```{figure-md} fig:olivine_fabrics
+
+
+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.
diff --git a/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/initial_random_cpo.png b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/initial_random_cpo.png
new file mode 100644
index 00000000000..76861c61dcd
Binary files /dev/null and b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/initial_random_cpo.png differ
diff --git a/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/model_setup_3D_box.png b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/model_setup_3D_box.png
new file mode 100644
index 00000000000..980c4506fa1
Binary files /dev/null and b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/model_setup_3D_box.png differ
diff --git a/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/prescribe_stokes.part.prm b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/doc/prescribe_stokes.part.prm
new file mode 100644
index 00000000000..649975b3511
--- /dev/null
+++ b/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
diff --git a/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/olivineA.prm b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/olivineA.prm
new file mode 100644
index 00000000000..d79c8699d86
--- /dev/null
+++ b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/olivineA.prm
@@ -0,0 +1,170 @@
+set CFL number = 0.01 ##Annotation CFL number
+set Timing output frequency = 1000
+set Dimension = 3
+set Pressure normalization = surface
+set Surface pressure = 0
+set Nonlinear solver scheme = single Advection, no Stokes
+set End time = 3e5
+set Use years in output instead of seconds = false
+set Output directory = out_drex_olivineA
+
+subsection Compositional fields
+ set Number of fields = 1
+ set Names of fields = water
+ set Compositional field methods = static
+end
+
+subsection Initial composition model
+ set Model name = function
+
+ subsection Function
+ set Variable names = x,y,z
+ set Function expression = 0
+ end
+end
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 0
+ end
+end
+
+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
+
+subsection Initial temperature model
+ set Model name = function
+
+ subsection Function
+ set Function expression = 1 ## Annotation: Temperature function
+ end
+end
+
+subsection Boundary temperature model
+ set List of model names = box
+ set Fixed temperature boundary indicators = bottom, top, left, right, front, back
+
+ subsection Box
+ set Bottom temperature = 1
+ set Left temperature = 1
+ set Right temperature = 1
+ set Top temperature = 1
+ set Front temperature = 1
+ set Back temperature = 1
+ end
+end
+
+subsection Boundary velocity model
+ subsection Function
+ set Variable names = x,y,z,t
+ set Function expression = z/1e5;0;0 ## Annotation set velocity condition
+ end
+end
+
+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
+
+subsection Material model
+ set Model name = simple
+
+ subsection Simple model
+ set Reference density = 1
+ set Reference specific heat = 1
+ set Reference temperature = 1
+ set Thermal conductivity = 1
+ set Thermal expansion coefficient = 1e-10
+ set Viscosity = 1
+ end
+end
+
+subsection Mesh refinement
+ set Initial global refinement = 1
+ set Initial adaptive refinement = 0
+ set Time steps between mesh refinement = 10000
+end
+
+subsection Postprocess
+ set List of postprocessors = velocity statistics, visualization, particles, crystal preferred orientation
+
+ subsection Visualization
+ set Time between graphical output = 1e5
+ set List of output variables = stress, viscosity, strain rate
+ end
+
+ subsection Crystal Preferred Orientation
+ set Time between data output = 1e5
+ 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 #, mineral 1: volume fraction, mineral 1: Euler angles
+ set Write out draw volume weighted cpo data = mineral 0: volume fraction, mineral 0: Euler angles, mineral 1: volume fraction, mineral 1: Euler angles
+ end
+
+ subsection Particles
+ set Time between data output = 1e4
+ set Data output format = vtu
+ set Exclude output properties = volume fraction, rotation_matrix
+ end
+end
+
+subsection Solver parameters
+ set Temperature solver tolerance = 1e-10
+end
+
+subsection Particles
+ set List of particle properties = integrated strain, crystal preferred orientation, cpo elastic tensor, cpo bingham average, integrated strain invariant, velocity, pT path
+
+ #set Load balancing strategy = add particles
+ #set Minimum particles per cell = 128
+ #set Maximum particles per cell = 512
+ set Particle generator name = ascii file
+
+ subsection Generator
+ subsection Ascii file
+ set Data directory = /Users/tian_bc/Documents/postdoc/postdoc_meetings_workshops/2024-ASPECT-hack/aspect-output/drex_olivine_branch_CPO_based_on_yijun/main240613/
+ set Data file name = particle_one.dat
+ end
+ end
+
+ subsection Crystal Preferred Orientation
+ 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
+ set Minerals = Olivine: A-fabric, Enstatite ## Annotation Minerals
+ set Volume fractions minerals = 0.7,0.3 ## Annotation Volume fractions minerals
+ end
+
+ 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
diff --git a/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/particle_one.dat b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/particle_one.dat
new file mode 100644
index 00000000000..d5e8e721f44
--- /dev/null
+++ b/cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/particle_one.dat
@@ -0,0 +1 @@
+0 0 0
\ No newline at end of file