title: "Simulating microstructures with gmGeostats"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{my-vignette}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
echo = TRUE,
comment = "#>",
fig.width = 7,
fig.height = 6
)
# ToDo: anistropoy: what means the ratio exactly?
```
```{r setup, echo=FALSE, message=FALSE}
library(gmGeostats)
library(ggplot2)
library(data.table) # only needed, because Solveig is faster with melting in data.table than with dplyr and Co. can be replaced by dplyr, if someone changes the code from data.table to dplyr later.
```
# Part 1: Granite/Gneiss
Granites or Gneiss consist mainly of Quartz, Feldspars and micas, e.g. Biotite and/or Muscovite.
Due to the brittle characteristics these types of rocks have often fractures or joint systems.
With `gmGeostats` mineral composition, anisotrpy of certain mineral phases and features like joints and fractures can be simulated, e.g. if simulations of crystalline rocks are required for testing code or work flows.
## General idea - Turning bands algorithm
here comes some nice background from Raimon
## Example for simulation
The simulation consists of five steps:
1. decide for the number of phases to be simulated. Fracture and joint systems count as a "phase", e.g. if there is one joint system, this would be one phase, if there are two joint systems this would be simulated by two additional phases.
2. define the characteristics of each phase, i.e. the crystal elongation (or void shape) and orientation/foliation
3. generate the variograms by using a function to define for each phase the model of the phase-variogram, the nugget, the sill, the range and the anistropy with respect to the range.
4. use the algorithm of Turning Bands to calculate the probability for each phase to occur.
5. Visualize the result.
In the following example, we will show the implementation for 2D.
Because the code expects 3D the third dimension is always "muted" by using zeros.
### Number of phases
For the example let's assume five phases, Quartz (qu), K-Feldspar (fsp), Biotite (bi), a fracture system (frac) and a joints system (joints).
Generate a template matrix with only zeros, which has the dimensions `n x n`:
In a later step when the five variograms are stacked, these five matrices define the sill for each phase.
Therefore, on the diagonal we replace the `0` with a `1` for the respective phase:
```{r}
s1 = template
s1[1,1] <-1
s2 = template
s2[2,2] <-1
s3 = template
s3[3,3] <-1
s4 = template
s4[4,4] <-1
s5 = template
s5[5,5] <-1
```
The example for Quartz looks like this:
```{r}
s1
```
### Characteristics of phases
The elongation and orientation of the elongation for each phase can be simulated by the function `anis2D.par2A`.
The function parameter `ratio` provides the ratio of ??? (length to width?).
E.g. `ratio = 1` is a sphere (isotrop), while `ratio = 2` would describe a mineral which is ??? (double as long as thick?).
The function parameter `angle` gives the orientation (???, does not really, the orientation seems to be somewhat different) of this anisotropy, counted counterclockwise from the east (E).
E.g. `angle = 45` with `ratio = 2` would result into North-South elongated crystals.
For Quartz, let's assume we have mostly roundish/quadratic crystals and only a minor preferred orientation of the crystals, so that they are slightly oriented along NNE-SSW:
```{r p-qu}
p1_qu = anis2D.par2A(ratio = .9, angle = 45)
```
The Feldspars have a more pronounced elongated shape, also with preferred orientation into the same direction as the Quartz grains:
```{r p-fsp}
p2_fsp = anis2D.par2A(ratio = .6, angle = 45)
```
If the anistropy of the mineral phases had been generated by the function `anis2D.par2A`, then we need to calculate the cross product of the matrices to gain the anistropy matrices `A`:
```{r}
A1_qu = tcrossprod(p1_qu)
A2_fsp = tcrossprod(p2_fsp)
```
This results into
```{r}
A1_qu
```
for Quartz and into
```{r}
A2_fsp
```
for Feldspar.
There is also the option to directly generate the anisotropy matrix.
In this case, the matrix can be written manually:
The Biotite occurs as thin plates and has a slightly different preferred orientation then the two main minerals:
For the fracture and joint systems the anisotropy is very pronounced, because they are very thin but passing through longer distances through the rock.
Accordingly, the anisotropy matrix should be chosen similar to this:
If single the phases should be changed then only certain steps have to be retraced.
For example, in the simulation above, the geologist would expect, that the fractures rather follow the biotite layers, instead of being orthogonal to the joints.
Changing the orientation of the fractures:
```{r}
p4_frac = anis2D.par2A(ratio = .01, angle = 95)
A4_frac = tcrossprod(p4_frac)
A4_frac
```
Then recalculate the variogram for the fractures and stack all variograms again:
