[petsc-users] Issues with the Variational Inequality solver
Barry Smith
bsmith at mcs.anl.gov
Thu Dec 3 16:26:13 CST 2015
> On Dec 3, 2015, at 3:59 PM, Justin Chang <jychang48 at gmail.com> wrote:
>
> 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.
"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 <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
>
>
More information about the petsc-users
mailing list