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

Justin Chang jychang48 at gmail.com
Thu Dec 3 15:59:21 CST 2015 

Sorry forgot to hit reply all

---------- Forwarded message ----------
From: Justin Chang <jychang48 at gmail.com>
Date: Thu, Dec 3, 2015 at 2:40 PM
Subject: Re: [petsc-users] Issues with the Variational Inequality solver
To: Matthew Knepley <knepley at 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.000000

cA: 6.297272e-01

cB: 4.616557e-16

cC: 1.705043e-08

Case 1:

psiA: 7.570078e-01

psiC: 1.754758e-08

k: 2.000000

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

cA: -7.684746e-17

cB: 5.362599e-09

cC: 4.616557e-16

Case 1:

psiA: 7.570078e-01

psiC: 1.754758e-08

k: 2.000000

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.

??

Justin

On Thu, Dec 3, 2015 at 2:26 PM, Matthew Knepley <knepley at gmail.com> wrote:

> On Thu, Dec 3, 2015 at 3:00 PM, Justin Chang <jychang48 at 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.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>
>>> >
>>> >
>>>
>>>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151203/3912cf61/attachment-0001.html>

More information about the petsc-users mailing list