[CIG-LONG] Convergence issues with long-term extension models

[John Naliboff]

Hello all,

I'm hoping to get some input on convergence problems I'm having with 2D long-term extension models.

The models are designed following Wijns et al. 2005 (EPSL), who modeled extension of a brittle upper crust overlying a viscous lower crust using Underworld.  The initial geometry is 160 km (width) x 60 km (depth) and the upper crust follows a yielding law similar to the DruckerPrager yield criterion.  The boundary conditions are free-slip on all sides, although a weak, low-density layer is present above the upper crust (i.e. top of lithosphere should behave in a similar fashion to a true free-surface).   Extension is driven by a fixed velocity of 3 cm/yr applied to the right side of the 2D box. The models in the paper are run up to roughly 80% extension (~ 4 Myr). 

The simulations I'm running mimic the models above very closely, although I'm using an upper free-surface rather than placing a weak, low-density layer above the upper crust and specifying free-slip on the top boundary.

In short, the problem I'm encountering reproducing the long-term extension models is that my simulations stop converging  after ~ 1-2 Myr (~ 25%-40% extension).  In terms of specific error messages, the problem is almost always related to the Stokes Uzawa Solver (ex: uzawaRhsScale has illegal value 'inf', etc).

In terms of the model output, by 1-2 Myr very large strain rate gradients have developed at the model surface roughly following horst and grabben structures.  I think this may be the source of the convergence issue, but at this point I'm really not certain.

After initially having convergence issues with a resolution of 1 km and 30 particles per cell (ppc), I systematically increased the resolution up 0.25 km with 60 ppc, but all the models still stopped converging between 1-2 Myr.
I also tried decreasing the time step by 2 (dtFactor = 0.5) and turning off the pressure and velocity interpolation onto the new mesh, neither of which had any significant effect on the total run time before the models stopped converging.  In all of the cases above I used the MUMPS solver, which performed better than the default solver.

The convergence issue is also not likely related to the specific upper crustal yielding rheology, as I've tried both DruckerPrager, VonMises and the yielding law used in the Wijns et al. models.

Finally, I tried using a free-slip top boundary with a weak, low-density layer above the upper crust, but this produced a whole different set of issues.

Thanks in advance to anyone who has some ideas/input on the convergence issues and I've attached a sample input file (resolution = 0.5 km, ppc = 60) if that helps at all.

Cheers,
John Naliboff
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://geodynamics.org/pipermail/cig-long/attachments/20100309/b3e8d6cb/attachment-0001.htm 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: wijns05epsl_fig2_res0.5_ppc60_dtf1_030210.xml
Type: application/xml
Size: 26848 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/cig-long/attachments/20100309/b3e8d6cb/attachment-0001.xml 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://geodynamics.org/pipermail/cig-long/attachments/20100309/b3e8d6cb/attachment-0002.htm