<div class="gmail_quote">On Fri, Mar 18, 2011 at 14:03, Gianluca Meneghello <span dir="ltr"><<a href="mailto:gianmail@gmail.com">gianmail@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div id=":dr">what do you mean by constant +1 and constant -1? Is it the<br>
starting/user-defined state vector +1 and -1?<br>
<br>
My problem should not be ill posed </div></blockquote><div><br></div><div>Not ill-posed, just ill-conditioned, which can be the result of poor scaling. For example, changing the units of velocity or pressure so that these parts of the Jacobian were "the same size" might improve the conditioning of the system. Depending on the solver, it may or may not make a convergence difference.</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;"><div id=":dr">
<br>
using -mat_fd_type_ds -snes_type test:<br>
<div class="im"><br>
Testing hand-coded Jacobian, if the ratio is<br>
O(1.e-8), the hand-coded Jacobian is probably correct.<br>
Run with -snes_test_display to show difference<br>
of hand-coded and finite difference Jacobian.<br>
</div> Norm of matrix ratio 9.17279e-09 difference 2.58883e-06<br>
Norm of matrix ratio 8.79304e-08 difference 2.59233e-05<br>
Norm of matrix ratio 8.87045e-08 difference 2.58875e-05<br></div></blockquote><div><br></div><div>This numbers are an order of magnitude smaller than before which indicates that the differencing is sensitive to choice of step.</div>
<div><br></div><div>In the best case, cond=1, this differencing can give you sqrt(machine_epsilon) which for double precision means you could get 7 or 8 digits correct. You lose precision when the conditioning of the operator degrades (because of the mesh, physics, or just poor relative scaling of fields) and if the step size is not chosen well. The default step size algorithm "wp" is cheaper to compute in parallel, "ds" is more robust.</div>
<div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;"><div id=":dr">
<br>
while using -ksp_rtol 1e-12 -pc_type lu -snes_monitor<br>
<br>
0 SNES Function norm 2.981436587230e+00<br>
1 SNES Function norm 1.487397936082e+00<br>
2 SNES Function norm 3.490136703474e-01<br>
3 SNES Function norm <a href="tel:1.3033970003">1.3033970003</a>44e-01<br>
4 SNES Function norm 3.181190030949e-02<br>
5 SNES Function norm 5.280858553562e-03<br>
6 SNES Function norm 3.168168082141e-04<br>
7 SNES Function norm <a href="tel:1.2177759654">1.2177759654</a>90e-06<br>
8 SNES Function norm <a href="tel:1.5593226240">1.5593226240</a>68e-11<br></div></blockquote><div><br></div><div>This looks like quadratic convergence so your Jacobian is probably correct.</div><div><br></div><div><div>
> is there a way to see/save/load into matlab the difference between the computed jacobian and the fd one, apart from output to screen?</div></div><div><br></div><div>Not as currently written, but it would be simple to add. If you are interested, look at src/snes/impls/test/snestest.c. You could add a runtime option -snes_test_display_matlab <filename> and then view to a Matlab file inside SNESSolve_Test(). Let us know if you get lost (preferably on petsc-dev).</div>
</div>