On Tue, Sep 18, 2012 at 7:42 PM, Jeffrey Thompson <span dir="ltr">&lt;<a href="mailto:jeffremt@gmail.com" target="_blank">jeffremt@gmail.com</a>&gt;</span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi Brad and Matt,<br><br>Thank you for your rapid responses and suggestions.  I&#39;ve spent the day looking into things, and here&#39;s what I&#39;ve managed to find out:<br><br>Brad, in terms of the mesh, I&#39;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.<br>

<br>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.  <br>
</blockquote><div><br></div><div>Okay, Brad explained to me more about your problem. The nonlinearity is coming from frictional constraints, rather than</div><div>material nonlinearities. I now think that we are destroying quadratic convergence, but it is hard to maintain. We are using</div>
<div>a simple gradient method which approximates the true Jacobian and does not handle the constraints exactly either. This</div><div>can result in linear convergence of the nonlinear system. We are investigating more sophisticated solvers (such as those</div>
<div>used for variational inequalities in contact problems), but those will take some doing to get working for this problem. Thus,</div><div>I think the best thing that can be done is to get as good an initial guess as possible right now.</div>
<div><br></div><div>   Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Thanks again for looking into this for me and for the suggestions,<br>Jeff<div class="HOEnZb">
<div class="h5"><br><br><div class="gmail_quote">On Tue, Sep 18, 2012 at 1:01 PM, Brad Aagaard <span dir="ltr">&lt;<a href="mailto:baagaard@usgs.gov" target="_blank">baagaard@usgs.gov</a>&gt;</span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Jeff,<br>
<br>
I suspect your slow convergence (in addition to any solver settings) is<br>
due to your spatial variation of discretization size. Why do you have<br>
such a fine discretization over the entire middle portion of the domain<br>
and then rapid coarsening right near the top and bottom (+y and -y)<br>
edges? You only need fine resolution where there are large strain gradients.<br>
<br>
Carrying around so many extra DOF unnecessarily increases the problem<br>
size and then coarsening right at the edges introduces distorted cells<br>
which decreases the rate of convergence.<br>
<br>
The length scale and shear modulus used in the nondimensionalization<br>
should be comparable to the discretization size and shear modulus in<br>
your problem. Picking ones that minimize the minimum residual is<br>
irrelevant. You want the nondimensionalized displacements and stresses<br>
to on the order of 1.0 in order to produce a well-conditioned system.<br>
<span><font color="#888888"><br>
Brad<br>
</font></span><div><div><br>
<br>
<br>
<br>
On 09/18/2012 11:30 AM, Jeffrey Thompson wrote:<br>
&gt; Hi Cig-Short,<br>
&gt;<br>
&gt; I&#39;m having an issue where I&#39;m only able to achieve a slow rate of linear<br>
&gt; convergence using the nonlinear solver, resulting in prohibitively long run<br>
&gt; times.  I&#39;ve attached a tarball with the relevant configuration, mesh, and<br>
&gt; spatialdb files to run my problem.<br>
&gt;<br>
&gt; The gist of the problem I&#39;m trying run is an extension of the dyke-opening<br>
&gt; example in the manual, with a horizontal opening under the action of a<br>
&gt; constant pressure (normal traction) distribution between a viscoelastic and<br>
&gt; elastic layer.  This is a two dimensional model.  I&#39;m using the following<br>
&gt; preconditioner scheme to insure rapid convergence of the linear solver (1-2<br>
&gt; iterations):<br>
&gt;<br>
&gt; [pylithapp.petsc]<br>
&gt; ksp_rtol = 1.0e-15<br>
&gt; pc_type = lu<br>
&gt; sub_pc_factor_shift_positive_definite = 0<br>
&gt; sub_pc_factor_shift_nonzero =<br>
&gt;<br>
&gt; fs_pc_type = fieldsplit<br>
&gt; fs_pc_fieldsplit_real_diagonal = true<br>
&gt; fs_pc_fieldsplit_type = schur<br>
&gt; fs_pc_fieldsplit_schur_factorization_type = full<br>
&gt;<br>
&gt; fs_fieldsplit_0_pc_type = lu<br>
&gt; fs_fieldsplit_0_ksp_type = preonly<br>
&gt; fs_fieldsplit_1_pc_type = jacobi<br>
&gt; fs_fieldsplit_1_ksp_type = gmres<br>
&gt; fs_fieldsplit_1_ksp_rtol = 1.0e-15<br>
&gt; fs_fieldsplit_2_pc_type = jacobi<br>
&gt; fs_fieldsplit_2_ksp_type = gmres<br>
&gt; fs_fieldsplit_2_ksp_rtol = 1.0e-15<br>
&gt; fs_fieldsplit_3_pc_type = jacobi<br>
&gt; fs_fieldsplit_3_ksp_type = gmres<br>
&gt; fs_fieldsplit_3_ksp_rtol = 1.0e-15<br>
&gt;<br>
&gt; log_summary = true<br>
&gt; ksp_max_it = 100<br>
&gt; ksp_gmres_restart = 50<br>
&gt; ksp_converged_reason = true<br>
&gt;<br>
&gt; snes_ls_monitor = true<br>
&gt; snes_rtol = 1.0e-8<br>
&gt; snes_atol = 1.0e-7<br>
&gt; snes_max_it = 10000<br>
&gt; snes_monitor = true<br>
&gt; snes_converged_reason = true<br>
&gt;<br>
&gt; However, as I mentioned earlier, the SNES residual only drops linearly with<br>
&gt; the default solving scheme past the first SNES iteration.  I&#39;ve tried using<br>
&gt; different PETSc options for SNES, but none of the scheme I&#39;ve tried had<br>
&gt; better results (and many were worse or never converged).  To summarize<br>
&gt; these attempts:<br>
&gt; NGMRES: did not converge (residual remained constant)<br>
&gt; NCG: Tried the 5 different flavours, all either diverged or never converged<br>
&gt; (bounced between a few residuals)<br>
&gt; other LINESEARCH options: none were as good as the default behavior<br>
&gt;<br>
&gt; I&#39;ve also tried different non-dimensionalization options and found a set<br>
&gt; that gives the optimal initial residual (i.e., the lowest value of SNES<br>
&gt; residual number 0).  It seems like for this problem, only the<br>
&gt; nondimensional length-scale directly impacts the SNES residual.<br>
&gt;<br>
&gt; The result of all of this is that I currently need around 800-1000<br>
&gt; non-linear iterations to converge on a single time-step at best.  I have a<br>
&gt; feeling that I must be doing something wrong, but I can&#39;t tell if this is a<br>
&gt; solver issue or some other type of problem (bad mesh, poorly defined<br>
&gt; boundary conditions, etc).  Any advice would be appreciated!<br>
&gt;<br>
&gt; Thanks for the help!<br>
&gt; Jeff Thompson<br>
&gt; Caltech Seismolab<br>
&gt;<br>
&gt;<br>
&gt;<br>
</div></div><div><div>&gt; _______________________________________________<br>
&gt; CIG-SHORT mailing list<br>
&gt; <a href="mailto:CIG-SHORT@geodynamics.org" target="_blank">CIG-SHORT@geodynamics.org</a><br>
&gt; <a href="http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short" target="_blank">http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short</a><br>
&gt;<br>
<br>
_______________________________________________<br>
CIG-SHORT mailing list<br>
<a href="mailto:CIG-SHORT@geodynamics.org" target="_blank">CIG-SHORT@geodynamics.org</a><br>
<a href="http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short" target="_blank">http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short</a><br>
</div></div></blockquote></div><br>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener<br>