[CIG-SHORT] Linear Convergence Using SNES

Brad Aagaard baagaard at usgs.gov
Fri Sep 21 11:15:59 PDT 2012 

Jeff,

Matt is correct that the rate of convergence with our current solution 
scheme is rather poor for this type of problem. When deformation is 
confined very close to the fault, it is much better. In the future we 
would like to switch to a more sophisticated scheme that is more 
efficient for enforcing inequality constraints. Such methods are much 
closer to the bleeding edge of computational science and will require 
significant effort.

In looking at your .cfg files, I noticed a few issues:

* You didn't set the PETSc settings for the friction sensitivity solve. 
See the PETSc settings beginning with friction_ in 
examples/3d/hex8/step20.cfg.

* Set ksp_rtol and snes_rtol to a very small number, like 1.0e-20. The 
FaultCohesiveDyn friction implementation really requires convergence at 
an absolute level, not a relative tolerance. Set the fault 
zero_tolerance property to a value 1-2 orders of magnitude greater than 
ksp_atol (*This is absolutely essential*).

* In generating your mesh, set the sizing function for the surface 
meshes for the central portions to bias so that you get a graded mesh.

surface SurfBottomMid sizing function type bias start curve 5 6 factor 1.15
surface SurfTopMid sizing function type bias start curve 5 6 factor 1.1


As for getting the answer you want more efficiently, my suggestions are:

(1) Create the smallest mesh (fewest DOF based on grading the cell size) 
that allows you to obtain a solution at the accuracy you need

   Use a coarser mesh and select the absolute tolerances to match the 
level of accuracy needed.

(2) Use prescribed slip distributions whenever possible

   You could use a static calculation to get the slip distribution 
associated with the traction distribution. Assuming the slip 
distribution would not be affected by the viscoelastic response, you 
could apply this slip distribution in a quasi-static simulation.

(3) Run in parallel on multiple cores to reduce runtime

   After starting with a very coarse model to work out all of the kinks 
and tuning the solver parameters, increase the resolution a couple of 
times to determine the resolution required to obtain the accuracy that 
you need and then let the computer do all of the work. Using multiple 
cores will speed up the linear solves and reduce runtime.


Regards,
Brad


On 09/19/2012 11:04 AM, Jeffrey Thompson wrote:
> Hi Brad and Matt,
>
> Thanks for the help and the offer to look over a diagram of my problem.
>   I've put together a one page project summary/description with conceptual
> and model diagrams of the problem.  If you need a format other than PDF,
> just let me know.
>
> Thanks again!
> Jeff
>
>
>
> On Wed, Sep 19, 2012 at 8:16 AM, Brad Aagaard <baagaard at usgs.gov> wrote:
>
>> Jeff,
>>
>> Can you send us a diagram of the problem you are trying to solve?
>>
>> There may be some simple tricks to improving the rate of convergence by
>> either adjusting the mesh, domain, and BC.
>>
>> Brad
>>
>>
>> On 09/18/2012 05:42 PM, Jeffrey Thompson wrote:
>>
>>> Hi Brad and Matt,
>>>
>>> Thank you for your rapid responses and suggestions.  I've spent the day
>>> looking into things, and here's what I've managed to find out:
>>>
>>> Brad, in terms of the mesh, I've run a few models with a uniform mesh size
>>> over a range of element lengths, and it seems like in each case, there is
>>> still only linear convergence, albeit with fewer iterations on each time
>>> step as the mesh is made coarser.  However, even at a coarseness beyond a
>>> useful level (100 meters spacing, which is not really enough nodes (5) to
>>> capture the nature of the applied pressure distribution), I still need 77
>>> iterations to converge in the nonlinear solver.  And this is running with
>>> the nondimensional scales set equal to the relaxation time, mesh spacing,
>>> and shear modulus.  This result leads me to think that the linear nature
>>> of
>>> the convergence is perhaps not related to the mesh.  However, your points
>>> about my mesh are good and the suggestions are certainly useful in terms
>>> of
>>> reducing the time it takes the model to iterate and ultimately reach
>>> convergence.
>>>
>>> Matt, I was able to import the Jacobian for my problem on the coarsest
>>> mesh
>>> into matlab to look for singularity.  The matrix is full-rank and
>>> invertible.  The calculated condition number from the COND function is
>>> 1.7409e+05.  And this problem had linear convergence in the nonlinear
>>> solver.
>>>
>>> Thanks again for looking into this for me and for the suggestions,
>>> Jeff
>>>
>>> On Tue, Sep 18, 2012 at 1:01 PM, Brad Aagaard <baagaard at usgs.gov> wrote:
>>>
>>>   Jeff,
>>>>
>>>> I suspect your slow convergence (in addition to any solver settings) is
>>>> due to your spatial variation of discretization size. Why do you have
>>>> such a fine discretization over the entire middle portion of the domain
>>>> and then rapid coarsening right near the top and bottom (+y and -y)
>>>> edges? You only need fine resolution where there are large strain
>>>> gradients.
>>>>
>>>> Carrying around so many extra DOF unnecessarily increases the problem
>>>> size and then coarsening right at the edges introduces distorted cells
>>>> which decreases the rate of convergence.
>>>>
>>>> The length scale and shear modulus used in the nondimensionalization
>>>> should be comparable to the discretization size and shear modulus in
>>>> your problem. Picking ones that minimize the minimum residual is
>>>> irrelevant. You want the nondimensionalized displacements and stresses
>>>> to on the order of 1.0 in order to produce a well-conditioned system.
>>>>
>>>> Brad
>>>>
>>>>
>>>>
>>>>
>>>> On 09/18/2012 11:30 AM, Jeffrey Thompson wrote:
>>>>
>>>>> Hi Cig-Short,
>>>>>
>>>>> I'm having an issue where I'm only able to achieve a slow rate of linear
>>>>> convergence using the nonlinear solver, resulting in prohibitively long
>>>>>
>>>> run
>>>>
>>>>> times.  I've attached a tarball with the relevant configuration, mesh,
>>>>>
>>>> and
>>>>
>>>>> spatialdb files to run my problem.
>>>>>
>>>>> The gist of the problem I'm trying run is an extension of the
>>>>>
>>>> dyke-opening
>>>>
>>>>> example in the manual, with a horizontal opening under the action of a
>>>>> constant pressure (normal traction) distribution between a viscoelastic
>>>>>
>>>> and
>>>>
>>>>> elastic layer.  This is a two dimensional model.  I'm using the
>>>>> following
>>>>> preconditioner scheme to insure rapid convergence of the linear solver
>>>>>
>>>> (1-2
>>>>
>>>>> iterations):
>>>>>
>>>>> [pylithapp.petsc]
>>>>> ksp_rtol = 1.0e-15
>>>>> pc_type = lu
>>>>> sub_pc_factor_shift_positive_**definite = 0
>>>>> sub_pc_factor_shift_nonzero =
>>>>>
>>>>> fs_pc_type = fieldsplit
>>>>> fs_pc_fieldsplit_real_diagonal = true
>>>>> fs_pc_fieldsplit_type = schur
>>>>> fs_pc_fieldsplit_schur_**factorization_type = full
>>>>>
>>>>> fs_fieldsplit_0_pc_type = lu
>>>>> fs_fieldsplit_0_ksp_type = preonly
>>>>> fs_fieldsplit_1_pc_type = jacobi
>>>>> fs_fieldsplit_1_ksp_type = gmres
>>>>> fs_fieldsplit_1_ksp_rtol = 1.0e-15
>>>>> fs_fieldsplit_2_pc_type = jacobi
>>>>> fs_fieldsplit_2_ksp_type = gmres
>>>>> fs_fieldsplit_2_ksp_rtol = 1.0e-15
>>>>> fs_fieldsplit_3_pc_type = jacobi
>>>>> fs_fieldsplit_3_ksp_type = gmres
>>>>> fs_fieldsplit_3_ksp_rtol = 1.0e-15
>>>>>
>>>>> log_summary = true
>>>>> ksp_max_it = 100
>>>>> ksp_gmres_restart = 50
>>>>> ksp_converged_reason = true
>>>>>
>>>>> snes_ls_monitor = true
>>>>> snes_rtol = 1.0e-8
>>>>> snes_atol = 1.0e-7
>>>>> snes_max_it = 10000
>>>>> snes_monitor = true
>>>>> snes_converged_reason = true
>>>>>
>>>>> However, as I mentioned earlier, the SNES residual only drops linearly
>>>>>
>>>> with
>>>>
>>>>> the default solving scheme past the first SNES iteration.  I've tried
>>>>>
>>>> using
>>>>
>>>>> different PETSc options for SNES, but none of the scheme I've tried had
>>>>> better results (and many were worse or never converged).  To summarize
>>>>> these attempts:
>>>>> NGMRES: did not converge (residual remained constant)
>>>>> NCG: Tried the 5 different flavours, all either diverged or never
>>>>>
>>>> converged
>>>>
>>>>> (bounced between a few residuals)
>>>>> other LINESEARCH options: none were as good as the default behavior
>>>>>
>>>>> I've also tried different non-dimensionalization options and found a set
>>>>> that gives the optimal initial residual (i.e., the lowest value of SNES
>>>>> residual number 0).  It seems like for this problem, only the
>>>>> nondimensional length-scale directly impacts the SNES residual.
>>>>>
>>>>> The result of all of this is that I currently need around 800-1000
>>>>> non-linear iterations to converge on a single time-step at best.  I have
>>>>>
>>>> a
>>>>
>>>>> feeling that I must be doing something wrong, but I can't tell if this
>>>>>
>>>> is a
>>>>
>>>>> solver issue or some other type of problem (bad mesh, poorly defined
>>>>> boundary conditions, etc).  Any advice would be appreciated!
>>>>>
>>>>> Thanks for the help!
>>>>> Jeff Thompson
>>>>> Caltech Seismolab
>>>>>
>>>>>
>>>>>
>>>>> ______________________________**_________________
>>>>> CIG-SHORT mailing list
>>>>> CIG-SHORT at geodynamics.org
>>>>> http://geodynamics.org/cgi-**bin/mailman/listinfo/cig-short<http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short>
>>>>>
>>>>>
>>>> ______________________________**_________________
>>>> CIG-SHORT mailing list
>>>> CIG-SHORT at geodynamics.org
>>>> http://geodynamics.org/cgi-**bin/mailman/listinfo/cig-short<http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short>
>>>>
>>>>
>>>
>>
>

More information about the CIG-SHORT mailing list