Skip to content

Commit 51a168e

Browse files
authored
Merge branch 'bk-geodynamics-tutorial_1' into bk-geodynamics-tutorial_1
2 parents bf05bca + 2f32f41 commit 51a168e

13 files changed

+51392
-23251
lines changed

docs/make.jl

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -103,16 +103,17 @@ makedocs(;
103103
"18 - Alpine data integration" => "man/Tutorial_AlpineData.md",
104104
"19 - Jura tutorial" => "man/Tutorial_Jura.md",
105105
"20 - 2D model setups" => "man/Tutorial_NumericalModel_2D.md",
106-
"21 - 3D model setups" => "man/Tutorial_NumericalModel_3D.md"
106+
"21 - 3D model setups" => "man/Tutorial_NumericalModel_3D.md",
107+
"22 - Build geometry from polygons" => "man/tutorial_Polygon_structures.md"
107108
],
108109
"User Guide" => Any[
109110
"Installation" => "man/installation.md",
110111
"Data Structures" => "man/datastructures.md",
111112
"Data Import" => "man/dataimport.md",
112113
"Projection" => "man/projection.md",
113-
"Surfaces" => "man/surfaces.md",
114114
"Paraview output" => "man/paraview_output.md",
115115
"Paraview collection" => "man/paraview_collection.md",
116+
"Surfaces" => "man/surfaces.md",
116117
"Tools" => "man/tools.md",
117118
"Visualisation" => "man/visualise.md",
118119
"Gravity code" => "man/gravity_code.md",
24 KB
Loading
79.2 KB
Loading
72.3 KB
Loading

docs/src/man/Tutorial_NumericalModel_3D.md

Lines changed: 44 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,8 @@ Temp = fill(1350.0, nx,ny,nz);
3838
Much of the options are explained in the 2D tutorial, which can directly be transferred to 3D.
3939
Therefore, we will start with a simple subduction setup, which consists of a horizontal part that has a mid-oceanic ridge on one explained
4040

41-
We use a lithospheric structure. Note that if the lowermost layer has the same phase as the mantle, you can define `Tlab` as the lithosphere-asthenosphere boundary which will automatically adjust the phase depending on temperature
41+
We use a lithospheric structure, which can be specified with the `LithosphericPhases` structure, where you can indicate `Layers` (the thickness of each lithospheric layer, starting from the top), and `Phases` the phases of the corresponding layers.
42+
Note that if the lowermost layer has the same phase as the mantle, you can define `Tlab` as the lithosphere-asthenosphere boundary which will automatically adjust the phase depending on temperature
4243

4344
```julia
4445
lith = LithosphericPhases(Layers=[15 45 10], Phases=[0 1 2], Tlab=1250)
@@ -86,26 +87,35 @@ Saved file: Grid3D_FreeSubduction.vts
8687

8788
#### More sophisticated setup
8889

89-
Next, lets consider a somewhat more complicated setup
90+
Next, lets consider a somewhat more complicated setup with curved slabs, an overriding plate and a thermal structure that transitions from half-space cooling to a slab that is heated from both sides
9091

9192
```julia
93+
nx,ny,nz = 512,512,128
94+
x = range(-1000,1000, nx);
95+
y = range(-1000,1000, ny);
96+
z = range(-660,0, nz);
97+
Grid = CartData(XYZGrid(x,y,z));
98+
9299
Phases = fill(2,nx,ny,nz);
93100
Temp = fill(1350.0, nx,ny,nz);
94101
```
95102

96-
Overriding plate
103+
Overriding plate with a 30 km crust and mantle lithosphere that where T<1250 celsius
97104

98105
```julia
99-
lith_cont = LithosphericPhases(Layers=[30 200 10], Phases=[3 4 2], Tlab=1250)
100-
AddBox!(Phases, Temp, Grid; xlim=(200,1000), ylim=(-1000, 0.0), zlim=(-240.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=150));
106+
lith_cont = LithosphericPhases(Layers=[30 200 50], Phases=[3 4 2], Tlab=1250)
107+
AddBox!(Phases, Temp, Grid; xlim=(400,1000), ylim=(-1000, 0.0), zlim=(-240.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=150));
108+
AddBox!(Phases, Temp, Grid; xlim=(200,1000), ylim=(-1000, 0.0), zlim=(-80.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=150));
101109

102110
lith_cont = LithosphericPhases(Layers=[30 200 10], Phases=[5 6 2], Tlab=1250)
103-
AddBox!(Phases, Temp, Grid; xlim=(200,1000), ylim=(0, 1000), zlim=(-240.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=200));
111+
AddBox!(Phases, Temp, Grid; xlim=(400,1000), ylim=(0, 1000), zlim=(-240.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=200));
112+
AddBox!(Phases, Temp, Grid; xlim=(200,1000), ylim=(0, 1000), zlim=( -80.0, 0.0), phase = lith_cont, T=HalfspaceCoolingTemp(Age=200));
104113
```
105114

106-
Define a ridge with transform faults
115+
Define an oceanic plate with ridge
107116

108117
```julia
118+
v_spread_cm_yr = 3 #spreading velocity
109119
lith = LithosphericPhases(Layers=[15 45 10], Phases=[0 1 2], Tlab=1250)
110120
AddBox!(Phases, Temp, Grid; xlim=(-800 , 200), ylim=(-1000, -400.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3));
111121
AddBox!(Phases, Temp, Grid; xlim=(-1000,-800), ylim=(-1000, -400.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3,MORside="right"));
@@ -117,21 +127,39 @@ AddBox!(Phases, Temp, Grid; xlim=(-650, 200), ylim=(200, 1000.0), zlim=(-80.0,
117127
AddBox!(Phases, Temp, Grid; xlim=(-1000,-650), ylim=(200, 1000.0), zlim=(-80.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=3,MORside="right"));
118128
```
119129

120-
Inclined part of slab
130+
Subducting parts of the oceanic plate
131+
132+
The starting thermal age at the trench is that of the horizontal part of the oceanic plate (which is different along-trench, given that we have 3 mid oceanic ridge segments!):
133+
We want to add a smooth transition from a halfspace cooling 1D thermal profile to a slab that is heated by the surrounding mantle below a decoupling depth `d_decoupling`.
121134

122135
```julia
123-
lith = LithosphericPhases(Layers=[15 15 45 40], Phases=[7 0 1 2], Tlab=1250)
136+
AgeTrench_Myrs = 1000*1e3/(v_spread_cm_yr/1e2)/1e6 #plate age @ trench
137+
trench1 = Trench(Start=(200.0,-1000.0), End=(200.0,-400.0), Thickness=90.0, θ_max=45.0, Length=600, Lb=200, WeakzoneThickness=15, WeakzonePhase=7, d_decoupling=175);
138+
T_slab = LinearWeightedTemperature( F1=HalfspaceCoolingTemp(Age=AgeTrench_Myrs), F2=McKenzie_subducting_slab(Tsurface=0,v_cm_yr=v_spread_cm_yr, Adiabat = 0.0))
139+
addSlab!(Phases, Temp, Grid, trench1, phase = lith, T=T_slab);
140+
```
124141

125-
AddBox!(Phases, Temp, Grid; xlim=(200,700), ylim=(-1000, 0.0), zlim=(-200.0, 15), phase = lith, T=ConstantTemp(T=1350), DipAngle=30);
126-
AddBox!(Phases, Temp, Grid; xlim=(200,700), ylim=(0, 100.0), zlim=(-200.0, 15), phase = lith, T=ConstantTemp(T=1350), DipAngle=35);
142+
```julia
143+
AgeTrench_Myrs = (900)*1e3/(v_spread_cm_yr/1e2)/1e6 #plate age @ trench
144+
trench1 = Trench(Start=(200.0,-400.0), End=(200.0,200.0), Thickness=90.0, θ_max=45.0, Length=600, Lb=200, WeakzoneThickness=15, WeakzonePhase=7, d_decoupling=175);
145+
T_slab = LinearWeightedTemperature( F1=HalfspaceCoolingTemp(Age=AgeTrench_Myrs), F2=McKenzie_subducting_slab(Tsurface=0,v_cm_yr=v_spread_cm_yr, Adiabat = 0.0))
146+
addSlab!(Phases, Temp, Grid, trench1, phase = lith, T=T_slab);
147+
```
127148

128-
AddBox!(Phases, Temp, Grid; xlim=(200,700), ylim=(-1000, 0.0), zlim=(-110.0, 15), phase = lith, T=McKenzie_subducting_slab(Tsurface=0,v_cm_yr=3), DipAngle=30);
129-
AddBox!(Phases, Temp, Grid; xlim=(200,700), ylim=(0, 1000), zlim=(-110.0, 15), phase = lith, T=McKenzie_subducting_slab(Tsurface=0,v_cm_yr=3), DipAngle=35);
149+
```julia
150+
AgeTrench_Myrs = 850e3/(v_spread_cm_yr/1e2)/1e6 #plate age @ trench
151+
trench1 = Trench(Start=(200.0,200.0), End=(200.0,1000.0), Thickness=90.0, θ_max=45.0, Length=600, Lb=200, WeakzoneThickness=15, WeakzonePhase=7, d_decoupling=175);
152+
T_slab = LinearWeightedTemperature( F1=HalfspaceCoolingTemp(Age=AgeTrench_Myrs), F2=McKenzie_subducting_slab(Tsurface=0,v_cm_yr=v_spread_cm_yr, Adiabat = 0.0))
153+
addSlab!(Phases, Temp, Grid, trench1, phase = lith, T=T_slab);
154+
```
130155

131-
ind=findall(isnan.(Temp)); Temp[ind] .= 0;
156+
Finally, it is often nice to see the deformation of the plate when it subducts. A simple way to do that is to put a `stripes` on top using `addStripes`, which has the same phase as the subducting crust.
157+
158+
```julia
159+
addStripes!(Phases, Grid; stripAxes = (1,1,0), phase = ConstantPhase(0), stripePhase = ConstantPhase(9), stripeWidth=50, stripeSpacing=200)
132160
```
133161

134-
Add them to the `CartData` dataset:
162+
Finally, we can add all this to the `CartData` dataset:
135163

136164
```julia
137165
Grid = addField(Grid,(;Phases, Temp))
@@ -143,6 +171,7 @@ Saved file: Grid3D_Ridges.vts
143171
144172
````
145173

174+
And the resulting image looks like
146175
![Mechanical3D_Tutorial_2](../assets/img/Mechanical3D_Tutorial_2.png)
147176

148177
---
Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
# Adding complex geometries to a model setup including sedimentary basins, lithospheric thinning and an accretionary prism
2+
3+
4+
## Goal
5+
6+
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.
7+
8+
## Steps
9+
10+
#### 1. Set up your simplified background model
11+
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.
12+
13+
```julia
14+
using GeophysicalModelGenerator
15+
16+
# number of cells in every direction
17+
nx = 100
18+
ny = 100
19+
nz = 200
20+
21+
# define domain size
22+
x = LinRange(0.0,800.0,nx)
23+
y = LinRange(0.0,800.0,ny)
24+
z = LinRange(-660,50,nz)
25+
Cart = CartData(XYZGrid(x, y, z))
26+
27+
# initialize phase and temperature matrix
28+
Phase = ones(Int32,size(X))
29+
Temp = ones(Float64,size(X))*1350
30+
31+
# add different phases: crust->2, Mantle Lithosphere->3 Mantle->1
32+
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))
33+
34+
# add air phase 0
35+
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))
36+
```
37+
38+
39+
#### 2. Add polygon structure
40+
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)).
41+
42+
43+
```julia
44+
# xlim: x-coordinates of the points, same ordering as zlim
45+
# zlim: z-coordinates of the points, same ordering as xlim
46+
# ylim: limits the object within the two ylim values
47+
# unlimited number of points possible to create the polygon
48+
49+
# add sediment basin
50+
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));
51+
52+
# add thinning of the continental crust attached to the slab and its thickness
53+
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));
54+
55+
# add accretionary prism
56+
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));
57+
```
58+
59+
#### 3. Export final model setup to a paraview file
60+
For visualisation and comparison to actual measured data, the mode setup is saved to a paraview file.
61+
62+
```julia
63+
# # Save data to paraview:
64+
Data_Final = CartData(X,Y,Z,(Phase=Phase,Temp=Temp));
65+
Write_Paraview(Data_Final, "Sedimentary_basin");
66+
```
67+
68+
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.
69+
70+
![Tutorial_Polygon_structures](../assets/img/Tutorial_Polygon_structures.png)
71+
72+
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)

src/Setup_geometry.jl

Lines changed: 79 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ import Base: show
1111
# These are routines that help to create input geometries, such as slabs with a given angle
1212
#
1313

14-
export AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!, addSlab!, addStripes!,
14+
export AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!, AddLayer!, addPolygon!, addSlab!, addStripes!,
1515
makeVolcTopo,
1616
ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp, LithosphericTemp,
1717
ConstantPhase, LithosphericPhases,
@@ -578,6 +578,81 @@ function Rot3D!(X,Y,Z, StrikeAngle, DipAngle)
578578
return nothing
579579
end
580580

581+
582+
583+
"""
584+
addPolygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Vector(), ylim=Vector(2), zlim=Vector(), phase = ConstantPhase(1), T=nothing )
585+
586+
Adds a polygon with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models
587+
588+
Parameters
589+
====
590+
- `Phase` - Phase array (consistent with Grid)
591+
- `Temp` - Temperature array (consistent with Grid)
592+
- `Grid` - Grid structure (usually obtained with ReadLaMEM_InputFile)
593+
- `xlim` - `x`-coordinate of the polygon points, same ordering as zlim, number of points unlimited
594+
- `ylim` - `y`-coordinate, limitation in length possible (two values (start and stop))
595+
- `zlim` - `z`-coordinate of the polygon points, same ordering as xlim, number of points unlimited
596+
- `phase` - specifies the phase of the box. See `ConstantPhase()`
597+
- `T` - specifies the temperature of the box. See `ConstantTemp()`,`LinearTemp()`,`HalfspaceCoolingTemp()`,`SpreadingRateTemp()`
598+
599+
Example
600+
========
601+
602+
Polygon with constant phase and temperature:
603+
604+
```julia
605+
julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
606+
LaMEM Grid:
607+
nel : (32, 32, 32)
608+
marker/cell : (3, 3, 3)
609+
markers : (96, 96, 96)
610+
x ϵ [-3.0 : 3.0]
611+
y ϵ [-2.0 : 2.0]
612+
z ϵ [-2.0 : 0.0]
613+
julia> Phases = zeros(Int32, size(Grid.X));
614+
julia> Temp = zeros(Float64, size(Grid.X));
615+
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))
616+
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
617+
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup") # Save model to paraview
618+
1-element Vector{String}:
619+
"LaMEM_ModelSetup.vts"
620+
```
621+
622+
"""
623+
function addPolygon!(Phase, Temp, Grid::AbstractGeneralGrid; # required input
624+
xlim::Vector=[], ylim::Vector=[], zlim::Vector=[], # limits of the box
625+
phase = ConstantPhase(1), # Sets the phase number(s) in the box
626+
T=nothing ) # Sets the thermal structure (various functions are available)
627+
628+
# Retrieve 3D data arrays for the grid
629+
X,Y,Z = coordinate_grids(Grid)
630+
631+
ind = zeros(Bool,size(X))
632+
ind_slice = zeros(Bool,size(X[:,1,:]))
633+
634+
# find points within the polygon, only in 2D
635+
for i = 1:size(Y)[2]
636+
if Y[1,i,1] >= ylim[1] && Y[1,i,1]<=ylim[2]
637+
inPolygon!(ind_slice, xlim,zlim, X[:,i,:], Z[:,i,:])
638+
ind[:,i,:] = ind_slice
639+
else
640+
ind[:,i,:] = zeros(size(X[:,1,:]))
641+
end
642+
end
643+
644+
645+
# Compute thermal structure accordingly. See routines below for different options
646+
if T != nothing
647+
Temp[ind] = Compute_ThermalStructure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T)
648+
end
649+
650+
# Set the phase. Different routines are available for that - see below.
651+
Phase[ind] = Compute_Phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase)
652+
653+
return nothing
654+
end
655+
581656
"""
582657
xrot, yrot, zrot = Rot3D(X::Number,Y::Number,Z::Number, cosStrikeAngle, sindStrikeAngle, cosDipAngle, sinDipAngle)
583658
@@ -830,7 +905,6 @@ function Compute_ThermalStructure(Temp, X, Y, Z, Phase, s::SpreadingRateTemp)
830905
SecYear = 3600*24*365
831906
dz = Z[end]-Z[1];
832907

833-
834908
MantleAdiabaticT = Tmantle .+ Adiabat*abs.(Z); # Adiabatic temperature of mantle
835909

836910
if MORside=="left"
@@ -1169,7 +1243,6 @@ end
11691243
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))
11701244

11711245

1172-
11731246
"""
11741247
McKenzie_subducting_slab
11751248
@@ -1200,15 +1273,16 @@ end
12001273
Compute the temperature field of a `McKenzie_subducting_slab`. Uses the analytical solution
12011274
of McKenzie (1969) ["Speculations on the consequences and causes of plate motions"]. The functions assumes
12021275
that the bottom of the slab is the coordinate Z=0. Internally the function shifts the coordinate.
1276+
12031277
Parameters
1278+
12041279
=============================
12051280
Temp Temperature array
12061281
- `X` X Array
12071282
- `Y` Y Array
12081283
- `Z` Z Array
12091284
- `Phase` Phase array
12101285
- `s` McKenzie_subducting_slab
1211-
12121286
"""
12131287
function Compute_ThermalStructure(Temp, X, Y, Z,Phase, s::McKenzie_subducting_slab)
12141288
@unpack Tsurface, Tmantle, Adiabat, v_cm_yr, κ, it = s
@@ -1668,5 +1742,4 @@ function addSlab!(Phase, Temp, Grid::AbstractGeneralGrid, trench::Trench;
16681742
end
16691743

16701744
return nothing
1671-
end
1672-
1745+
end

test/test_setup_geometry.jl

Lines changed: 23 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,8 @@ AddBox!(Phases,Temp,Grid, xlim=(2,4), zlim=(-15,-10), phase=ConstantPhase(3), Di
1414
AddEllipsoid!(Phases,Temp,Grid, cen=(4,15,-17), axes=(1,2,3), StrikeAngle=90, DipAngle=45, phase=ConstantPhase(2), T=ConstantTemp(600))
1515
@test sum(Temp[1,1,:]) 14850.0
1616

17+
18+
1719
# CartData
1820
X,Y,Z = XYZGrid(1.0:1:10.0, 11.0:1:20.0, -20:1:-10);
1921
Data = zeros(size(X));
@@ -213,6 +215,27 @@ AddBox!(Phase, Temp, Cart; xlim=(0.0,600.0),ylim=(0.0,600.0), zlim=(-80.0, 0.0),
213215

214216
Data_Final = addField(Cart,"Temp",Temp)
215217

218+
219+
# test polygon structure
220+
221+
x = LinRange(0.0,1200.0,64);
222+
y = LinRange(0.0,1200.0,64);
223+
z = LinRange(-660,50,64);
224+
Cart = CartData(XYZGrid(x, y, z));
225+
226+
# initialize phase and temperature matrix
227+
Phase = ones(Int32,(length(x),length(y),length(z)));
228+
Temp = ones(Float64,(length(x),length(y),length(z)))*1350;
229+
230+
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));
231+
232+
# add accretionary prism
233+
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))
234+
235+
@test maximum(Phase) == 8
236+
@test minimum(Temp) == 21.40845070422536
237+
@test sum(Phase) == 292736
238+
216239
# Test the Bending slab geometry
217240

218241
# Create CartGrid struct
@@ -329,4 +352,3 @@ Phases[ind] .= 0;
329352
#Grid2D = CartData(Grid2D.x.val,Grid2D.y.val,Grid2D.z.val, (;Phases, Temp))
330353
#Write_Paraview(Grid2D,"Grid2D_SubductionCurvedOverriding");
331354

332-

0 commit comments

Comments
 (0)