Hi all,
I’m trying to create a 3D continental extension model that includes passive plate tension and active mantle convection. However, I’m having some problems with independent models for both two conditions.
For all cases, I selected the geometry model of the box with the lithosphere boundary indicators. And I chose the visco plastic material model.
subsection Geometry model
set Model name = box with lithosphere boundary indicators
subsection Box with lithosphere boundary indicators
set X extent = 1200e3
set Y extent = 1200e3
set Z extent = 400e3
set Lithospheric thickness = 120e3
set X repetitions = 30
set Y repetitions = 30
set Z repetitions = 7
set Z repetitions lithosphere = 3
end
end
subsection Mesh refinement
set Initial global refinement = 1
set Initial adaptive refinement = 1
set Time steps between mesh refinement = 0
set Strategy = minimum refinement function
subsection Minimum refinement function
set Coordinate system = cartesian
set Variable names = x,y,z
set Function expression = if ( z>=280.e3, 2, 1)
end
set Run postprocessors on initial refinement = false
end
subsection Material model
set Model name = visco plastic
set Material averaging = pick largest
subsection Visco Plastic
# Reference temperature
set Reference temperature = 273
set Minimum strain rate = 1.e-20
set Reference strain rate = 1.e-15
# Limit the viscosity with minimum and maximum values
set Minimum viscosity = 1e18
set Maximum viscosity = 1e26
set Define thermal conductivities = true
# Upper_c Lower_c Lith_m Mantle Weak_z
set Thermal conductivities = 3.3, 2.5, 2.5, 3.0, 3.3, 3.0
set Heat capacities = 1250, 750, 750, 1250, 1250, 1250
#
set Densities = 3300, 2800, 2900, 3300, 3300, 3300
set Thermal expansivities = 3.0e-5, 2.5e-5, 2.5e-5, 3.0e-5, 3.0e-5, 3.0e-5
# Harmonic viscosity averaging
set Viscosity averaging scheme = harmonic
set Viscous flow law = composite
# Dislocation creep parameters for Upper_c Lower_c Lith_m Mantle Weak_z
set Prefactors for dislocation creep = 1.58e-17, 5.07e-18, 6.31e-22, 1.58e-17, 1.58e-17, 1.58e-17
set Stress exponents for dislocation creep = 3.5, 2.0, 3.2, 3.5, 3.5, 3.5
set Activation energies for dislocation creep = 530.e3, 167.e3, 238.e3, 530.e3, 530.e3, 530.e3
set Activation volumes for dislocation creep = 20.e-6, 0., 0., 20.e-6, 20.e-6, 20.e-6
set Prefactors for diffusion creep = 2.46e-16, 5.07e-18, 6.31e-22, 2.46e-16, 2.46e-16, 2.46e-16
set Stress exponents for diffusion creep = 1.0, 2.0, 3.2, 1.0, 1.0, 1.0
set Activation energies for diffusion creep = 375.e3, 167.e3, 238.e3, 375.e3, 375.e3, 375.e3
set Activation volumes for diffusion creep = 10.e-6, 0., 0., 10.e-6, 10.e-6, 10.e-6
set Grain size = 1e-3
set Grain size exponents for diffusion creep = 3.0, 0.0, 0.0, 0.0, 3.0, 0.0
set Angles of internal friction = 20., 20., 20., 20., 20., 3
set Cohesions = 20.e6, 20.e6, 20.e6, 20.e6, 20.e6, 2.e6
end
end
I set up a weak zone 10km wide in the lithospheric mantle. In my study area, the direction of weak zone changes.
In the passive model, I want to keep the left boundary fixed, and set the outflow and inflow velocity at the right lithosphere boundary and mantle boundary. I set the top boundary to be a free surface.
# 0:left, 1:right, 2:front, 3:back, 4:bottom, 5:top, 6:left lithosphere, 7:right lithosphere, 8:front lithosphere, 9:back lithosphere
subsection Boundary velocity model
set Tangential velocity boundary indicators = 0, 2, 3, 4, 6, 8, 9
set Prescribed velocity boundary indicators = 1 x:function, 7 x:function
subsection Function
set Variable names = x,y,z
set Function constants = v=0.005, h=120.e3, H=260.e3, d=20.e3, W=1200.e3
set Function expression = if(z>=H+d, v, \
if(z<=H, -v*(h+d/2)/(H+d/2), ((-v*(h+d/2)/(H+d/2)-v)/d)*(H+d-z)+v)); \
0; \
0
end
end
I set up the GMG solver based on the discussion of the relevant posts in the forum, and set the solver parameters as follows.
############### Global parameters
set Dimension = 3
set Start time = 0
set End time = 30e6
set Use years in output instead of seconds = true
set Maximum time step = 0.25e5
set CFL number = .1
set Output directory = test_motion_nomantle_wz10_3D_1
set Nonlinear solver scheme = single Advection, iterated defect correction Stokes
set Nonlinear solver tolerance = 1e-2
set Max nonlinear iterations = 100
set Pressure normalization = no
set Timing output frequency = 1
subsection Solver parameters
subsection Stokes solver parameters
set Number of cheap Stokes solver steps = 2000
set Linear solver tolerance = 1e-3
set Use full A block as preconditioner = true
set GMRES solver restart length = 100
set Stokes solver type = block GMG
end
subsection Newton solver parameters
set Maximum linear Stokes solver tolerance = 1e-2
set Use Eisenstat Walker method for Picard iterations = true
end
subsection Advection solver parameters
set GMRES solver restart length = 50
end
end
Problems of the passvie model:
1.The result cannot converge when the nonlinear solver tolerance or linear Stokes solver tolerance is small(<=10^-3).
2.After increasing the tolerance, initially, each step required only around ten iterations. However, after computing for several million years, even after a full 100 iterations, the convergence doesn’t reach 10^-2, making the computation very slow.
############################################################################
In the active model, I set all boundaries except the top boundary to free slip. I set up a thermal anomaly with a radius of 100km in the center of the bottom boundary. The solver and mesh parameters are the same as above.
subsection Initial temperature model
set List of model names = function
subsection Function
set Variable names = x,y,z
set Function constants = h=400.e3,H=280.e3,ts1=273,ts2=633,ts3=893,ts4=1593, \
A1=1.00e-6,A2=0.25e-6,A3=0.0, \
k1=2.5,k2=2.5,k3=3.0, \
qs1=0.055,qs2=0.035,qs3=0.02625
set Function expression = if( (h-z)<=20.e3, \
ts1 + (qs1/k1)*(h-z) - (A1*(h-z)*(h-z))/(2.0*k1), \
if( (h-z)>20.e3 && (h-z)<=40.e3, \
ts2 + (qs2/k2)*(h-z-20.e3) - (A2*(h-z-20.e3)*(h-z-20.e3))/(2.0*k2), \
if( (h-z)>40.e3 && (h-z)<=120.e3, \
ts3 + (qs3/k3)*(h-z-40.e3) - (A3*(h-z-40.e3)*(h-z-40.e3))/(2.0*k3) ,\
if( (h-z)>120.e3 && z>0, \
ts4 + 140 * (1.0-z/H), \
if( z==0 && sqrt((x-600.e3)^2+(y-600.e3)^2)>100.e3 ,\
ts4 + 140 * (1.0-z/H), \
ts4 + 140 * (1.0-z/H) + 150 ) ) ) ) );
end
end
Problems of the active model:
- The model cannot operate normally; it interrupts and displays an out-of-memory error after running for two to three hours.
I do not know whether the above problems are due to the complexity of my model settings, or the solver and grid settings are inappropriate. I hope you can give me some advice and help.
Sincerely grateful!
Yingjie
