[petsc-users] Issues with the Variational Inequality solver

Thu Dec 3 15:26:44 CST 2015

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/7d3c0fe2/attachment.html>