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

Wed Mar 10 18:03:21 PST 2010

Here is my reply to John.

________________________________
From: Duclaux, Guillaume (CESRE, Kensington)
Sent: Thursday, 11 March 2010 9:06 AM
To: 'John Naliboff'
Cc: nicolas.riel at free.fr
Subject: RE: [CIG-LONG] Convergence issues with long-term extension models

Hi John,

There's a hidden magical compressibility function in Gale/Underworld.

In the component part of your script, create a new shape (let say called "compressibleShape") at the top of the model, then define it's rheology and mechanical properties.

You can create a new viscosity. As you suggested, a low viscosity works. You may have to do some testing. In case your model resolution is too low, you'll find an instability between the compressible layer and the upper crust with particles "flying' upward.

    <struct name="compressibleViscous">
      <param name="Type">MaterialViscosity</param>
      <param name="eta0">1e-4</param>
    </struct>

Then define the compressibility function:

<struct name="compresibility">
      <param name="Type">Compressible</param>
      <param name="GeometryMesh">mesh-linear</param>
      <param name="StiffnessMatrix">c_matrix</param>
      <param name="Swarm">picIntegrationPoints</param>
      <param name="oneOnLambda">10.0</param>             <!-this is the 'compressibility factor, 1/lambda. It has to be set between 0 and 10. As a rule of thumb, the larger 1/lambda (ie 10), the more compressible -->
    </struct>

and apply the compressible material properties to the compressibleShape  material:

    <struct name="compressibleRheology">
      <param name="Type">RheologyMaterial</param>
      <param name="Shape">compressibleShape</param>
      <param name="density">XX</param>  <!-do some testing, but 0 works usually fine -->
      <list name="Rheology">
        <param>compressibleViscous</param>
        <!-- <param>storeViscosity</param>
        <param>storeStress</param> --> <!-this is optional, dunno if you're outputting Stress and Viscosity for the particles -->
      </list>
      <param name="Compressible">compresibility</param>
    </struct>

Ge it a shot, and let me know if you have any trouble with this function.

Cheers

Gilly
________________________________________________

Dr Guillaume Duclaux
CSIRO Earth Science and Resource Engineering
Visiting address: ARRC, 26 Dick Perry Av., Kensington WA 6151
Postal address: PO Box 1130, Bentley WA 6102, Australia
Ph: + 61 8 6436 8728    Fax: + 61 8 6436 8555    Web: www.csiro.au<http://www.csiro.au>

________________________________
From: John Naliboff [mailto:jbnaliboff at gmail.com] On Behalf Of John Naliboff
Sent: Thursday, 11 March 2010 8:29 AM
To: Duclaux, Guillaume (CESRE, Kensington)
Cc: nicolas.riel at free.fr
Subject: Re: [CIG-LONG] Convergence issues with long-term extension models

Hi again Gilly,

Apologies for the second email.

When implementing the compressible layer, did you simply specify a low-viscosity/low-density region and then not subtract the hydrostatic term from this region (i.e. hydrostatic term still subtracted from the lithosphere)?

This is the only way I could think of to create a quasi-compressible layer with the current version of Gale (1.4.1).  If used another method to create the compressible layer, however, any tips would be great!

Cheers,
John

On Mar 10, 2010, at 9:32 AM, John Naliboff wrote:

Hi Gilly,

Thanks for the suggestion and I'll give the compressible layer a shot today!

John

On Mar 10, 2010, at 12:20 AM, <Guillaume.Duclaux at csiro.au<mailto:Guillaume.Duclaux at csiro.au>> wrote:

Hi John,

Have you tried adding a compressible layer at the surface of the model?

Nicolas Riel and I ran some extension models recently and found out that the only way to achieve large deformation is to use such compressible layer.

The initial geometry of the model is as follow with the compressible layer in red:

<image003.jpg>

so you can extend the model quite extensively (the compressible has been removed for visualisation)

<image004.jpg>

Give it a go!

Cheers

gilly

________________________________________________

Dr Guillaume Duclaux
CSIRO Earth Science and Resource Engineering
Visiting address: ARRC, 26 Dick Perry Av., Kensington WA 6151
Postal address: PO Box 1130, Bentley WA 6102, Australia
Ph: + 61 8 6436 8728    Fax: + 61 8 6436 8555    Web: www.csiro.au<http://www.csiro.au/>

________________________________
From: cig-long-bounces at geodynamics.org<mailto:cig-long-bounces at geodynamics.org> [mailto:cig-long-bounces at geodynamics.org] On Behalf Of John Naliboff
Sent: Wednesday, 10 March 2010 3:03 PM
To: cig-long at geodynamics.org<mailto:cig-long at geodynamics.org>
Subject: [CIG-LONG] Convergence issues with long-term extension models
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/20100311/0419f5e2/attachment-0001.htm 