Merge pull request #73 from TatjanaWeiler/main

Add Polygon structure in pseudo3D

Author: boriskaus

docs/make.jl

@@ -101,16 +101,17 @@ makedocs(;
             "16 - Create movies" =>  "man/tutorial_time_Seismicity.md",
             "17 - Fault Density Map" =>  "man/Tutorial_FaultDensity.md",
             "18 - Alpine data integration" =>  "man/Tutorial_AlpineData.md",
-            "19 - Jura tutorial" =>  "man/Tutorial_Jura.md"
+            "19 - Jura tutorial" =>  "man/Tutorial_Jura.md", 
+            "20 - Build geometry from polygons" =>  "man/tutorial_Polygon_structures.md"
         ],
         "User Guide" => Any[
             "Installation" =>  "man/installation.md",
             "Data Structures" =>  "man/datastructures.md",
             "Data Import" =>  "man/dataimport.md",
             "Projection" =>  "man/projection.md",
-            "Surfaces" =>  "man/surfaces.md",
             "Paraview output" => "man/paraview_output.md",
             "Paraview collection" => "man/paraview_collection.md",
+            "Surfaces" =>  "man/surfaces.md",
             "Tools" => "man/tools.md",
             "Visualisation" => "man/visualise.md",
             "Gravity code" => "man/gravity_code.md",

docs/src/assets/img/Tutorial_Polygon_structures.png

@@ -0,0 +1,72 @@
+# Adding complex geometries to a model setup including sedimentary basins, lithospheric thinning and an accretionary prism
+
+
+## Goal
+
+This tutorial visualizes simplified geological as it is done here for a passive margin where lithospheric thinning, sedimentary basin and accretionary prism occur. The simplification is based on a polygon structure for a pseudo-3D model. While the structure can have a random shape in the `x`- and `z`-direction, in the `y`-direction only the extent is variable. 
+
+## Steps
+
+#### 1. Set up your simplified background model
+Before adding specific geological features, a general simplified model setup is necessary. The construction is made by using the `addBox!` function. For the model the discontinuities are in 15, 45, 145, and 945 km depth.
+
+```julia
+using GeophysicalModelGenerator
+
+# number of cells in every direction
+nx = 100
+ny = 100
+nz = 200
+
+# define domain size
+x        = LinRange(0.0,800.0,nx)
+y        = LinRange(0.0,800.0,ny)
+z        = LinRange(-660,50,nz)
+Cart     = CartData(XYZGrid(x, y, z))
+
+# initialize phase and temperature matrix
+Phase   = ones(Int32,size(X))
+Temp    = ones(Float64,size(X))*1350
+
+# add different phases: crust->2, Mantle Lithosphere->3 Mantle->1
+AddBox!(Phase, Temp, Cart; xlim=(0.0,800.0),ylim=(0.0,800.0), zlim=(-800.0,0.0), phase = LithosphericPhases(Layers=[15 30 100 800], Phases=[2 3 1 5], Tlab=1300 ), T=LinearTemp(Ttop=20, Tbot=1600))
+
+# add air phase 0
+AddBox!(Phase, Temp, Cart; xlim=(0.0,800.0),ylim=(0.0,800.0), zlim=(0.0,50.0), phase = ConstantPhase(0), T=ConstantTemp(20.0))
+```
+
+
+#### 2. Add polygon structure
+To include the geological structures of a passive margin into the model, we use polygons for depths of up to 150 km. In the example, the sediment basin shows a more trapezoidal (2D in `x`-/`z`-direction) shape, while the thinning of the plate has a more triangular (2D in `x`-/`z`-direction). More complex structures can be build using arbitrarily sized polygons in `x`- and `z`-direction, wheraes in the `y`-direction only the length can be varied (specified by two values). The `x`- and `z`-values of the points need to be in the same order for selecting the correct point (P1(1/3), P2(2/2) --> xlim(1,2), ylim(3,2)).
+
+
+```julia
+# xlim: x-coordinates of the points, same ordering as zlim
+# zlim: z-coordinates of the points, same ordering as xlim
+# ylim: limits the object within the two ylim values
+# unlimited number of points possible to create the polygon
+
+# add sediment basin 
+addPolygon!(Phase, Temp, Cart; xlim=[0.0,0.0, 160.0, 200.0],ylim=[100.0,300.0], zlim=[0.0,-10.0,-20.0,0.0], phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30));
+
+# add thinning of the continental crust attached to the slab and its thickness 
+addPolygon!(Phase, Temp, Cart; xlim=[0.0, 200.0, 0.0],ylim=[500.0,800.0], zlim=[-100.0,-150.0,-150.0], phase = ConstantPhase(5), T=LinearTemp(Ttop=1000, Tbot=1100));
+
+# add accretionary prism 
+addPolygon!(Phase, Temp, Cart; xlim=[800.0, 600.0, 800.0],ylim=[100.0,800.0], zlim=[0.0,0.0,-60.0], phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30));
+```
+
+#### 3. Export final model setup to a paraview file
+For visualisation and comparison to actual measured data, the mode setup is saved to a paraview file.
+
+```julia
+# # Save data to paraview:
+Data_Final      =   CartData(X,Y,Z,(Phase=Phase,Temp=Temp)); 
+Write_Paraview(Data_Final, "Sedimentary_basin");
+```
+
+After importing and looking at the file to paraview, some unresolved areas might be visible as they are visible in this model. That is due to the resolution and shape of the polygon. To reduce those artefacts an increase in resolution or a change of the polygon angle might help.
+
+![Tutorial_Polygon_structures](../assets/img/Tutorial_Polygon_structures.png)
+
+If you want to run the entire example, you can find the .jl code [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/blob/main/tutorial/Tutorial_polygon_geometry.jl)

docs/src/man/tutorial_Polygon_structures.md

@@ -11,7 +11,7 @@ import Base: show
 # These are routines that help to create input geometries, such as slabs with a given angle
 #
 
-export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!, addSlab!, addStripes!,
+export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!, addPolygon!, addSlab!, addStripes!,
         makeVolcTopo,
         ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp,
         ConstantPhase, LithosphericPhases,
@@ -571,6 +571,81 @@ function Rot3D!(X,Y,Z, StrikeAngle, DipAngle)
     return nothing
 end
 
+
+
+"""
+        addPolygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Vector(), ylim=Vector(2), zlim=Vector(), phase = ConstantPhase(1), T=nothing )   
+
+Adds a polygon with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models
+
+Parameters
+====
+- `Phase` - Phase array (consistent with Grid)
+- `Temp`  - Temperature array (consistent with Grid)
+- `Grid`  - Grid structure (usually obtained with ReadLaMEM_InputFile)
+- `xlim`  - `x`-coordinate of the polygon points, same ordering as zlim, number of points unlimited
+- `ylim`  - `y`-coordinate, limitation in length possible (two values (start and stop))
+- `zlim`  - `z`-coordinate of the polygon points, same ordering as xlim, number of points unlimited
+- `phase` - specifies the phase of the box. See `ConstantPhase()`
+- `T` - specifies the temperature of the box. See `ConstantTemp()`,`LinearTemp()`,`HalfspaceCoolingTemp()`,`SpreadingRateTemp()`
+
+Example
+========
+
+Polygon with constant phase and temperature:
+
+```julia
+julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
+LaMEM Grid:
+  nel         : (32, 32, 32)
+  marker/cell : (3, 3, 3)
+  markers     : (96, 96, 96)
+  x           ϵ [-3.0 : 3.0]
+  y           ϵ [-2.0 : 2.0]
+  z           ϵ [-2.0 : 0.0]
+julia> Phases = zeros(Int32,   size(Grid.X));
+julia> Temp   = zeros(Float64, size(Grid.X));
+julia> addPolygon!(Phase, Temp, Cart; xlim=(0.0,0.0, 1.6, 2.0),ylim=(0.0,0.8), zlim=(0.0,-1.0,-2.0,0.0), phase = ConstantPhase(8), T=ConstantTemp(30))
+julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
+julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
+1-element Vector{String}:
+ "LaMEM_ModelSetup.vts"
+```
+
+"""
+function addPolygon!(Phase, Temp, Grid::AbstractGeneralGrid;    # required input
+    xlim::Vector=[], ylim::Vector=[], zlim::Vector=[],          # limits of the box
+    phase = ConstantPhase(1),                                   # Sets the phase number(s) in the box
+    T=nothing )                                                 # Sets the thermal structure (various functions are available)
+
+# Retrieve 3D data arrays for the grid
+X,Y,Z = coordinate_grids(Grid)
+
+ind = zeros(Bool,size(X))
+ind_slice = zeros(Bool,size(X[:,1,:]))
+
+# find points within the polygon, only in 2D
+for i = 1:size(Y)[2]
+    if Y[1,i,1] >= ylim[1] && Y[1,i,1]<=ylim[2] 
+        inPolygon!(ind_slice, xlim,zlim, X[:,i,:], Z[:,i,:])
+        ind[:,i,:] = ind_slice
+    else
+        ind[:,i,:] = zeros(size(X[:,1,:]))
+    end
+end
+
+
+# Compute thermal structure accordingly. See routines below for different options
+if T != nothing
+   Temp[ind] = Compute_ThermalStructure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
+end
+
+# Set the phase. Different routines are available for that - see below.
+Phase[ind] = Compute_Phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase)
+
+return nothing
+end
+
 """
     xrot, yrot, zrot = Rot3D(X::Number,Y::Number,Z::Number, cosStrikeAngle, sindStrikeAngle, cosDipAngle, sinDipAngle)
 
@@ -823,7 +898,6 @@ function Compute_ThermalStructure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp)
     SecYear     =   3600*24*365
     dz          =   Z[end]-Z[1];
 
-
     MantleAdiabaticT    =   Tmantle .+ Adiabat*abs.(Z);   # Adiabatic temperature of mantle
 
     if MORside=="left"
@@ -1162,7 +1236,6 @@ end
 Compute_Phase(Phase, Temp, Grid::LaMEM_grid, s::LithosphericPhases) = Compute_Phase(Phase, Temp, Grid.X, Grid.Y, Grid.Z, s::LithosphericPhases, Ztop=maximum(Grid.coord_z))
 
 
-
 """
     McKenzie_subducting_slab
 
@@ -1193,15 +1266,16 @@ end
 Compute the temperature field of a `McKenzie_subducting_slab`. Uses the analytical solution
 of McKenzie (1969) ["Speculations on the consequences and causes of plate motions"]. The functions assumes
 that the bottom of the slab is the coordinate Z=0. Internally the function shifts the coordinate. 
+
 Parameters
+
 =============================
 Temp Temperature array
 - `X`    X Array 
 - `Y`    Y Array 
 - `Z`    Z Array 
 - `Phase` Phase array 
 - `s`    McKenzie_subducting_slab
-
 """
 function Compute_ThermalStructure(Temp, X, Y, Z,Phase, s::McKenzie_subducting_slab)
     @unpack Tsurface, Tmantle, Adiabat, v_cm_yr, κ, it = s
@@ -1661,5 +1735,4 @@ function addSlab!(Phase, Temp, Grid::AbstractGeneralGrid,  trench::Trench;
     end
 
     return nothing
-end
-
+end

src/Setup_geometry.jl

@@ -14,6 +14,8 @@ AddBox!(Phases,Temp,Grid, xlim=(2,4), zlim=(-15,-10), phase=ConstantPhase(3), Di
 AddEllipsoid!(Phases,Temp,Grid, cen=(4,15,-17), axes=(1,2,3), StrikeAngle=90, DipAngle=45, phase=ConstantPhase(2), T=ConstantTemp(600))
 @test sum(Temp[1,1,:]) ≈ 14850.0
 
+
+
 # CartData 
 X,Y,Z               =   XYZGrid(1.0:1:10.0, 11.0:1:20.0, -20:1:-10);
 Data                =   zeros(size(X));
@@ -213,6 +215,26 @@ AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0),
 
 Data_Final =   AddField(Cart,"Temp",Temp)
 
+
+# test polygon structure
+
+x        = LinRange(0.0,1200.0,64);
+y        = LinRange(0.0,1200.0,64);
+z        = LinRange(-660,50,64);
+Cart     = CartData(XYZGrid(x, y, z));
+
+# initialize phase and temperature matrix
+Phase   = ones(Int32,(length(x),length(y),length(z)));
+Temp    = ones(Float64,(length(x),length(y),length(z)))*1350;
+
+AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0), phase = ConstantPhase(5), T=T=ConstantTemp(120.0));
+
+# add accretionary prism 
+addPolygon!(Phase, Temp, Cart; xlim=[500.0, 200.0, 500.0],ylim=[100.0,400.0], zlim=[0.0,0.0,-60.0], phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
+
+@test maximum(Phase) == 8
+@test minimum(Temp) == 21.40845070422536
+
 # Test the Bending slab geometry
 
 # Create CartGrid struct
@@ -329,4 +351,3 @@ Phases[ind] .= 0;
 #Grid2D = CartData(Grid2D.x.val,Grid2D.y.val,Grid2D.z.val, (;Phases, Temp))
 #Write_Paraview(Grid2D,"Grid2D_SubductionCurvedOverriding");
 
-

test/test_setup_geometry.jl

@@ -0,0 +1,45 @@
+using GeophysicalModelGenerator
+
+
+# number of cells in every direction
+nx = 100
+ny = 100
+nz = 200
+
+# define domain size
+x        = LinRange(0.0,800.0,nx)
+y        = LinRange(0.0,800.0,ny)
+z        = LinRange(-660,50,nz)
+X,Y,Z    = XYZGrid(x, y, z);
+Cart     = CartData(X,Y,Z, (Data=Z,))
+
+# initialize phase and temperature matrix
+Phase   = ones(Int32,size(X))
+Temp    = ones(Float64,size(X))*1350
+
+# add different phases: crust->2, Mantle Lithosphere->3 Mantle->1
+AddBox!(Phase, Temp, Cart; xlim=(0.0,800.0),ylim=(0.0,800.0), zlim=(-800.0,0.0), phase = LithosphericPhases(Layers=[15 30 100 800], Phases=[2 3 1 5], Tlab=1300 ), T=LinearTemp(Ttop=20, Tbot=1600) )#T=HalfspaceCoolingTemp(Tsurface=20.0, Tmantle=1350, Age=120, Adiabat=0.4)
+
+
+# xlim: x-coordinates of the points, same ordering as zlim
+# zlim: z-coordinates of the points, same ordering as xlim
+# ylim: limits the object within the two ylim values
+# unlimited number of points possible to create the polygon
+# add sediment basin # depending on the resolution and angle if it the edge is visible in paraview
+addPolygon!(Phase, Temp, Cart; xlim=[0.0,0.0, 160.0, 200.0],ylim=[100.0,300.0], zlim=[0.0,-10.0,-20.0,0.0], phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
+
+# add thinning of the continental crust attached to the slab and its thickness 
+addPolygon!(Phase, Temp, Cart; xlim=[0.0, 200.0, 0.0],ylim=[500.0,800.0], zlim=[-100.0,-150.0,-150.0], phase = ConstantPhase(5), T=LinearTemp(Ttop=1000, Tbot=1100))
+
+# add accretionary prism 
+addPolygon!(Phase, Temp, Cart; xlim=[800.0, 600.0, 800.0],ylim=[100.0,800.0], zlim=[0.0,0.0,-60.0], phase = ConstantPhase(8), T=LinearTemp(Ttop=20, Tbot=30))
+
+
+# add air phase 0
+AddBox!(Phase, Temp, Cart; xlim=(0.0,800.0),ylim=(0.0,800.0), zlim=(0.0,50.0), phase = ConstantPhase(0), T=ConstantTemp(20.0))
+
+# # Save data to paraview:
+Data_Final      =   CartData(X,Y,Z,(Phase=Phase,Temp=Temp)) 
+Write_Paraview(Data_Final, "Sedimentary_basin")
+
+