-
Notifications
You must be signed in to change notification settings - Fork 0
/
test1.jl
182 lines (153 loc) · 6.64 KB
/
test1.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
using LaMEM, GeophysicalModelGenerator, Plots
include("subindshapes.jl") # Include the file containing shape definitions and functions
include("custom_gmg.jl") # Include the file containing custom functions for GeophysicalModelGenerator
modelHeight = 2880.0km
T_mantle = 1350 # in Celcius
refDensity = 3200.0kg / m^3
deltaRhoMax = 80.0kg / m^3
gravity = 9.8m / s^2
refViscosity = 5e20Pas
bodyForce = deltaRhoMax * gravity
K_tau = bodyForce * modelHeight
CharUnits = GEO_units(length=modelHeight, stress=K_tau, viscosity=refViscosity)
modelHeight = modelHeight.val
resFactor = 1.0
model = Model(
Grid(
coord_x=[0, 5040, 6192, 11520],
nel_x=round.(Int, [365, 480, 435] / resFactor),
# nel_x=round.(Int, [365, 480*2, 435] / resFactor),
bias_x=[0.03125, 1.0, 32],
y=[-2.5, 2.5],
nel_y=1,
# coord_z=[-2880, -384, 192],
coord_z=[-2880, -384, 0],
nel_z=round.(Int, [320, 320] / resFactor),
bias_z=[0.125, 1.0]
),
# BoundaryConditions(noslip=[0, 0, 1, 1, 0, 0]), # (left right front back bottom top)
Scaling(CharUnits)
)
refViscosity = refViscosity.val
opViscFactor = 1.0
modelMaterials = [
# Phase(ID=0, Name="UpperMantle", eta=refViscosity, rho=3200.0),
Phase(ID=0, Name="UpperMantle", eta=refViscosity, eta0=refViscosity, n=3.5, e0=2.5e-15, rho=3200.0, eta_st=0.1 * refViscosity, eta_vp=0.1 * refViscosity),
Phase(ID=1, Name="UppperCrustIndoAustralianPlate", eta=1e3 * refViscosity, rho=3280.0, ch=6e6), #,eta_st=0.1*refViscosity,eta_vp=0.1*refViscosity),
Phase(ID=2, Name="LowerCrustIndoAustralianPlate", eta=1e3 * refViscosity, rho=3280.0, ch=50e6),
Phase(ID=3, Name="LithosphericMantleCrustIndoAustralianPlate", eta=1e3 * refViscosity, rho=3280.0, ch=350e6),
Phase(ID=4, Name="LithosphericMantleIndoAustralianPlate", eta=5e1 * refViscosity, rho=3280.0, ch=30e6),
Phase(ID=5, Name="UpperCrustIndianIndentor", eta=1e3 * refViscosity, rho=2800.0, ch=06e6),
Phase(ID=6, Name="UpperCrustIndianIndentor", eta=1e3 * refViscosity, rho=2800.0, ch=100e6),
Phase(ID=7, Name="WeakULCrustIndianIndentor", eta=1e3 * refViscosity, rho=3000.0, ch=06e6),
Phase(ID=8, Name="LowerCrustIndianIndentor", eta=1e3 * refViscosity, rho=3100.0, ch=200e6),
Phase(ID=9, Name="WeakLLCrustIndianIndentor", eta=1e3 * refViscosity, rho=3100.0, ch=06e6),
Phase(ID=10, Name="LithosphericMantleIndianIndentor", eta=1e3 * refViscosity, rho=3100.0, ch=350e6),
Phase(ID=11, Name="LithosphericMantleIndianIndentor", eta=5e1 * refViscosity, rho=3200.0, ch=30e6),
Phase(ID=12, Name="CrustEurasianPlateForeArc", eta=1e3 * refViscosity * opViscFactor, rho=2800.0),
Phase(ID=13, Name="LithosphericMantleEurasianPlateForeArc", eta=5e2 * refViscosity * opViscFactor, rho=3220.0),
Phase(ID=14, Name="CrustEurasianPlateBackArc", eta=2.5e2 * refViscosity * opViscFactor, rho=2800.0),
Phase(ID=15, Name="LithosphericMantleEurasianPlateBackArc", eta=1.25e2 * refViscosity * opViscFactor, rho=3220.0),
Phase(ID=16, Name="CrustEurasianPlate", eta=2e3 * refViscosity * opViscFactor, rho=2800.0),
Phase(ID=17, Name="LithosphericMantleEurasianPlate", eta=1e3 * refViscosity * opViscFactor, rho=3220.0),
]
slabshapes = Slab2D(
top_x=2082.0, top_y=0.0, length=3672.1, taper=30.0, dip=29.0, depth=180.0,
thickness_array=[7.5, 22.5, 30.0, 30.0])
indentorshapes = Indentor2D(
top_x=864.0, top_y=0.0, length=2448.0, taper=18.0, taper2=12.0,
thickness_array=[5.0, 7.5, 7.5, 12.5, 7.5, 30.0, 30.0,])
overRidingShapesForeArc = OverRidingPlate2D(
top_x=11520.0, top_y=0.0, length=-5760.1, taper=90.0, dip=29.0,
thickness_array=[30.0, 30.0])
overRidingShapesBackArc = OverRidingPlate2D(
top_x=11520.0, top_y=0.0, length=-5586.0, taper=90.0, dip=90.0,
thickness_array=[30.0, 30.0])
overRidingShapes = OverRidingPlate2D(
top_x=11520.0, top_y=0.0, length=-5040.0, taper=90.0, dip=17.0,
thickness_array=[40.0, 80.0])
modelshapePolygons = vcat(
slabshapes.polygons,
indentorshapes.polygons,
overRidingShapesForeArc.polygons,
overRidingShapesBackArc.polygons,
overRidingShapes.polygons)
add_phase!(model, modelMaterials[1])
model.Grid.Phases .= 0;
for (index, polygon) in enumerate(modelshapePolygons)
AddPoly!(model, vertices=polygon, phase=ConstantPhase(index))
add_phase!(model, modelMaterials[index+1])
end
model.Solver = Solver(
SolverType="direct",
DirectSolver="mumps",
DirectPenalty=1e6,
# SolverType = "multigrid",
MGLevels=4,
# MGCoarseSolver = "mumps",
# DirectSolver = "mumps",
PETSc_options=[
# Eisenstatt-Walker (sometimes converges faster for nonlinear cases)
"-snes_ksp_ew",
"-snes_ksp_ew_rtolmax 1e-1",
"-snes_max_linear_solve_fail 10000", # maximum number of linear solve failures before we give up
# "-snes_rtol 1e-3",
# "-snes_atol 1e-4",
"-snes_max_it 100", #
"-pcmat_pgamma 1e5",
# Newton/picard options
"-snes_PicardSwitchToNewton_rtol 1e0", # relative tolerance to switch to Newton (1e-2)
"-snes_NewtonSwitchToPicard_it 20", # number of Newton iterations after which we switch back to Picard
# Linesearch options
"-snes_linesearch_type l2", #Linesearch type (one of) shell basic l2 bt cp (SNESLineSearchSetType)
"-snes_linesearch_maxstep 1.0",# very important to prevent the code from "blowing up"
# Jacobian solver
"-js_ksp_type fgmres",
"-js_ksp_max_it 50",
# "-js_ksp_atol 1e-5",
# "-js_ksp_rtol 1e-4",
"-da_refine_y 1"
]
)
# model.Solver = Solver(
# SolverType="direct",
# DirectSolver="mumps",
# DirectPenalty=1e7,
# # SolverType = "multigrid",
# # MGLevels = 4,
# # MGCoarseSolver = "mumps",
# # DirectSolver = "mumps",
# PETSc_options=[
# "-snes_ksp_ew",
# "-snes_ksp_ew_rtolmax 1e-4",
# "-snes_rtol 1e-3",
# "-snes_atol 1e-4",
# "-snes_max_it 200",
# "-snes_PicardSwitchToNewton_rtol 5e-2",
# "-snes_NewtonSwitchToPicard_it 20",
# "-js_ksp_type fgmres",
# "-js_ksp_max_it 20",
# "-js_ksp_atol 1e-5",
# "-js_ksp_rtol 1e-4",
# "-snes_linesearch_type l2",
# "-snes_linesearch_maxstep 10",
# "-da_refine_y 1"
# ]
# )
model.SolutionParams = SolutionParams(
eta_min=1e-1 * refViscosity,
eta_ref=refViscosity,
eta_max=2e3 * refViscosity,
min_cohes=1e6,
# init_guess=0 # initial guess flag
# DII=1e-18
)
model.Time = Time(time_end=2000.0,
dt=0.001,
dt_min=0.00001,
dt_max=0.1,
nstep_max=400,
nstep_out=10
)
prepare_lamem(model, 64)
# run_lamem(model, 4)