<div class="gmail_extra">On Thu, Apr 26, 2012 at 9:24 AM, Barry Smith <span dir="ltr">&lt;<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</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 class="HOEnZb"><div class="h5"><br>
On Apr 26, 2012, at 6:28 AM, Dominik Szczerba wrote:<br>
<br>
&gt; Jed, Barry, Matt:<br>
&gt;<br>
&gt; Below a collection of all pending answers in one mail.<br></div></div></blockquote><div><br></div><div>You can use -pc_type lu to test for a rank deficient preconditioning matrix.</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"><div class="HOEnZb"><div class="h5">
&gt;&gt; How do you define &quot;right results&quot;?<br>
&gt;<br>
&gt; I run the exact same problem on in a commercial FEM package and<br>
&gt; compare the results. The test problem is diffusion-reaction, with a<br>
&gt; nonlinear reaction term of a polynomial form (tested up to order 5). I<br>
&gt; Form the Jacobian as follows. We have a problem Ax = b(x), with A<br>
&gt; constant and b depending pointwise on x (i.e. db/dxn does not depend<br>
&gt; on xm, m!=n). Therefore I provide a simple Jacobian equal to A with<br>
&gt; the diagonal augmented with -db/dx. The relevant code is below. I do<br>
&gt; hope I do something wrong and you can point it out.<br>
&gt;<br>
&gt;       ierr = MatCopy(that-&gt;system.A, *B, SAME_NONZERO_PATTERN); CHKERRQ(ierr);<br>
&gt;<br>
&gt;       ierr = VecScale(that-&gt;dbdx, -1.0); CHKERRQ(ierr);<br>
&gt;       ierr = MatDiagonalSet(*B, that-&gt;dbdx, ADD_VALUES); CHKERRQ(ierr);<br>
&gt;<br>
&gt;       ierr = MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
&gt;       ierr = MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
&gt;<br>
&gt;       if(*J != *B)<br>
&gt;       {<br>
&gt;               ierr = MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
&gt;               ierr = MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
&gt;       }<br>
&gt;<br>
&gt;       *flag = SAME_NONZERO_PATTERN;<br>
&gt;<br>
&gt;<br>
&gt;&gt; Is SNES converging quadratically when you<br>
&gt;&gt; do not use -snes_mf_operator?<br>
&gt;<br>
&gt; Without -snes_mf_operator, providing analytically derived Jacobian:<br>
&gt;<br>
&gt;  0 SNES Function norm 1.747536282885e+01<br>
&gt;  1 SNES Function norm 1.918464419679e-03<br>
&gt;  2 SNES Function norm 1.086746779744e-03<br>
&gt;  3 SNES Function norm 2.904556511511e-06<br>
&gt;  4 SNES Function norm 7.804808098453e-08<br>
&gt;<br>
&gt; Without -snes_mf_operator, passing A as the jacobian (which you call<br>
&gt; Picard), and slightly underrelaxing the linear solves (otherwise it<br>
&gt; would not converge):<br>
&gt;<br>
&gt;  0 SNES Function norm 1.747886139999e+01<br>
&gt;  1 SNES Function norm 3.495796117011e+00<br>
&gt;  2 SNES Function norm 6.991806248330e-01<br>
&gt;  3 SNES Function norm 1.398879936054e-01<br>
&gt;  4 SNES Function norm 2.809047018752e-02<br>
&gt;  5 SNES Function norm 5.843745932752e-03<br>
&gt;  6 SNES Function norm 1.594583643926e-03<br>
&gt;  7 SNES Function norm 7.712066092135e-04<br>
&gt;  8 SNES Function norm 4.932447251451e-04<br>
&gt;  9 SNES Function norm 3.217945635817e-04<br>
&gt; 10 SNES Function norm 2.129867244811e-04<br>
&gt; 11 SNES Function norm 1.399983127264e-04<br>
&gt; 12 SNES Function norm 9.246882776024e-05<br>
&gt; 13 SNES Function norm 6.088284618615e-05<br>
&gt; 14 SNES Function norm 4.017029247381e-05<br>
&gt; 15 SNES Function norm 2.646770266042e-05<br>
&gt; 16 SNES Function norm 1.745497080382e-05<br>
&gt; 17 SNES Function norm 1.150437737331e-05<br>
&gt; 18 SNES Function norm 7.585434690336e-06<br>
&gt; 19 SNES Function norm 5.000147553531e-06<br>
&gt; 20 SNES Function norm 3.296560754852e-06<br>
&gt; 21 SNES Function norm 2.173152296533e-06<br>
&gt; 22 SNES Function norm 1.432685503660e-06<br>
&gt; 23 SNES Function norm 9.444717211995e-07<br>
&gt; 24 SNES Function norm 6.226464090397e-07<br>
&gt; 25 SNES Function norm 4.104757979681e-07<br>
&gt; 26 SNES Function norm 2.706068836265e-07<br>
&gt; 27 SNES Function norm 1.783963930887e-07<br>
&gt; 28 SNES Function norm 1.176079057620e-07<br>
&gt;<br>
&gt; Dominik<br>
&gt;<br>
&gt;&gt; Run with -snes_type test and see what it reports.<br>
&gt;<br>
&gt; Running with the options:<br>
&gt;<br>
&gt; -ksp_type bcgs -pc_type jacobi -ksp_rtol 1e-4 -ksp_max_it 1000<br>
&gt; -ksp_converged_use_initial_residual_norm -ksp_norm_type<br>
&gt; unpreconditioned -ksp_monitor_true_residual -ksp_converged_reason<br>
&gt; -ksp_view -snes_monitor -snes_converged_reason -snes_ls_monitor<br>
&gt; -snes_view -snes_type test -snes_test_display<br>
&gt;<br>
&gt; I get the output:<br>
&gt;<br>
&gt; Testing hand-coded Jacobian, if the ratio is<br>
&gt; O(1.e-8), the hand-coded Jacobian is probably correct.<br>
&gt;<br>
&gt; and the procedure stagnates, calling FormJacobian only once, and then<br>
&gt; FormFunction infinitely, with my own computation of ||F|| indicating<br>
&gt; stagnation.<br>
<br>
</div></div>   It is not calling it infinitely. Run with the -snes_type test option for a small problem, with say 10 unknowns to start.<br>
<br>
   Barry<br>
<br>
&gt;<br>
&gt; Regards,<br>
&gt; Dominik<br>
<br>
</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>
</div>