[petsc-users] Issues with the Variational Inequality solver

Thu Dec 3 15:00:07 CST 2015

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.000000

cA: 9.400282e-01

cB: 1.555428e-16

cC: 6.045961e-09

Case 1:

psiA: 8.472901e-01

psiC: 7.425271e-09

k: 2.000000

cA: 8.472901e-01

cB: 2.602870e-16

cC: 7.425270e-09

Case 2:

psiA: 8.337199e-01

psiC: 7.831614e-09

k: 2.000000

cA: 8.337199e-01

cB: 2.942675e-16

cC: 7.831613e-09

Case 3:

psiA: 8.268140e-01

psiC: 7.912911e-09

k: 2.000000

cA: 8.268140e-01

cB: 3.029178e-16

cC: 7.912910e-09

Case 4:

psiA: 9.852477e-01

psiC: 7.992223e-09

k: 2.000000

cA: 9.852477e-01

cB: 2.593282e-16

cC: 7.992222e-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.000000

cA: 1.318866e-16

cB: 3.097570e-08

cC: 6.045961e-09

Case 1:

psiA: 8.472901e-01

psiC: 7.425271e-09

k: 2.000000

cA: 4.624944e-17

cB: 3.015792e-08

cC: 7.425271e-09

Case 2:

psiA: 8.337199e-01

psiC: 7.831614e-09

k: 2.000000

cA: 1.733276e-16

cB: 3.079856e-08

cC: 7.831614e-09

Case 3:

psiA: 8.268140e-01

psiC: 7.912911e-09

k: 2.000000

cA: 1.721078e-16

cB: 3.061785e-08

cC: 7.912911e-09

Case 4:

psiA: 9.852477e-01

psiC: 7.992223e-09

k: 2.000000

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 <bsmith at mcs.anl.gov> wrote:

>
> > On Dec 2, 2015, at 1:56 PM, Justin Chang <jychang48 at gmail.com> 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 <bsmith at mcs.anl.gov> wrote:
> >
> > > On Dec 1, 2015, at 5:13 PM, Justin Chang <jychang48 at gmail.com> 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 <module>
> > >     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,&isequal);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
> > > <testVI.py><psiA_1><psiB_1>
> >
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151203/52f1b186/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: testVI.py
Type: text/x-python-script
Size: 3071 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151203/52f1b186/attachment-0001.bin>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: testA
Type: application/octet-stream
Size: 64 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151203/52f1b186/attachment-0004.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: testC
Type: application/octet-stream
Size: 64 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151203/52f1b186/attachment-0005.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: newton_output
Type: application/octet-stream
Size: 11308 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151203/52f1b186/attachment-0006.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: vi_output
Type: application/octet-stream
Size: 43977 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151203/52f1b186/attachment-0007.obj>