[CIG-SHORT] Friction Solver optimization

Fri Aug 7 10:39:20 PDT 2015

Thanks, I will try the fieldsplit again.

For the boundary conditions, if I wait long enough (until steady state), it becomes negligible.

R

Sent with my phone, sorry for typos...

> Le 7 août 2015 à 17:25, Brad Aagaard <baagaard at usgs.gov> a écrit :
> 
> Romain,
> 
> The spatial variations in the tractions/stresses where the fault meets a free surface are a result of the discretization and BC not the solver. Usually the only way to remove these effects is to apply a shear traction to the boundary; obviously this doesn't make sense if you want a free surface.
> 
> Fixing the node on the edge should eliminate the NANORINF you get with the fieldsplit solver.
> 
> Regards,
> Brad
> 
> 
>> On 08/07/2015 01:01 AM, Romain Jolivet wrote:
>> 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 <mailto: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
>>> <mailto: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>
>>>> <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>
>>>> <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> <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>
>>>>>> <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>
>>>>> <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 <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