[CIG-SHORT] Friction Solver optimization

Fri Aug 7 01:01:45 PDT 2015

Brad, 

I fixed one node of the edge to the right of my domain in the Y direction but it didn’t make much difference. I guess that is because gravity is turned on and Tractions balancing gravity are applied at the bottom (that should act the same, I guess).

After playing around to try to speed things up, I ended up finding the initial fault traction distribution I wanted. However, I still have a little issue (see figures attached).
For the upper and lower nodes of the fault (which belong to the free surface and the bottom boundary of the domain, respectively), normal stress deviates from the target value at the elastic pre-step. 
I would really start with homogeneous friction on the fault...
Why is that happening? Do you think it is more related to my set up (boundary conditions, initial stress conditions, integration schemes, etc) or to the solver itself? 

Cheers, 
Romain 

Romain Jolivet
Cambridge, UK
jolivetinsar at gmail.com
UK: +44 7 596 703 148
France: +33 6 52 91 76 39
US: +1 (626) 560 6356

> On Aug 5, 2015, at 18: 01, Brad Aagaard <baagaard at usgs.gov> wrote:
> 
> Romain,
> 
> It sounds like you don't have any Dirichlet constraints in the y direction. You need to constrain the y DOF of at least one node on each side of the fault to prevent rigid body motion.
> 
> What do you want to happen on the bottom boundary as you compress the domain and the fault slips? Do you want rigid body motion on each side of the fault or do you want the bottom boundary to remain where it is? This is going to dictate whether you constrain just the corners or impose some vertical motion on one side relative to the other.
> 
> Regards,
> Brad
> 
> On 08/05/2015 09:53 AM, Romain Jolivet wrote:
>> Ok, I tried this and it works… as long as I do not impose any motion on the
>> boundaries of my model.
>> Please find attached a sketch of what I am trying to do.
>> 
>> What happens is:
>> - I impose a DirichletBC on the left edge of the domain (motion toward the right
>> with 1 cm/yr velocity). Right edge is fixed horizontally with DirichletBC (it
>> can move along y). Bottom boundary has a NeumannBC with normal traction equal to
>> \rho*g*h (see sketch). Upper boundary is left free. Material is
>> elasticplanestrain with quite canonical values. Domain is prestressed to
>> withstand lithostatic pressure. Fault rheology is static friction with no
>> cohesion (\mu=0.4, here).
>> 
>> - Until friction is lower than critical friction, nothing happens on the fault,
>> the domain is compressed, the solver is blazing fast
>> 
>> - Critical friction on the fault is always reached first at the node that
>> touches the bottom of the domain, hence this node starts to slip and there, the
>> solver gets really really slow and when it gets there, only the deepest node
>> slips, not any other node (hence, I end up with crazy stress values inside my
>> domain after a few iterations).
>> 
>> I tried to implement some fieldsplit to speed things up (as you showed me a few
>> years back) and KSP gives a NANORINF norm (so I must be doing something wrong here).
>> I tried to get to critical friction faster by imposing an initial shear stress
>> equal to \mu*\sigma_0 (see sketch), but the first step does not converge (norm
>> goes up and down around 1e-1…).
>> 
>> My guess is there is something physically inconsistent in the way I set up my
>> experiment, but after a few days of tinkering, I need some external person
>> looking at it…
>> If nothing seems obvious, I’ll send my cfg files and the mesh.
>> Cheers,
>> Romain
>> 
>> Romain Jolivet
>> Cambridge, UK
>> jolivetinsar at gmail.com <mailto:jolivetinsar at gmail.com>
>> UK: +44 7596 703 148
>> France: +33 6 52 91 76 39
>> USA: +1 (626) 560 6356
>> 
>> On 4 Aug 2015, at 18:37, Brad Aagaard <baagaard at usgs.gov
>> <mailto:baagaard at usgs.gov>> wrote:
>> 
>>> Romain,
>>> 
>>> You should be able to put Neumann BC on any boundary. For example, even if you
>>> overlap a Neumann BC with a Dirichlet BC, the Neumann BC just doesn't
>>> contribute. You do need sufficient Dirichlet constraints on each side of the
>>> fault to prevent rigid body motion. Instead of a Neumann BC you can also try a
>>> Dirichlet BC with just the node on the fault removed from the Dirichlet BC.
>>> 
>>> If the above suggestions don't resolve the issue, please send a diagram of
>>> your geometry and BC.
>>> 
>>> Regards,
>>> Brad
>>> 
>>> 
>>> On 08/04/2015 06:43 AM, Romain Jolivet wrote:
>>>> Dear Brad, Matt, Charles and other Pylith developers,
>>>> 
>>>> I am currently trying to optimise the solver in an experiment I am
>>>> running and I am having some difficulties getting my run to go
>>>> faster...
>>>> 
>>>> The problem is 2D and has a fault cutting from top to bottom (with
>>>> some angle).  Since I gravity is “on”, I have to prescribe adequate
>>>> boundary conditions, particularly at the bottom of the model. As the
>>>> fault cuts the entire domain, one of its node belongs as well to the
>>>> bottom set of nodes. On this set of nodes, I therefore cannot impose
>>>> DirichletBC (otherwise, I am fixing to many degrees of freedom on the
>>>> node that also belongs to the fault). I therefore have imposed
>>>> Neumann BC, with normal traction opposing gravity forces at the
>>>> bottom of the domain.
>>>> 
>>>> When I try to use fieldsplit (to speed things up, following one of
>>>> your tutorials, a few years back), the linear solver does not run
>>>> since PETSc issues some NANORINF norm. When I remove the fault from
>>>> the model, it runs fine. I guess this is because one node of my
>>>> domain has 2 types of constraints (one Neumann BC and some frictional
>>>> BC). Is that correct? Should I keep going that way or is it a bad
>>>> idea in general (to mix Neumann and friction on the same node)?
>>>> 
>>>> Cheers, Romain
>>>> 
>>>> —————————————————————————————————————
>>>> ————————————————————————————————————— Romain Jolivet Postdoctoral
>>>> Fellow
>>>> 
>>>> University of Cambridge Department of Earth Sciences Bullard Labs
>>>> Madingley Rise Madingley Road Cambridge CB3 0EZ United Kingdom
>>>> 
>>>> email: rpj29 at cam.ac.uk <mailto:rpj29 at cam.ac.uk> Phone: +44 1223 748 938
>>>> Mobile: +44 7596 703
>>>> 148
>>>> 
>>>> France: +33 6 52 91 76 39 US: +1 (626) 560 6356
>>>> —————————————————————————————————————
>>>> —————————————————————————————————————
>>>> 
>>>> _______________________________________________ CIG-SHORT mailing
>>>> list CIG-SHORT at geodynamics.org <mailto:CIG-SHORT at geodynamics.org>
>>>> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>>>> 
>>> 
>>> _______________________________________________
>>> CIG-SHORT mailing list
>>> CIG-SHORT at geodynamics.org <mailto:CIG-SHORT at geodynamics.org>
>>> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>> 
>> 
>> 
>> _______________________________________________
>> CIG-SHORT mailing list
>> CIG-SHORT at geodynamics.org
>> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>> 
> 
> _______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/cig-short/attachments/20150807/d62fcd65/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Stresses.png
Type: image/png
Size: 37593 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/cig-short/attachments/20150807/d62fcd65/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Friction.png
Type: image/png
Size: 26720 bytes
Desc: not available
URL: <http://lists.geodynamics.org/pipermail/cig-short/attachments/20150807/d62fcd65/attachment-0003.png>