SimulatingMicrostructures.Rmd

---
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`:
```{r phases}
n_phases = 5

template = matrix(0, nrow = n_phases, ncol = n_phases)
phases = c("qu", "fsp", "bi", "frac", "joints")
colnames(template) <- rownames(template) <- phases
```
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:
```{r}
A3_bi = matrix(data = c(1, .6, 0, .6, 1, 0, 0, 0, 1), ncol = 3)
A3_bi
```

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:
```{r}
A4_frac = matrix(data = c(1, -.999, 0, -.999, 1, 0, 0, 0, 1), ncol = 3)
A4_frac
```

And let's assume, that the joint system is orthogonal to the fracture system:
```{r}
A5_joints = matrix(data = c(1, .999, 0, .999, 1, 0, 0, 0, 1), ncol = 3)
A5_joints
```

### Generate variograms

Fore each phase, variogram functions will predict the probability to occur.
The function `setCgram` generates these variograms.
It has four parameters: 

- The `type` provides the variogram model, e.g. Gaussian or Exponential.
Gaussian models provide smooth transitions to other minerals, and may be used if bands of minerals should be modelled. 
Exponential models provide rather abrupt changes between mineral phases and should be preferred model for sharp edges.
In this example, Quarts and Feldspar are modelled by Gaussion types, and Biotite and the fractures by Exponential type.
- The `nugget` gives the nugget effect, but in the case of microstructures it can be omitted.
- For the `sill`, the partial sill of the correlation function, the microstructure require simply the matrices with one `1` in the diagonal.
We produced these matrices in the first part, where the number of phases was defined.
- The `anisRages` use the information from the anistropy matrices generated in the previous part.
This anisotropy range has to be downscaled a bit, therefore a factor is applied prior to the the anisotropy matrices:

```{r}
v1_qu = setCgram(type = vg.Gau, sill = s1, anisRanges = 0.015*A1_qu)
v2_fsp = setCgram(type = vg.Gau, sill = s2, anisRanges = 0.015*A2_fsp)
v3_bi = setCgram(type = vg.Exp, sill = s3, anisRanges = 0.015*A3_bi)
v4_frac = setCgram(type = vg.Exp, sill = s4, anisRanges = 0.05*A4_frac)
v5_joints = setCgram(type = vg.Exp, sill = s5, anisRanges = 0.05*A5_joints)
```

Plotting one of these variograms shows this:
```{r}
plot(v3_bi)
```

The factor applied to the anisotropy matrix is can be changed according to target patterns, which should be achieved by the simulation.

The variograms have now to be stacked:
```{r}
vm = v5_joints + v4_frac + v3_bi  + v2_fsp + v1_qu
```
Please note, that the order of the stacking is important.
The first phase will have least influence, the last phase will have the main influence.

The visualisation of the stacked variograms looks as follows:

```{r}
plot(vm)
```

### Apply Turning Bands algorithm

In a first step, a grid of nodes has to be generated, to determine the number and extend of the pixels of the finalized "image".
This can be done by `expand.grid`:
```{r}
# make grid
resol=200
Xgrid = as.matrix(expand.grid(x = seq(0,1, length.out = resol),
                              y = seq(0,1, length.out = resol),
                              z = 0))
```

Pleas note: `Xgrid` needs to have the class `matrix`, and the third dimension is muted by using only zeros for the z-direction.

The data are generated by:
```{r, message=FALSE}
set.seed(122487) # 122487
rock1 = gsi.TurningBands(X = Xgrid, 
                         vgram = vm, 
                         nBands = 1000)
colnames(rock1)[-(1:3)] = phases # apply the names
```

This gives back a large matrix which looks like this:
```{r}
rock = data.frame(rock1)
head(rock)
```

For each 'row' in the resulting image, i.e. keep a y-value constant, and look along the x-axis, the five phases have a certain probability to appear.
For example for y = 0.10025 :
```{r, echo=F}
rock_melt = melt(data.table(rock), id.vars = c("x", "y", "z"), 
                            variable.name = "phase", 
                            value.name = "relation")
rock_melt[, phase := factor(phase, levels = phases)]
ggplot(rock_melt[y > 0.1 & y < 0.101], aes(x = x, y = relation, color = phase)) + 
  theme_bw() +
  #geom_hline(yintercept = 0, color = "red") +
  geom_line() +
  scale_color_brewer(type = "qual", palette = "Set1")
```

For each pixel, i.e. for each x/y combination respectively each row in `rock`, the phase with maximum value is chosen:
```{r}
rocksims = apply(rock1[, -(1:3)],1, function(x) which.max(x))
dim(rocksims) = c(resol, resol)
```

### Visualization

The resulting image is:
```{r, echo=F}
# somehow ggplot doesnt work...
# phase_info = data.table(n_phase = c(1:5), phase = factor(phases, levels = phases))
# rocksims_melt = melt(data.table(data.frame(rocksims), keep.rownames = T), id.vars = "rn", value.name = "n_phase")
# rocksims_melt = rocksims_melt[phase_info, on = "n_phase"]
# # plot
# ggplot(rocksims_melt, aes(y = rn, x = variable, fill = phase)) + 
#   geom_tile() + 
#   coord_equal() + 
#   scale_fill_brewer(type = "qual", palette = "Set1") + 
#   theme(axis.text = element_blank(), axis.title = element_blank(), axis.ticks = element_blank())
```
```{r}
image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)
```




If one of the phases should be less pronounced, e.g. here in the example the joints are very frequent, then the phases can be regulated by weights.
The weights simply increase or diminish the values of the `rock` matrix for the phase.
For example:
```{r}
ph_weights = c(0,0,0,0,-1) # diminish only the joints
rocksims = apply(rock1[, -(1:3)],1, function(x) which.max(x + ph_weights))
dim(rocksims) = c(resol, resol)
```

Then the the line is pushed down:
```{r, echo=F}
rock_melt = melt(data.table(rock)[, joints := joints -1], id.vars = c("x", "y", "z"), 
                            variable.name = "phase", 
                            value.name = "relation")
rock_melt[, phase := factor(phase, levels = phases)]
ggplot(rock_melt[y > 0.1 & y < 0.101], aes(x = x, y = relation, color = phase)) + 
  theme_bw() +
  #geom_hline(yintercept = 0, color = "red") +
  geom_line() +
  scale_color_brewer(type = "qual", palette = "Set1")
```

And then the pattern looks like:

```{r}
image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)
```

### Example with more realistic fractures

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:
```{r}
v4_frac = setCgram(type = vg.Exp, sill = s4, anisRanges = 0.05*A4_frac)
vm = v5_joints + v4_frac + v3_bi  + v2_fsp + v1_qu
```

Plot:
```{r}
plot(vm)
```

Applying the algorithm of the turning bands:
```{r}
set.seed(384952) 
rock2 = gsi.TurningBands(X = Xgrid, vgram = vm, nBands = 1000)
colnames(rock2)[-(1:3)] = phases
```

Get the phases per pixel, including the weights for the joints:
```{r}
ph_weights = c(0,0,0,0,-1) # diminish only the joints
rocksims = apply(rock2[, -(1:3)],1, function(x) which.max(x + ph_weights))
dim(rocksims) = c(resol, resol)
```

and then visualize:
```{r}
image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)
```

In this plot, the fracture seem to be too big.
Again, this can be easily resolved by downweighing the fracturesa (and joints) and increase a bit the Quartz:
```{r}
ph_weights = c(.1,0,0,-1.2,-1) # diminish joints and fracture, increase quartz
rocksims = apply(rock2[, -(1:3)],1, function(x) which.max(x + ph_weights))
dim(rocksims) = c(resol, resol)
```

```{r}
image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)
```