[CIG-SHORT] Friction Solver optimization

Fri Aug 7 09:25:58 PDT 2015

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
>