<br><br><div class="gmail_quote">On Mon, Feb 20, 2012 at 3:17 PM, Matthew Knepley <span dir="ltr">&lt;<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
On Mon, Feb 20, 2012 at 4:36 PM, Max Rudolph <span dir="ltr">&lt;<a href="mailto:rudolph@berkeley.edu" target="_blank">rudolph@berkeley.edu</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">

<div style="word-wrap:break-word"><pre style="line-height:normal;text-indent:0px;letter-spacing:normal;text-align:-webkit-auto;font-variant:normal;text-transform:none;font-style:normal;font-weight:normal;word-spacing:0px">
Matt, </pre><pre style="line-height:normal;text-indent:0px;letter-spacing:normal;text-align:-webkit-auto;font-variant:normal;text-transform:none;font-style:normal;font-weight:normal;word-spacing:0px">Thank you for your help.</pre>

<div><div><pre style="line-height:normal;text-indent:0px;letter-spacing:normal;text-align:-webkit-auto;font-variant:normal;text-transform:none;font-style:normal;font-weight:normal;word-spacing:0px"><blockquote type="cite">

On Mon, Feb 20, 2012 at 2:05 PM, Max Rudolph &lt;<a href="https://lists.mcs.anl.gov/mailman/listinfo/petsc-users" target="_blank">rudolph at berkeley.edu</a>&gt; wrote:

&gt;<i> Hi Dave,
</i>&gt;<i> Thanks for your help.
</i>&gt;<i>
</i>&gt;<i> Max
</i>&gt;<i>
</i>&gt;<i> Hey Max,
</i>&gt;<i>
</i>&gt;<i> Without knowing anything about the specific application related to
</i>&gt;<i> your Stokes problem, or information about the mesh you are using, I
</i>&gt;<i> have a couple of questions and suggestions which might help.
</i>&gt;<i>
</i>&gt;<i>
</i>&gt;<i> The test case that I am working with is isoviscous convection, benchmark
</i>&gt;<i> case 1a from Blankenbach 1989.
</i>&gt;<i>
</i>&gt;<i> 1) If  A, is your stokes operator A = ( K,B ; B^T, 0 ), what is your
</i>&gt;<i> precondition operator?
</i>&gt;<i> Specifically, what is in the (2,2) slot in the precondioner? - i.e.
</i>&gt;<i> what matrix are you you applying -stokes_fieldsplit_1_pc_type jacobi
</i>&gt;<i> -stokes_fieldsplit_1_ksp_type preonly to?
</i>&gt;<i> Is it the identity as in the SpeedUp notes?
</i>&gt;<i>
</i>&gt;<i>
</i>&gt;<i> I think that this is the problem. The (2,2) slot in the LHS matrix is all
</i>&gt;<i> zero (pressure does not appear in the continuity equation), so I think that
</i>&gt;<i> the preconditioner is meaningless. I am still confused as to why this
</i>&gt;<i> choice of preconditioner was suggested in the tutorial, and what is a
</i>&gt;<i> better choice of preconditioner for this block? Should I be using one of
</i>&gt;<i> the Schur complement methods instead of the additive or multiplicative
</i>&gt;<i> field split?
</i>&gt;<i>
</i>
Its not suggested, it is demonstrated. Its the first logical choice, since
Jacobi gives the identity for a 0 block (see
<a href="http://www.jstor.org/pss/2158202" target="_blank">http://www.jstor.org/pss/2158202</a>). Its
not meaningless. All the better preconditioners involve either a Schur
complement (also shown in the tutorial), or an auxiliary operator which is
more
difficult to setup and thus not shown.<br></blockquote><br></pre></div></div><pre style="line-height:normal;text-indent:0px;letter-spacing:normal;text-align:-webkit-auto;font-variant:normal;text-transform:none;font-style:normal;font-weight:normal;word-spacing:0px">
Thank you for clarifying this.</pre><pre style="line-height:normal;text-indent:0px;letter-spacing:normal;text-align:-webkit-auto;font-variant:normal;text-transform:none;font-style:normal;font-weight:normal;word-spacing:0px">
<div><div><blockquote type="cite">
&gt;<i> 2) This choice
</i>&gt;<i> -stokes_fieldsplit_0_pc_type ml -stokes_fieldsplit_0_ksp_type preonly
</i>&gt;<i> may simply not be a very effective and degrade the performance of the
</i>&gt;<i> outer solver.
</i>&gt;<i> I&#39;d make the solver for the operator in the (1,1) slot much stronger,
</i>&gt;<i> for example
</i>&gt;<i>  -stokes_fieldsplit_0_ksp_type gmres
</i>&gt;<i>  -stokes_fieldsplit_0_ksp_rtol 1.0e-4
</i>&gt;<i>  -stokes_fieldsplit_0_mg_levels_ksp_type gmres
</i>&gt;<i>  -stokes_fieldsplit_0_mg_levels_pc_type bjacobi
</i>&gt;<i>  -stokes_fieldsplit_0_mg_levels_ksp_max_it 4
</i>&gt;<i>
</i>&gt;<i> Add a monitor on this solver (-stokes_fieldsplit_0_ksp_XXX) to see how
</i>&gt;<i> ML is doing.
</i>&gt;<i>
</i>&gt;<i> 3) Using -stokes_pc_fieldsplit_type MULTIPLICATIVE should reduce the
</i>&gt;<i> number of outer iterations by a factor of two, but it will use more
</i>&gt;<i> memory.
</i>&gt;<i>
</i>&gt;<i> 4) You should use a flexible Krylov method on the outer most solve
</i>&gt;<i> (-stokes_ksp_XXX) as the preconditioner is varying between each outer
</i>&gt;<i> iteration. Use -stokes_ksp_type fgmres or -stokes_ksp_type gcr
</i>&gt;<i>
</i>&gt;<i>
</i>&gt;<i> Thanks for pointing this out. I made that change.
</i>&gt;<i>
</i>&gt;<i> 5) Depending on how the physical problem is scaled
</i>&gt;<i> (non-dimensionalised), the size of the residuals associated with the
</i>&gt;<i> momentum and continuity equation make be quite different. You are
</i>&gt;<i> currently use the entire residual from (u,p) to determine when to stop
</i>&gt;<i> iterating. You might want to consider writing a monitor which examines
</i>&gt;<i> the these residuals independently.
</i>&gt;<i>
</i>&gt;<i>
</i>&gt;<i> I think that I have scaled the problem correctly. I (slowly) obtain a
</i>&gt;<i> sufficiently accurate solution using as options only:
</i>&gt;<i> -stokes_ksp_atol 1e-5 -stokes_ksp_rtol 1e-5
</i>&gt;<i> -stokes_ksp_monitor_true_residual -stokes_ksp_norm_type UNPRECONDITIONED
</i>&gt;<i>
</i>
How do you know the problem is scaled correctly? Have you looked at norms
of the residuals for the two systems</blockquote></div></div><blockquote type="cite">  Thanks,

     Matt


&gt;<i> Cheers,
</i>&gt;<i>  Dave</i></blockquote></pre><div><br></div>Yes, here are the norms computed for the P, X, and Y components, following the last residual that ksp_monitor_true_residual returned:<br><br>383 KSP unpreconditioned resid norm 1.121628211019e-03 true resid norm 1.121628224178e-03 ||r(i)||/||b|| 9.626787321554e-10<br>

P, X, Y residual norms 5.340336e-02, 4.463404e-02, 2.509621e-02<br></div></blockquote><div><br></div><div>I am more interested in the initial residuals.</div></div></blockquote><div><br></div><div>Is this the information that you&#39;re looking for? I ran with ksp_max_it -1. The Y-residual below is much larger than the X-or P-residuals, presumably because of the initial density perturbation and body force. Thanks again for helping me.</div>
<div><br></div><div><div>  Residual norms for stokes_ solve.</div><div>  0 KSP unpreconditioned resid norm 1.165111661413e+06 true resid norm 1.165111661413e+06 ||r(i)||/||b|| 1.000000000000e+00</div><div>KSP Object:(stokes_) 2 MPI processes</div>
<div>  type: gmres</div><div>    GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement</div><div>    GMRES: happy breakdown tolerance 1e-30</div><div>  maximum iterations=-1, initial guess is zero</div>
<div>  tolerances:  relative=1e-09, absolute=1e-09, divergence=10000</div><div>  right preconditioning</div><div>  using UNPRECONDITIONED norm type for convergence test</div><div>PC Object:(stokes_) 2 MPI processes</div><div>
  type: none</div><div>  linear system matrix = precond matrix:</div><div>  Matrix Object:  CmechLHS   2 MPI processes  </div><div>    type: mpiaij</div><div>    rows=2883, cols=2883</div><div>    total: nonzeros=199809, allocated nonzeros=199809</div>
<div>    total number of mallocs used during MatSetValues calls =0</div><div>      Matrix Object:      CmechLHS       2 MPI processes      </div><div>        type: mpiaij</div><div>        rows=2883, cols=2883</div><div>        total: nonzeros=199809, allocated nonzeros=199809</div>
<div>        total number of mallocs used during MatSetValues calls =0</div><div>          using I-node (on process 0) routines: found 496 nodes, limit used is 5</div><div>converged reason -3</div><div>stokes residual...</div>
<div>converged reason -3</div><div>stokes residual...</div><div>X, Y, P residual norms 0.000000e+00, 1.165112e+06, 0.000000e+00</div></div><div><br></div><div><br></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div><br></div><div>  Thanks</div><div><br></div><div>    Matt </div><span class="HOEnZb"><font color="#888888">
</font></span></div><span class="HOEnZb"><font color="#888888"><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>

</font></span></blockquote></div><br>