minicorrections vignette microstructure simulation

vignettes/SimulatingMicrostructures.Rmd

@@ -42,11 +42,11 @@ 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.
+4. use the algorithm of Turning Bands to calculate the potential 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.
+Because the code expects 3D the third dimension is always "muted" by setting it as a constant.
 
 ### Number of phases
 
@@ -83,25 +83,23 @@ 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.
+The anisotropy of the elongation for each phase can be obtained by the function `anis2D.par2A`.
+The function parameter `ratio` provides the ratio of shortest to longest range, and should therefore be smaller than 1.
+E.g. `ratio = 1` is a sphere (isotrop), while `ratio = 0.5` would describe a mineral which grains are typicall double as long as thick.
+The function parameter `angle` gives the orientation of this anisotropy, measured in degrees counterclockwise from the east (E).
+E.g. `angle = 45` with `ratio = 0.5` would result into NE-SW 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:
+For Quartz, let's assume we have mostly roundish/quadratic grains and only a minor preferred orientation of the crystals, so that they are slightly oriented along NE-SW:
 ```{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:
-
+The Feldspars have a more pronounced elongated shape, also with preferred orientation along 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`:
-
+Function `anis2D.par2A` returns the Cholesky decomposition of the anisotropy matrix. So, if the anistropy of the mineral phases had been generated by this function, then we need to calculate the cross product of the matrices to obtain the anistropy matrices `A`:
 ```{r}
 A1_qu = tcrossprod(p1_qu)   
 A2_fsp = tcrossprod(p2_fsp)
@@ -142,41 +140,38 @@ A5_joints
 
 ### Generate variograms
 
-Fore each phase, variogram functions will predict the probability to occur.
+Fore each phase, variogram functions will construct the potentials that control which phase will actually 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.
+Gaussian models provide smooth transitions to other minerals, and may be used if bands of minerals in equilibrium should be modelled. 
+Exponential models provide rather abrupt changes between mineral phases and should be preferred model for fractal contacts, such as non-equilibrium contacts, dissolution/porosity, alteration, etc.
+In this example, Quarts and Feldspar are modelled by Gaussion types, and Biotite and the fractures by the Exponential type.
+- The `nugget` gives the nugget effect, but in the case of microstructures it should 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:
+This anisotropy range has to be downscaled, 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)
+v1_qu = setCgram(type = vg.Gau, sill = s1, anisRanges = 0.005*A1_qu)
+v2_fsp = setCgram(type = vg.Gau, sill = s2, anisRanges = 0.005*A2_fsp)
+v3_bi = setCgram(type = vg.Exp, sill = s3, anisRanges = 0.005*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:
+This factor controls the typical size of each phase, but in a complex, indirect way. 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:
+The variograms need now to be nested:
 ```{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.
+Please note, that the order of the stacking is not important.
 
 The visualisation of the stacked variograms looks as follows:
 
@@ -196,11 +191,11 @@ Xgrid = as.matrix(expand.grid(x = 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.
+Pleas note: `Xgrid` needs to have the class `matrix`, and the third dimension is set to constant by using only zeros for the z-direction. However, the whole algorithm here presented is directly suitable for 3D simulations: the complexity rather lies in the representation of the outcomes.
 
 The data are generated by:
 ```{r, message=FALSE}
-set.seed(122487) # 122487
+set.seed(122487)
 rock1 = gsi.TurningBands(X = Xgrid, 
                          vgram = vm, 
                          nBands = 1000)
@@ -213,7 +208,7 @@ 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 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 potential to appear.
 For example for y = 0.10025 :
 ```{r, echo=F}
 rock_melt = melt(data.table(rock), id.vars = c("x", "y", "z"), 
@@ -227,7 +222,7 @@ ggplot(rock_melt[y > 0.1 & y < 0.101], aes(x = x, y = relation, color = phase))
   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:
+For each pixel, i.e. for each x/y combination respectively each row in `rock`, the phase with maximum potential is chosen:
 ```{r}
 rocksims = apply(rock1[, -(1:3)],1, function(x) which.max(x))
 dim(rocksims) = c(resol, resol)
@@ -256,11 +251,11 @@ 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.
+If one of the phases should be less pronounced, e.g. here in the example the joints are very frequent, then the potentisal can be regulated by weights.
+The weights simply increase or diminish the potential of each phase.
 For example:
 ```{r}
-ph_weights = c(0,0,0,0,-1) # diminish only the joints
+ph_weights = c(1,1,0,0,-1) # inurease qu and fsp, diminish joints
 rocksims = apply(rock1[, -(1:3)],1, function(x) which.max(x + ph_weights))
 dim(rocksims) = c(resol, resol)
 ```
@@ -278,7 +273,7 @@ ggplot(rock_melt[y > 0.1 & y < 0.101], aes(x = x, y = relation, color = phase))
   scale_color_brewer(type = "qual", palette = "Set1")
 ```
 
-And then the pattern looks like:
+and the resulting pattern looks like:
 
 ```{r}
 image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)