Hi John,

@jbnaliboff - I made couple of test, ASPECT 2.5 and 2.6 really on the same output (please see my answer to Wolfgang). Thank you very much for your detailed explanation and help regarding parameters. After making the edits you suggested, the code looks like below,

```
# Global parameters
set Dimension = 2
set Start time = 0
set End time = 150e6
set Use years in output instead of seconds = true
set Nonlinear solver scheme = single Advection, iterated Stokes
set Nonlinear solver tolerance = 1e-4
set Max nonlinear iterations = 120
set CFL number = 0.1
set Output directory = new_aspect26_output
set Pressure normalization = no
set Timing output frequency = 1
# Model geometry (200x60 km)
subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 100 # for test
set Y repetitions = 20 # for test
set X extent = 200e3
set Y extent = 60e3
end
end
# No mesh refinement, but the coarse mesh is 200x100
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
set Time steps between mesh refinement = 0
end
# The nonlinear solver typically converges more quickly
# when no cheap Stokes solver steps are used for
# problems with Drucker-Prager plasticity.
subsection Solver parameters
subsection Stokes solver parameters
set Stokes solver type = block AMG
set Number of cheap Stokes solver steps = 500
set Linear solver tolerance = 1e-7
set GMRES solver restart length = 100
set Use full A block as preconditioner = true
end
end
# Free surface boundary classifications.
# Advecting the free surface vertically rather than
# in the surface normal direction can result in a
# more stable mesh when the deformation is large
subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface
subsection Free surface
set Surface velocity projection = vertical
end
end
# The top boundary is open (zero traction), which allows the sticky air to
# flow freely through it as topography develops along the wedge. Additional
# testing revealed that using a true free surface leads to significant mesh
# distortion and associated numerical instabilities.
subsection Boundary traction model
set Prescribed traction boundary indicators = top: zero traction
end
# Velocity on boundaries characterized by functions
# Total extension rate is 1 cm/yr (0.5 cm/yr on each side)
subsection Boundary velocity model
set Prescribed velocity boundary indicators = left x: function, right x:function
set Tangential velocity boundary indicators = bottom
subsection Function
set Variable names = x,y
set Function constants = m=-0.05, year=1
set Function expression = if (x<100e3 , -1*m/year, 1*m/year); m*2*100e3/200e3 # Include influx
end
end
# Number and name of compositional fields
# The field plastic_strain is used for tracking the plastic finite strain invariant
# upper: brittle upper crust; seed: 'weaker' brittle region
# lower: viscous lower crust
subsection Compositional fields
set Number of fields = 4
set Names of fields = plastic_strain, upper, seed, lower
end
# Spatial domain of different compositional fields
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = 0; \
if((y>=50e3 && x<=98.0e3) || (y>=50e3 && x>=102.0e3) || \
(y>=51e3 && x>=96.0e3 && x<=104.0e3), 1, 0); \
if(y>=50e3 && y<51e3 && x>95.0e3 && x<100.0e3, 1, 0); \
if(y<50e3, 1, 0);
end
end
# Composition boundary conditions
subsection Boundary composition model
set List of model names = initial composition
end
# Use discontinuous composition bound preserving limiter
subsection Discretization
set Composition polynomial degree = 2
set Stokes velocity polynomial degree = 2
set Temperature polynomial degree = 2
set Use discontinuous composition discretization = true
subsection Stabilization parameters
set Use limiter for discontinuous composition solution = true # apply the limiter to the DG solutions
set Global composition maximum = 100.0, 1.0, 1.0, 1.0
set Global composition minimum = 0.0, 0.0, 0.0, 1.0
end
end
# Temperature boundary conditions
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 1400
set Top temperature = 550
end
end
# Temperature initial conditions (linear)
subsection Initial temperature model
set Model name = function
subsection Function
set Function expression = 550 + (1400 - 550) * y / 60e3
end
end
# Material model
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
set Reference temperature = 273
set Minimum strain rate = 1.e-20
set Reference strain rate = 1.e-16
set Minimum viscosity = 1e18
set Maximum viscosity = 1e26
set Thermal diffusivities = 1.8e-6
set Heat capacities = 750.
set Densities = 2800
set Thermal expansivities = 0.
set Viscosity averaging scheme = harmonic
set Viscous flow law = dislocation
set Prefactors for dislocation creep = 5.e-25, 5.e-25, 5.e-26, 5.e-25, 5.e-11
set Stress exponents for dislocation creep = 1.0
set Activation energies for dislocation creep = 0.
set Activation volumes for dislocation creep = 0.
set Angles of internal friction = 30., 30., 30., 4., 30.
set Cohesions = 20.e6, 20.e6, 20.e6, 2.e6, 1.e11
set Strain weakening mechanism = plastic weakening with plastic strain only
set Start plasticity strain weakening intervals = 0.5, 0.5, 0.5, 0.5, 0.5
set End plasticity strain weakening intervals = 1.5, 1.5, 1.5, 1.5, 1.5
set Cohesion strain weakening factors = 0.1, 0.1, 0.1, 1.0, 1.0
set Friction strain weakening factors = 0.13333, 0.1333, 0.1333, 1.0, 1.0
end
end
# Gravity model
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 8.87
end
end
# Post processing
subsection Postprocess
set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization
subsection Visualization
set List of output variables = density, viscosity, strain rate, error indicator, partition, surface stress, shear stress, surface dynamic topography
set Time between graphical output = 1e5
set Interpolate output = false
end
end
```

and I can run it without getting an error. However, as you can see in the attached gif file, from the *first timestep*, the mesh deforms from the edges and starts to rise.

I hope you can help me with this.

GMRES solver restart length = 150, works better.

In post-processing, when I wanted to output the particles for the material (which is not set above), I saw an incredible increase in the calculation process. I wonder if this is normal.

Thanks,

T.