(Petsc SNES) Reevaluating Function right after evaluating the Jacobian

[Barry Smith]

On Dec 1, 2008, at 3:21 PM, Edson Tadeu wrote:

>  Barry,
>
> It's not that I'm changing the definition of the function... the  
> problem is that the function have some "ifs" on it. This comes  
> ultimately from the physics of the problem... there is some change  
> of dependency between the variables, like this:
>
> A is an unknown, B is an unknown, C is dependent (not on the system):
>
>    C = 1 - A - B + delta_C
>
> (where delta_C is also dependent from some other variables).
>
> But if C == 0 and delta_C == 0, then the dependency changes:
>
> A is an unknown, but B is now dependent:
>
>    B = 1 - A
>
> This is why, when this happens, the Jacobian row of B becomes  
> linearly dependent with the row of A. So, actually, B should be  
> taken out of the system! [See that there is a real change here,  
> because now the derivative dB/dA == -1, but before, A and B were  
> unrelated.]
>
> Well... this scheme was working before using full Newton steps... ;)

    Yes, but it may explain why Newton's method was converging slowly  
for you.
>
>
> But thanks for the tip, I'll also try to treat accordingly the null  
> space that appears with the linear dependence. Hmm... but by  
> removing the null space of J, I'll be effectively removing the row  
> corresponding to the B unknown... so, I would need to define delta_B  
> = -delta_A later, anyway, wouldn't I (or maybe PETSc would define it  
> automatically, given the basis vector [1 -1 0 0 0...] for the null  
> space)?

   If you know the null space then you can use KSPSetNullSpace() to  
set the null space just before the linear
solve.  It is unfortunately a bit of a hack to do that. In the end of  
your FormJacobian() do a SNESGetKSP() and then
create your null space object with MatNullSpaceCreate() and then  
attach it to the KSP with KSPSetNullSpace().
If you do this then YOU don't need to monkey with the changing the  
search direction manually.

    Barry

BTW: I don't totally know what I am talking about here.


>
>
>   Thanks,
> Edson
>
>
> On Mon, Dec 1, 2008 at 6:29 PM, Barry Smith <bsmith at mcs.anl.gov>  
> wrote:
>
>   The mathematician side of me cringes (actually every side of me  
> cringes) when I read this.
>
>   You are changing the definition of the function you are finding  
> the zero of WHILE you are
> finding the zero of it. Yuck. How do you know what you get at the  
> end is anything but garbage?
>
>  It seems to me that either
> 1) there is a bug in your Jacobian or
> 2) near (or at) the solution (when F(x) = 0) the Jacobian is singular
>
> For 1) you can run with -snes_ls_type test to do some simple tests  
> on the Jacobian.
>
> For 2) I submit that what you really want to do (instead of changing  
> the "Jacobian" to something
> that is not the Jacobian) is to properly solve the singular linear  
> system. How to do this
> depends on the nature of the null space of the Jacobian and if you  
> have a handle on what
> the null space is). In PETSc the simplest thing to try  
> (sequentially) is -pc_type lu -pc_factor_shift_nonzero
> and see what happens.
>
>   Barry
>
>
>
>
> On Dec 1, 2008, at 2:04 PM, Edson Tadeu wrote:
>
> Barry,
>
>  Thanks for the reply;
>
> On Mon, Dec 1, 2008 at 4:56 PM, Barry Smith <bsmith at mcs.anl.gov>  
> wrote:
>
> On Dec 1, 2008, at 12:28 PM, Edson Tadeu wrote:
>
>  Hi,
>
>  I'm solving non-linear problems with Line Search, and sometimes I  
> need to change an entire row of the Jacobian together with the  
> corresponding element in the Function, but I only know which rows/ 
> elements need to be changed when I'm computing the Jacobian. So,  
> what I'd like to know first is:
>
>  Maybe if you explain why you need to "change an entire row of the  
> Jacobian together with the corresponding element in the Function" we  
> could come up
> with some suggestions on how to proceed.
>
> It's because, depending on some conditions, I need to make some  
> change of variables in the original equations, otherwise the  
> Jacobian will have linear dependent rows. For example, I need to  
> tell that x(4) = 1.0-x(3) (instead of the original equation for  
> x(4)). I would do it by changing the Jacobian like this:
>
> J(4,3) = 1; J(4,4) = 1; all other columns of row 4 of J are zero.  
> F(4) = 0.
>
> This way, I would have delta_x(4) = -delta_x(3). I've solved this  
> now by directly changing the delta solution vector (delta_x) in the  
> pre check phase... it seems to be working, but it doesn't seem to be  
> the right way, I don't know...
>
>
>
>
>
>  1) Am I allowed to change the Function vector from inside  
> FormJacobian? (I've looked inside SNESSolve_LS, and it seems that  
> this is not allowed, because "fnorm" is evaluated before  
> SNESComputeJacobian, and if I change F inside FormJacobian, fnorm  
> would be wrong).
>
>  or
>
>  2) Am I allowed to change the Function vector from PreCheck routine?
>
>  If none of that is allowed, another solution that I thought was to  
> compute the Jacobian while computing the Function. The FAQ says that  
> "You are free to have your 'FormFunction' compute as much of the  
> Jacobian at that point as you like, keep the information in the user  
> context (the final argument to FormFunction and FormJacobian) and  
> then retreive the information in your FormJacobian() function". The  
> problem is that FormFunction is called many times more than  
> FormJacobian, and Jacobian calculation is slow, so:
>
>  3) Is there any way to know when I should really recompute the  
> Jacobian from inside FormFunction? (I don't think so, but...)
>
>  No. But if you do not use any line search at all then each function  
> evaluation will have an associated Jacobian so you can compute them
> together. (Of course you lose the globalization of the line search  
> so SNES may not converge).  You can turn off all line searches with - 
> snes_ls basic
> from the command line or from the code  
> SNESLineSearchSet(snes,SNESLineSeachNo,PETSC_NULL);
>
>  Barry
>
>  Ok, but the whole reason that I'm using SNES is to speed up my  
> computations by using the Line Search algorithms (or maybe Trust  
> Region). I had coded an algorithm before that manually computed  
> using full steps, and it were working, but were slow!
>
>  Thanks,
> Edson
>
>
>