# Problem with running 3D continental extension models

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:

1. 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

Can someone help me?

@egg_wang I don’t know enough about the continental extension problem to help with that, but the out of memory error is suspicious because your mesh resolution is so coarse. Can you provide the complete input file that leads to the out of memory error?
Best
W.

Hi @bangerth , thank you for your reply! I know it’s weird to run out of memory under the coarse mesh. I originally wanted to use a coarse mesh to save computation time and test if the model can compute correctly. But due to trying and modifying in the past few days, I can’t find the version that initially out of memory. The model still fails to converge. I suspect this may be related to my mesh settings. And Here is my input and error file.
test_convection_wz10_3D.prm (11.0 KB)
job-err.txt (12.5 KB)

@egg_wang I’m afraid someone else will have to help with the extension problem

@bangerth Well, thank you for your reply anyway!

Hi @egg_wang,

Thanks for posting to the forum and sorry it has taken me some time to get to this question. I saw it when it was originally posted, but did not have time until now to go properly go through and think about the setup and reported issues.

Some responses below:

1.The result cannot converge when the nonlinear solver tolerance or linear Stokes solver tolerance is small(<=10^-3).

My experience is that for this type of complex setup, the GMG solver has trouble converging when the mix/max viscosity contrasts is more than 4-5 orders of magnitude. For rifting problems, it is often the lower viscosity in the asthenosphere (1e18 or 1e19 Pa s) that seems to present the most problems, rather than the upper viscosity limit.

To proceed, I suggest trying the following:

1. Decrease the minimum viscosity to 1e20 Pas and the maximum viscosity to 1e24 Pa s.
2. Decrease the nonlinear solver tolerance to 1e-4 or 1e-5 (the stricter, the better)
3. Switch the nonlinear solver scheme to `single Advection, iterated Stokes` (see this forum post for more info).
4. Decrease the linear solver tolerance to 1e-7
5. Increasing the CFL value to 0.25 or 0.5

Run a model using the changes above and see if it produces more stable convergence behavior (it should).

From there, try expanding the min-max viscosity range (ex: 1e20-1e25, 1e19-1e24, etc) to see if you can still achieve stable convergence behavior with the same settings.

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.

I think this is due to both the constraints outlined above, and a compound effect where the system is not being solved accurately enough through time due the high non-linear solver tolerance.

1. The model cannot operate normally; it interrupts and displays an out-of-memory error after running for two to three hours.

I’m not sure what is happening here either, but let’s start with the passive model adjustments first.

Another general suggestion - if you have not already done so, try simplifying the material settings such that the viscosity in constant everywhere (min viscosity = max viscosity), and use that simplified setup to make sure the boundary and initial conditions (temperature, composition) are giving the expected results.

Can you try the suggestions above and report back?

Thanks,
John

Hi @jbnaliboff ,

Thanks for your suggestions! I modified the prm files and spend several days testing them. Although some test cases of passive model have not finished computing yet, there has been no failure to converge so far. The range of viscosity I tested included:
1. 1e20-1e24
2. 1e20-1e25
3. 1e19-1e24
4. 1e19-1e25
5. 1e18-1e24
6. 1e19-1e26
7. 1e18-1e25
8. 1e18-1e26

According to the current calculation situation, increasing the min-max viscosity range has little effect on the calculation speed of the passive model (The average calculation time for 80 cores is about 24 hours). And I also made the same parameter settings and modifications to the active model. In the viscosity range of 1e20-1e24, the model can converge normally, but high viscosity is not conducive to convection generation. When I use 1e19-1e25 , the model does not converge.I increase the Number of cheap Stokes solver steps and it doesn’t work either.

Hi @egg_wang,

Glad to hear the suggested changes allowed most of the passive and active models to converge. Do you have a sense for which parameter or combination of solver parameters enabled convergence?

When I use 1e19-1e25 , the model does not converge.I increase the Number of cheap Stokes solver steps and it doesn’t work either.

Some additional things to try:

1. Decrease the CFL / maximum time step values
2. Use material averaging over the entire cell in the material model section (ex: `set Material averaging = harmonic average only viscosity`)
3. Switch from fields to active particles, with a cell average interpolation scheme

For some examples of implementing the above features in a rift-inversion model with a 1e19-1e25 Pa s viscosity range, see the following repository that accompanies this paper. We also had a quite a bit of trouble getting stable convergence in the linear Stokes solver when the asthenosphere viscosity was reduced from 1e20 to 1e19 Pa s.

Cheers,
John

Hi @jbnaliboff ,

Since I didn’t test your suggestions one by one, I can’t pinpoint exactly which parameter contributes to convergence. Based on my experience so far,I think it should be the decrease of nonlinear/linear solver tolerance and the increase of CFL that play a rule. Because I have tried to change the solver scheme to `single Advection, iterated Stokes` before but failed to converge. I decrease the CFL and maximum time step values of the active model and it can be calculated normally at present.

I read the paper and one question I was curious about because I haven’t done it yet. How do you get the figures of the compositional fields, as shown in your paper? From the files in the repository, it appears that this is done through some python tool. I wonder if this is also possible with GMT.

Sincerely grateful,
Yingjie

Hi @egg_wang,

Completely understand and glad that the bulk updates settings enabled convergence in most of the models.

I read the paper and one question I was curious about because I haven’t done it yet. How do you get the figures of the compositional fields, as shown in your paper? From the files in the repository, it appears that this is done through some python tool. I wonder if this is also possible with GMT.

The lead author (Dylan Vasey) generated those figures with pyvista, and the linked repository has a few examples of how to generate those figures (`vtk_plot.py,` `plotting_scripts` folder). This is also possible with GMT and pyGMT, and my group is now mainly using pyGMT. We don’t currently have anything that is ready to share for pyGMT, but hope to add a few examples to the ASPECT repository in the coming months.

Cheers,
John