<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Dec 3, 2015 at 3:00 PM, Justin Chang <span dir="ltr"><<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Thanks Barry,<div><br></div><div>Next issue:</div></div></blockquote><div><br></div><div>Barry, when Justin described this problem to me, it sounded like there might be a bug in the active set method which</div><div>put in the lower value when it was supposed to be the upper value.</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Consider the following problem:</div><div><br></div><div>psiA = cA - 2*cB (1a)</div><div>psiC = cC + 2*cB (1b)</div><div>cC^2 = k*cA^2*cB (1c)</div><div><br></div><div>where psiA, psiC, and k are given. If I solve these five problems using the standard Newton method:</div><div><br></div><div><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">Case 0:</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiA: 9.400282e-01</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiC: 6.045961e-09</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>k: 2.000000</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cA: 9.400282e-01</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cB: 1.555428e-16</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cC: 6.045961e-09</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">Case 1:</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiA: 8.472901e-01</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiC: 7.425271e-09</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>k: 2.000000</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cA: 8.472901e-01</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cB: 2.602870e-16</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cC: 7.425270e-09</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">Case 2:</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiA: 8.337199e-01</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiC: 7.831614e-09</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>k: 2.000000</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cA: 8.337199e-01</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cB: 2.942675e-16</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cC: 7.831613e-09</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">Case 3:</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiA: 8.268140e-01</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiC: 7.912911e-09</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>k: 2.000000</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cA: 8.268140e-01</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cB: 3.029178e-16</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cC: 7.912910e-09</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">Case 4:</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiA: 9.852477e-01</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiC: 7.992223e-09</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>k: 2.000000</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cA: 9.852477e-01</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cB: 2.593282e-16</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">
</p><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cC: 7.992222e-09</p><br>These solutions are "correct", that is if I plug the c's back into equations (1a)-(1b), i get the original psi's. </div><div><br></div><div>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:</div><div><br></div><div><p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">Case 0:</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiA: 9.400282e-01</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiC: 6.045961e-09</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>k: 2.000000</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cA: 1.318866e-16</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cB: 3.097570e-08</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cC: 6.045961e-09</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">Case 1:</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiA: 8.472901e-01</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiC: 7.425271e-09</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>k: 2.000000</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cA: 4.624944e-17</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cB: 3.015792e-08</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cC: 7.425271e-09</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">Case 2:</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiA: 8.337199e-01</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiC: 7.831614e-09</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>k: 2.000000</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cA: 1.733276e-16</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cB: 3.079856e-08</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cC: 7.831614e-09</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">Case 3:</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiA: 8.268140e-01</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiC: 7.912911e-09</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>k: 2.000000</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cA: 1.721078e-16</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cB: 3.061785e-08</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cC: 7.912911e-09</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)">Case 4:</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiA: 9.852477e-01</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>psiC: 7.992223e-09</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>k: 2.000000</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cA: 2.666770e-16</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cB: 4.610822e-08</p>
<p style="margin:0px;font-size:12px;line-height:normal;font-family:Courier;color:rgb(76,47,45);background-color:rgb(223,219,196)"><span style="white-space:pre-wrap"> </span>cC: 7.992223e-09</p><br>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?</div><div><br></div><div>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:</div><div><br></div><div>python testVI.py testA testC 2 <0/1></div><div><br></div><div>where 0 is Newton, 1 is VI.</div><div><br></div><div>Know what's going on here? </div><div><br></div><div>Thanks,</div><div>Justin</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Dec 2, 2015 at 1:36 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span><br>
> On Dec 2, 2015, at 1:56 PM, Justin Chang <<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>> wrote:<br>
><br>
> Barry,<br>
><br>
> 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.<br>
><br>
> I haven't run the C debugger or done any of that yet, but could this be a plausible explanation for my error?<br>
<br>
</span> 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.<br>
<span><br>
> Originally I was thinking about creating/destroying SNES for each problem but I was wondering if that might affect performance.<br>
<br>
</span> Using SNESReset() is the way you should do it. Creating and destroying it each time should only be a tiny, immeasurably slower.<br>
<span><font color="#888888"><br>
Barry<br>
</font></span><div><div><br>
><br>
> Thanks,<br>
> Justin<br>
><br>
> On Tue, Dec 1, 2015 at 5:58 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
> > On Dec 1, 2015, at 5:13 PM, Justin Chang <<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>> wrote:<br>
> ><br>
> > Hi all,<br>
> ><br>
> > 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.<br>
> ><br>
> > First issue, when I run the program like this:<br>
> ><br>
> > python testVI.py psiA_1 psiB_1 2 1<br>
> ><br>
> > I get this error:<br>
> ><br>
> > Traceback (most recent call last):<br>
> > File "testVI.py", line 103, in <module><br>
> > snes.solve(None,X)<br>
> > File "PETSc/SNES.pyx", line 520, in petsc4py.PETSc.SNES.solve (src/petsc4py.PETSc.c:169677)<br>
> > petsc4py.PETSc.Error: error code 60<br>
> > [0] SNESSolve() line 3992 in /Users/justin/Software/petsc/src/snes/interface/snes.c<br>
> > [0] SNESSolve_VINEWTONRSLS() line 507 in /Users/justin/Software/petsc/src/snes/impls/vi/rs/virs.c<br>
> > [0] KSPSetOperators() line 544 in /Users/justin/Software/petsc/src/ksp/ksp/interface/itcreate.c<br>
> > [0] PCSetOperators() line 1170 in /Users/justin/Software/petsc/src/ksp/pc/interface/precon.c<br>
> > [0] Nonconforming object sizes<br>
> > [0] Cannot change local size of Amat after use old sizes 2 2 new sizes 3 3<br>
> ><br>
> > 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?<br>
><br>
> Justin,<br>
><br>
> 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.<br>
><br>
> In the rs solver we have<br>
><br>
> ierr = ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal);CHKERRQ(ierr);<br>
> if (!isequal) {<br>
> ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr);<br>
> }<br>
><br>
> 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.<br>
><br>
> 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.<br>
><br>
> Barry<br>
><br>
> ><br>
> > More issues to come :)<br>
> ><br>
> > Thanks,<br>
> > Justin<br>
> > <testVI.py><psiA_1><psiB_1><br>
><br>
><br>
<br>
</div></div></blockquote></div><br></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>