Re: [petsc-users] Issues with the Variational Inequality solver

[Barry Smith]

> On Dec 3, 2015, at 3:59 PM, Justin Chang <jychan...@gmail.com> wrote:
> 
> Sorry forgot to hit reply all 
> 
> -- Forwarded message --
> From: Justin Chang <jychan...@gmail.com>
> Date: Thu, Dec 3, 2015 at 2:40 PM
> Subject: Re: [petsc-users] Issues with the Variational Inequality solver
> To: Matthew Knepley <knep...@gmail.com>
> 
> 
> Matt, for these problems I set an upper bound of 10, so I am not sure that's 
> related to the problem at hand. However, that is the *next* issue I will 
> bring up if this current one is resolved ;)
> 
> Another thing I found out...
> 
> The five case studies I listed used no preconditioner. When I use the jacobi 
> preconditioner, the Newton and VI solutions are identical.
> 
> Whereas these problems now fail:
> 
> with Newton:
> Case 0:
>   psiA: 6.297272e-01
>   psiC: 1.705043e-08
>   k: 2.00
>   cA: 6.297272e-01
>   cB: 4.616557e-16
>   cC: 1.705043e-08
> Case 1:
>   psiA: 7.570078e-01
>   psiC: 1.754758e-08
>   k: 2.00
>   cA: 7.570078e-01
>   cB: 4.067561e-16
>   cC: 1.754758e-08
> 
> with VI:
> 
> Case 0:
>   psiA: 6.297272e-01
>   psiC: 1.705043e-08
>   k: 2.00
>   cA: -7.684746e-17
>   cB: 5.362599e-09
>   cC: 4.616557e-16
> Case 1:
>   psiA: 7.570078e-01
>   psiC: 1.754758e-08
>   k: 2.00
>   cA: 3.781787e-16
>   cB: 1.152521e-08
>   cC: 4.067561e-16
> 
> For the first set of problems, Jacobi works but no preconditioner fails.
> For the second set of problems, Jacobi fails but no preconditioner works.

   "Fails" is not very informative. What happens if you use LU? When debugging 
a nonlinear solver always use a direct solver for the linear solver (and make 
sure the direct solver actually worked).


> 
> ??
> 
> Justin
> 
> On Thu, Dec 3, 2015 at 2:26 PM, Matthew Knepley <knep...@gmail.com> wrote:
> On Thu, Dec 3, 2015 at 3:00 PM, Justin Chang <jychan...@gmail.com> wrote:
> Thanks Barry,
> 
> Next issue:
> 
> Barry, when Justin described this problem to me, it sounded like there might 
> be a bug in the active set method which
> put in the lower value when it was supposed to be the upper value.
> 
>Matt
>  
> Consider the following problem:
> 
> psiA = cA - 2*cB (1a)
> psiC = cC + 2*cB(1b)
> cC^2 = k*cA^2*cB  (1c)
> 
> where psiA, psiC, and k are given. If I solve these five problems using the 
> standard Newton method:
> 
> Case 0:
>   psiA: 9.400282e-01
>   psiC: 6.045961e-09
>   k: 2.00
>   cA: 9.400282e-01
>   cB: 1.555428e-16
>   cC: 6.045961e-09
> Case 1:
>   psiA: 8.472901e-01
>   psiC: 7.425271e-09
>   k: 2.00
>   cA: 8.472901e-01
>   cB: 2.602870e-16
>   cC: 7.425270e-09
> Case 2:
>   psiA: 8.337199e-01
>   psiC: 7.831614e-09
>   k: 2.00
>   cA: 8.337199e-01
>   cB: 2.942675e-16
>   cC: 7.831613e-09
> Case 3:
>   psiA: 8.268140e-01
>   psiC: 7.912911e-09
>   k: 2.00
>   cA: 8.268140e-01
>   cB: 3.029178e-16
>   cC: 7.912910e-09
> Case 4:
>   psiA: 9.852477e-01
>   psiC: 7.992223e-09
>   k: 2.00
>   cA: 9.852477e-01
>   cB: 2.593282e-16
>   cC: 7.99e-09
> 
> These solutions are "correct", that is if I plug the c's back into equations 
> (1a)-(1b), i get the original psi's. 
> 
> Now suppose I use the Variational Inequality such that cA/B/C >= 0, I would 
> expect to get the exact same solution (since the c's would be non-negative 
> regardless). But instead I get these solutions:
> 
> Case 0:
>   psiA: 9.400282e-01
>   psiC: 6.045961e-09
>   k: 2.00
>   cA: 1.318866e-16
>   cB: 3.097570e-08
>   cC: 6.045961e-09
> Case 1:
>   psiA: 8.472901e-01
>   psiC: 7.425271e-09
>   k: 2.00
>   cA: 4.624944e-17
>   cB: 3.015792e-08
>   cC: 7.425271e-09
> Case 2:
>   psiA: 8.337199e-01
>   psiC: 7.831614e-09
>   k: 2.00
>   cA: 1.733276e-16
>   cB: 3.079856e-08
>   cC: 7.831614e-09
> Case 3:
>   psiA: 8.268140e-01
>   psiC: 7.912911e-09
>   k: 2.00
>   cA: 1.721078e-16
>   cB: 3.061785e-08
>   cC: 7.912911e-09
> Case 4:
>   psiA: 9.852477e-01
>   psiC: 7.992223e-09
>   k: 2.00
>   cA: 2.666770e-16
>   cB: 4.610822e-08
>   cC: 7.992223e-09
> 
> Basically, the cA's shoot down to zero and the cB's are slightly greater 

Re: [petsc-users] Issues with the Variational Inequality solver

[Justin Chang]

Thanks Barry,

Next issue:

Consider the following problem:

psiA = cA - 2*cB (1a)
psiC = cC + 2*cB(1b)
cC^2 = k*cA^2*cB  (1c)

where psiA, psiC, and k are given. If I solve these five problems using the
standard Newton method:

Case 0:

psiA: 9.400282e-01

psiC: 6.045961e-09

k: 2.00

cA: 9.400282e-01

cB: 1.555428e-16

cC: 6.045961e-09

Case 1:

psiA: 8.472901e-01

psiC: 7.425271e-09

k: 2.00

cA: 8.472901e-01

cB: 2.602870e-16

cC: 7.425270e-09

Case 2:

psiA: 8.337199e-01

psiC: 7.831614e-09

k: 2.00

cA: 8.337199e-01

cB: 2.942675e-16

cC: 7.831613e-09

Case 3:

psiA: 8.268140e-01

psiC: 7.912911e-09

k: 2.00

cA: 8.268140e-01

cB: 3.029178e-16

cC: 7.912910e-09

Case 4:

psiA: 9.852477e-01

psiC: 7.992223e-09

k: 2.00

cA: 9.852477e-01

cB: 2.593282e-16

cC: 7.99e-09

These solutions are "correct", that is if I plug the c's back into
equations (1a)-(1b), i get the original psi's.

Now suppose I use the Variational Inequality such that cA/B/C >= 0, I would
expect to get the exact same solution (since the c's would be non-negative
regardless). But instead I get these solutions:

Case 0:

psiA: 9.400282e-01

psiC: 6.045961e-09

k: 2.00

cA: 1.318866e-16

cB: 3.097570e-08

cC: 6.045961e-09

Case 1:

psiA: 8.472901e-01

psiC: 7.425271e-09

k: 2.00

cA: 4.624944e-17

cB: 3.015792e-08

cC: 7.425271e-09

Case 2:

psiA: 8.337199e-01

psiC: 7.831614e-09

k: 2.00

cA: 1.733276e-16

cB: 3.079856e-08

cC: 7.831614e-09

Case 3:

psiA: 8.268140e-01

psiC: 7.912911e-09

k: 2.00

cA: 1.721078e-16

cB: 3.061785e-08

cC: 7.912911e-09

Case 4:

psiA: 9.852477e-01

psiC: 7.992223e-09

k: 2.00

cA: 2.666770e-16

cB: 4.610822e-08

cC: 7.992223e-09

Basically, the cA's shoot down to zero and the cB's are slightly greater
than zero. When I plug the c's into equations (1a) - (1b), I do not get the
correct solution at all. I would expect discrepancies if my c's were
originally negative, but for these five problems, shouldn't VI give the
same answer as the Newton's method?

Attached is the petsc4py code, the two input files containing psiA and
psiC, and the outputs from '-snes_monitor -snes_view
-ksp_monitor_true_residual'. Run as:

python testVI.py testA testC 2 <0/1>

where 0 is Newton, 1 is VI.

Know what's going on here?

Thanks,
Justin

On Wed, Dec 2, 2015 at 1:36 PM, Barry Smith  wrote:

>
> > On Dec 2, 2015, at 1:56 PM, Justin Chang  wrote:
> >
> > Barry,
> >
> > So I do not believe there is any problem with the ISEqual() per se,
> because what I am doing is I am solving 5151 different VI problem using the
> same SNES/KSP/PC/Mat. That is, I am not "resetting" these data structures
> once I am done with one problem and move on to the next. If I ran each case
> individually, I get no error; the error comes when I run these problems
> sequentially without resetting the SNES after each problem.
> >
> > I haven't run the C debugger or done any of that yet, but could this be
> a plausible explanation for my error?
>
>This is exactly the issue. Independent of VIs if you change the size of
> a problem you pass to SNES you need to call SNESReset() between each call
> of a different sizes.
>
> > Originally I was thinking about creating/destroying SNES for each
> problem but I was wondering if that might affect performance.
>
>Using SNESReset() is the way you should do it. Creating and destroying
> it each time should only be a tiny, immeasurably slower.
>
>   Barry
>
> >
> > Thanks,
> > Justin
> >
> > On Tue, Dec 1, 2015 at 5:58 PM, Barry Smith  wrote:
> >
> > > On Dec 1, 2015, at 5:13 PM, Justin Chang  wrote:
> > >
> > > Hi all,
> > >
> > > So I am running into some issues with the SNESVINEWTONRSLS solver. I
> originally had this done in firedrake, but have ported it to petsc4py.
> Attached is the code as well as the two required input files.
> > >
> > > First issue, when I run the program like this:
> > >
> > > python testVI.py psiA_1 psiB_1 2 1
> > >
> > > I get this error:
> > >
> > > Traceback (most recent call last):
> > >   File "testVI.py", line 103, in 
> > > snes.solve(None,X)
> > >   File "PETSc/SNES.pyx", line 520, in petsc4py.PETSc.SNES.solve
> (src/petsc4py.PETSc.c:169677)
> > > petsc4py.PETSc.Error: error code 60
> > > [0] SNESSolve() line 3992 in
> /Users/justin/Software/petsc/src/snes/interface/snes.c
> > > [0] SNESSolve_VINEWTONRSLS() line 507 in
> /Users/justin/Software/petsc/src/snes/impls/vi/rs/virs.c
> > > [0] KSPSetOperators() line 544 in
> /Users/justin/Software/petsc/src/ksp/ksp/interface/itcreate.c
> > > [0] PCSetOperators() line 1170 in
> /Users/justin/Software/petsc/src/ksp/pc/interface/precon.c
> > > [0] Nonconforming object sizes
> > > [0] Cannot change local size of Amat after use old sizes 2 2 new sizes
> 3 3
> > >
> > > No such error occurred when this was ran in firedrake, though I did
> notice that -snes_view did indicate that some of the iterations 

Re: [petsc-users] Issues with the Variational Inequality solver

[Justin Chang]

Barry,

So I do not believe there is any problem with the ISEqual() per se, because
what I am doing is I am solving 5151 different VI problem using the same
SNES/KSP/PC/Mat. That is, I am not "resetting" these data structures once I
am done with one problem and move on to the next. If I ran each case
individually, I get no error; the error comes when I run these problems
sequentially without resetting the SNES after each problem.

I haven't run the C debugger or done any of that yet, but could this be a
plausible explanation for my error? Originally I was thinking about
creating/destroying SNES for each problem but I was wondering if that might
affect performance.

Thanks,
Justin

On Tue, Dec 1, 2015 at 5:58 PM, Barry Smith  wrote:

>
> > On Dec 1, 2015, at 5:13 PM, Justin Chang  wrote:
> >
> > Hi all,
> >
> > So I am running into some issues with the SNESVINEWTONRSLS solver. I
> originally had this done in firedrake, but have ported it to petsc4py.
> Attached is the code as well as the two required input files.
> >
> > First issue, when I run the program like this:
> >
> > python testVI.py psiA_1 psiB_1 2 1
> >
> > I get this error:
> >
> > Traceback (most recent call last):
> >   File "testVI.py", line 103, in 
> > snes.solve(None,X)
> >   File "PETSc/SNES.pyx", line 520, in petsc4py.PETSc.SNES.solve
> (src/petsc4py.PETSc.c:169677)
> > petsc4py.PETSc.Error: error code 60
> > [0] SNESSolve() line 3992 in
> /Users/justin/Software/petsc/src/snes/interface/snes.c
> > [0] SNESSolve_VINEWTONRSLS() line 507 in
> /Users/justin/Software/petsc/src/snes/impls/vi/rs/virs.c
> > [0] KSPSetOperators() line 544 in
> /Users/justin/Software/petsc/src/ksp/ksp/interface/itcreate.c
> > [0] PCSetOperators() line 1170 in
> /Users/justin/Software/petsc/src/ksp/pc/interface/precon.c
> > [0] Nonconforming object sizes
> > [0] Cannot change local size of Amat after use old sizes 2 2 new sizes 3
> 3
> >
> > No such error occurred when this was ran in firedrake, though I did
> notice that -snes_view did indicate that some of the iterations had 2x2
> instead of 3x3 matrices. Why is this happening? I thought I should be
> solving a 3x3 system each time so why would a 2x2 ever arise?
>
>   Justin,
>
>For VI's the rs (reduced space solver) solves the linearized problem on
> a subset of the variables, thus the size of the linear system may change
> from iteration of Newton to the next iteration of Newton.
>
>In the rs solver we have
>
> ierr = ISEqual(vi->IS_inact_prev,vi->IS_inact,);CHKERRQ(ierr);
> if (!isequal) {
>   ierr =
> SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr);
> }
>
> so what __should__ happen is that if the size of the reduced space changes
> KSPReset() is called on the KSP allowing the KSP/PC to handle a different
> sized system then before. But why doesn't it work? Why isn't KSPReset()
> called? Somehow the ISEqual is not working.
>
>Can you run the C debugger and put a breakpoint in the ISEqual() then
> check the inputs and see if it is correctly setting the isequal to false
> when it should?  Each time the vi->IS_inact changes the KSPReset() should
> get called. You can check this in the debugger.
>
>Barry
>
> >
> > More issues to come :)
> >
> > Thanks,
> > Justin
> > 
>
>

Re: [petsc-users] Issues with the Variational Inequality solver

[Barry Smith]

> On Dec 2, 2015, at 1:56 PM, Justin Chang  wrote:
> 
> Barry,
> 
> So I do not believe there is any problem with the ISEqual() per se, because 
> what I am doing is I am solving 5151 different VI problem using the same 
> SNES/KSP/PC/Mat. That is, I am not "resetting" these data structures once I 
> am done with one problem and move on to the next. If I ran each case 
> individually, I get no error; the error comes when I run these problems 
> sequentially without resetting the SNES after each problem.
> 
> I haven't run the C debugger or done any of that yet, but could this be a 
> plausible explanation for my error?

   This is exactly the issue. Independent of VIs if you change the size of a 
problem you pass to SNES you need to call SNESReset() between each call of a 
different sizes. 

> Originally I was thinking about creating/destroying SNES for each problem but 
> I was wondering if that might affect performance.

   Using SNESReset() is the way you should do it. Creating and destroying it 
each time should only be a tiny, immeasurably slower.

  Barry

> 
> Thanks,
> Justin
> 
> On Tue, Dec 1, 2015 at 5:58 PM, Barry Smith  wrote:
> 
> > On Dec 1, 2015, at 5:13 PM, Justin Chang  wrote:
> >
> > Hi all,
> >
> > So I am running into some issues with the SNESVINEWTONRSLS solver. I 
> > originally had this done in firedrake, but have ported it to petsc4py. 
> > Attached is the code as well as the two required input files.
> >
> > First issue, when I run the program like this:
> >
> > python testVI.py psiA_1 psiB_1 2 1
> >
> > I get this error:
> >
> > Traceback (most recent call last):
> >   File "testVI.py", line 103, in 
> > snes.solve(None,X)
> >   File "PETSc/SNES.pyx", line 520, in petsc4py.PETSc.SNES.solve 
> > (src/petsc4py.PETSc.c:169677)
> > petsc4py.PETSc.Error: error code 60
> > [0] SNESSolve() line 3992 in 
> > /Users/justin/Software/petsc/src/snes/interface/snes.c
> > [0] SNESSolve_VINEWTONRSLS() line 507 in 
> > /Users/justin/Software/petsc/src/snes/impls/vi/rs/virs.c
> > [0] KSPSetOperators() line 544 in 
> > /Users/justin/Software/petsc/src/ksp/ksp/interface/itcreate.c
> > [0] PCSetOperators() line 1170 in 
> > /Users/justin/Software/petsc/src/ksp/pc/interface/precon.c
> > [0] Nonconforming object sizes
> > [0] Cannot change local size of Amat after use old sizes 2 2 new sizes 3 3
> >
> > No such error occurred when this was ran in firedrake, though I did notice 
> > that -snes_view did indicate that some of the iterations had 2x2 instead of 
> > 3x3 matrices. Why is this happening? I thought I should be solving a 3x3 
> > system each time so why would a 2x2 ever arise?
> 
>   Justin,
> 
>For VI's the rs (reduced space solver) solves the linearized problem on a 
> subset of the variables, thus the size of the linear system may change from 
> iteration of Newton to the next iteration of Newton.
> 
>In the rs solver we have
> 
> ierr = ISEqual(vi->IS_inact_prev,vi->IS_inact,);CHKERRQ(ierr);
> if (!isequal) {
>   ierr = 
> SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr);
> }
> 
> so what __should__ happen is that if the size of the reduced space changes 
> KSPReset() is called on the KSP allowing the KSP/PC to handle a different 
> sized system then before. But why doesn't it work? Why isn't KSPReset() 
> called? Somehow the ISEqual is not working.
> 
>Can you run the C debugger and put a breakpoint in the ISEqual() then 
> check the inputs and see if it is correctly setting the isequal to false when 
> it should?  Each time the vi->IS_inact changes the KSPReset() should get 
> called. You can check this in the debugger.
> 
>Barry
> 
> >
> > More issues to come :)
> >
> > Thanks,
> > Justin
> > 
> 
> 

Re: [petsc-users] Issues with the Variational Inequality solver

[Barry Smith]

> On Dec 1, 2015, at 5:13 PM, Justin Chang  wrote:
> 
> Hi all,
> 
> So I am running into some issues with the SNESVINEWTONRSLS solver. I 
> originally had this done in firedrake, but have ported it to petsc4py. 
> Attached is the code as well as the two required input files.
> 
> First issue, when I run the program like this:
> 
> python testVI.py psiA_1 psiB_1 2 1
> 
> I get this error:
> 
> Traceback (most recent call last):
>   File "testVI.py", line 103, in 
> snes.solve(None,X)
>   File "PETSc/SNES.pyx", line 520, in petsc4py.PETSc.SNES.solve 
> (src/petsc4py.PETSc.c:169677)
> petsc4py.PETSc.Error: error code 60
> [0] SNESSolve() line 3992 in 
> /Users/justin/Software/petsc/src/snes/interface/snes.c
> [0] SNESSolve_VINEWTONRSLS() line 507 in 
> /Users/justin/Software/petsc/src/snes/impls/vi/rs/virs.c
> [0] KSPSetOperators() line 544 in 
> /Users/justin/Software/petsc/src/ksp/ksp/interface/itcreate.c
> [0] PCSetOperators() line 1170 in 
> /Users/justin/Software/petsc/src/ksp/pc/interface/precon.c
> [0] Nonconforming object sizes
> [0] Cannot change local size of Amat after use old sizes 2 2 new sizes 3 3
> 
> No such error occurred when this was ran in firedrake, though I did notice 
> that -snes_view did indicate that some of the iterations had 2x2 instead of 
> 3x3 matrices. Why is this happening? I thought I should be solving a 3x3 
> system each time so why would a 2x2 ever arise?

  Justin,

   For VI's the rs (reduced space solver) solves the linearized problem on a 
subset of the variables, thus the size of the linear system may change from 
iteration of Newton to the next iteration of Newton. 

   In the rs solver we have 

ierr = ISEqual(vi->IS_inact_prev,vi->IS_inact,);CHKERRQ(ierr);
if (!isequal) {
  ierr = 
SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr);
}

so what __should__ happen is that if the size of the reduced space changes 
KSPReset() is called on the KSP allowing the KSP/PC to handle a different sized 
system then before. But why doesn't it work? Why isn't KSPReset() called? 
Somehow the ISEqual is not working. 

   Can you run the C debugger and put a breakpoint in the ISEqual() then check 
the inputs and see if it is correctly setting the isequal to false when it 
should?  Each time the vi->IS_inact changes the KSPReset() should get called. 
You can check this in the debugger.

   Barry

> 
> More issues to come :)
> 
> Thanks,
> Justin
>