<div dir="ltr">Hi Dmitry,<div><br></div><div>Thanks so much, that's great! I've updated the pull request now to include your patch:</div><div><a href="https://github.com/dknez/libmesh/commit/247cc12f7536a5f848fd5c63765d97af307d9fa7">https://github.com/dknez/libmesh/commit/247cc12f7536a5f848fd5c63765d97af307d9fa7</a><br></div><div>We can update the code later, once MatReset() is available in PETSc.</div><div><br></div><div>Regarding the convergence issues: I think the main issue is that here we're using the penalty method for contact, which gives a non-smooth term in the PDE, and hence we don't get nice convergence with the nonlinear solver. In particular, the SNES solver terminates with "DIVERGED_LINE_SEARCH", but I don't think that's very surprising here. The simulation results seem to be good though (e.g. see the screenshots in the PR). I'd be interested in your thoughts on this though?</div><div><br></div><div>I'm not sure why the KSP isn't converging fully. I didn't put any thought into the choice of default solver, happy to change that (to CG or whatever). But I also used a direct solver for comparison, and SuperLU gives the same overall result that I get with GMRES (though MUMPS crashes with "Numerically singular matrix" for some reason).</div><div><br></div><div>Regarding the VI solver: That sounds very good, I'll be interested to try the overhauled VI solver out for contact problems, once it's available.</div><div class="gmail_extra"><div class="gmail_quote"><br></div><div class="gmail_quote">Thanks,</div><div class="gmail_quote">David</div><div class="gmail_quote"><br></div><div class="gmail_quote"><br></div><div class="gmail_quote"><br></div><div class="gmail_quote">On Wed, Feb 25, 2015 at 9:56 AM, Dmitry Karpeyev <span dir="ltr"><<a href="mailto:karpeev@mcs.anl.gov" target="_blank">karpeev@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr">The previous patch didn't clean out the Mat data structure thoroughly enough. Please try the attached patch.  It runs for me, but there is no convergence.  In particular, KSP doesn't seem to converge, although that's with the default GMRES+BJACOBI/ILU(0) options. So I don't think <span style="font-size:13.1999998092651px">resetting the Mat is the cause.</span><div><br><div>Let me know if this works for you. <span style="font-size:13.1999998092651px">MatReset() will be in the next release as well as, hopefully, an overhauled VI solver </span><span style="font-size:13.1999998092651px">that should be able to handle these contact problems.</span></div><span class=""><font color="#888888"><div><br></div><div>Dmitry.</div></font></span><div><div class="h5"><br><div class="gmail_quote">On Mon Feb 23 2015 at 10:17:19 AM David Knezevic <<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.com</a>> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">OK, sounds good. Let me know if I can help with the digging.</div></div></div><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><br></div><div class="gmail_quote">David</div></div></div><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><br></div><div class="gmail_quote"><br></div><div class="gmail_quote"><br></div><div class="gmail_quote">On Mon, Feb 23, 2015 at 11:10 AM, Dmitry Karpeyev <span dir="ltr"><<a href="mailto:karpeev@mcs.anl.gov" target="_blank">karpeev@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr">I just tried building against petsc/master, but there needs to be more work on libMesh before it can work with petsc/master: <br>the new VecLockPush()/Pop() stuff isn't respected by vector manipulation in libMesh.<br><div>I put a hack equivalent to MatReset() into your branch (patch attached, in case you want to take a look at it), </div><div>but it generates the same error in MatCreateColmap <span style="line-height:1.5;font-size:13.1999998092651px">that you reported earlier.  It's odd that it occurs </span></div><div><span style="line-height:1.5;font-size:13.1999998092651px">on the second nonlinear iteration.   </span><span style="line-height:1.5;font-size:13.1999998092651px">I'll have to dig a bit deeper  </span><span style="line-height:1.5;font-size:13.1999998092651px">to see what's going on.</span></div><span><font color="#888888"><div><span style="line-height:1.5;font-size:13.1999998092651px"><br></span></div><div><span style="line-height:1.5;font-size:13.1999998092651px">Dmitry.</span></div></font></span><div><div><div><br></div><br><div class="gmail_quote">On Mon Feb 23 2015 at 10:03:33 AM David Knezevic <<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.com</a>> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr">Hi Dmitry,<div><br></div><div>OK, good to hear we're seeing the same behavior for the example.</div><div><br></div><div>Regarding this comment:</div><div><br></div><div><br></div><div class="gmail_extra"><div class="gmail_quote"></div></div></div><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div><div>libMesh needs to change the way configure extracts PETSc information -- configuration data were moved:</div><div>conf --> lib/petsc-conf</div><div>${PETSC_ARCH}/conf --> ${PETSC_ARCH}/lib/petsc-conf</div><div><br></div><div>At one point I started looking at m4/petsc.m4, but that got put on the back burner.  For now making the relevant symlinks by hand lets you configure and build libMesh with petsc/master.</div></div></div></blockquote><div><br></div><div><br></div></div></div></div><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div>So do you suggest that the next step here is to build libmesh against petsc/master so that we can try your PETSc pull request that implements MatReset() to see if that gets this example working?</div></div></div></div><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div><br></div><div>David</div></div></div></div><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div><br></div><div><br></div><div><br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div><div><div><div class="gmail_quote">On Mon Feb 23 2015 at 9:15:44 AM David Knezevic <<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.com</a>> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr">Hi Dmitry,<div><br></div><div>Thanks very much for testing out the example.</div><div><br></div><div><span style="font-size:12.8000001907349px">examples/systems_of_equations/</span><span style="font-size:12.8000001907349px"><u></u>ex8 works fine for me in serial, but it fails for me if I run with more than 1 MPI process. Can you try it with, say, 2 or 4 MPI processes?</span><br></div><div><span style="font-size:12.8000001907349px"><br></span></div><div><span style="font-size:12.8000001907349px">If we need access to MatReset in libMesh to get this to work, I'll be happy to work on a libMesh pull request for that.</span></div><div><span style="font-size:12.8000001907349px"><br></span></div><div><span style="font-size:12.8000001907349px">David</span></div><div><span style="font-size:12.8000001907349px"><br></span></div></div><div class="gmail_extra"><br clear="all"><div><div><div dir="ltr"><span style="color:rgb(0,0,0);font-family:Calibri,sans-serif;font-size:12px"><pre cols="72">--</pre><pre cols="72">David J. Knezevic | CTO
Akselos | 17 Bay State Road | Boston, MA | 02215
Phone (office): <a href="tel:%2B1-857-265-2238" value="+18572652238" target="_blank">+1-857-265-2238</a><br>Phone (mobile): <a href="tel:%2B1-617-599-4755" value="+16175994755" target="_blank">+1-617-599-4755</a>
Web: <a href="http://www.akselos.com" target="_blank">http://www.akselos.com</a></pre></span></div></div></div></div><div class="gmail_extra">
<br><div class="gmail_quote">On Mon, Feb 23, 2015 at 10:08 AM, Dmitry Karpeyev <span dir="ltr"><<a href="mailto:karpeev@mcs.anl.gov" target="_blank">karpeev@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr">David,<br><br><div>What code are you running when you encounter this error?  I'm trying to reproduce it and</div><div>I tried examples/systems_of_equations/<u></u>ex8, but it ran for me, ostensibly to completion.</div><div><br></div><div>I have a small PETSc pull request that implements MatReset(), which passes a small PETSc test,</div><div>but libMesh needs some work to be able to build against petsc/master because of some recent</div><div>changes to PETSc. </div><span><font color="#888888"><div><br></div><div>Dmitry.</div></font></span><div><div><br><div class="gmail_quote">On Mon Feb 23 2015 at 7:17:06 AM David Knezevic <<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.com</a>> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr">Hi Barry, hi Dmitry,<div><br></div><div>I set the matrix to BAIJ and back to AIJ, and the code got a bit further. But I now run into the error pasted below (Note that I'm now using "--with-debugging=1"):</div><div><br></div><div><div>PETSC ERROR: --------------------- Error Message ------------------------------<u></u>------------------------------<u></u>--</div></div></div><div dir="ltr"><div><div>PETSC ERROR: Petsc has generated inconsistent data</div></div></div><div dir="ltr"><div><div>PETSC ERROR: MPIAIJ Matrix was assembled but is missing garray</div></div></div><div dir="ltr"><div><div>PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/<u></u>documentation/faq.html</a> for trouble shooting.</div><div>PETSC ERROR: Petsc Release Version 3.5.2, Sep, 08, 2014 </div><div>PETSC ERROR: ./example-dbg on a arch-linux2-c-debug named david-Lenovo by dknez Mon Feb 23 08:05:44 2015</div><div>PETSC ERROR: Configure options --with-shared-libraries=1 --with-debugging=1 --download-suitesparse=1 --download-parmetis=1 --download-blacs=1 --download-scalapack=1 --download-mumps=1 --download-metis --download-superlu_dist --prefix=/home/dknez/software/<u></u>libmesh_install/dbg_real/petsc --download-hypre</div><div>PETSC ERROR: #1 MatCreateColmap_MPIAIJ_<u></u>Private() line 361 in /home/dknez/software/petsc-3.<u></u>5.2/src/mat/impls/aij/mpi/<u></u>mpiaij.c</div><div>PETSC ERROR: #2 MatSetValues_MPIAIJ() line 538 in /home/dknez/software/petsc-3.<u></u>5.2/src/mat/impls/aij/mpi/<u></u>mpiaij.c</div><div>PETSC ERROR: #3 MatSetValues() line 1136 in /home/dknez/software/petsc-3.<u></u>5.2/src/mat/interface/matrix.c</div><div>PETSC ERROR: #4 add_matrix() line 765 in /home/dknez/software/libmesh-<u></u>src/src/numerics/petsc_matrix.<u></u>C</div><div>------------------------------<u></u>------------------------------<u></u>--------------</div></div><div><br></div><div>This occurs when I try to set some entries of the matrix. Do you have any suggestions on how I can resolve this?</div><div><br></div><div>Thanks!</div></div><div dir="ltr"><div>David</div></div><div dir="ltr"><div><br></div><div><br></div><div><br></div><div><br></div><div>On Sun, Feb 22, 2015 at 10:22 PM, Dmitry Karpeyev <span dir="ltr"><<a href="mailto:dkarpeev@gmail.com" target="_blank">dkarpeev@gmail.com</a>></span> wrote:<br></div><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><br><br><div class="gmail_quote"><span>On Sun Feb 22 2015 at 9:15:22 PM Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><br>
> On Feb 22, 2015, at 9:09 PM, David Knezevic <<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.com</a>> wrote:<br>
><br>
> Hi Dmitry,<br>
><br>
> Thanks for the suggestion. I tried MatSetType(mat,MATMPIAIJ) followed by MatXAIJSetPreallocation(...), but unfortunately this still gives me the same error as before: "nnz cannot be greater than row length: local row 168 value 24 rowlength 0".<br>
><br>
> I gather that the idea here is that MatSetType builds a new matrix object, and then I should be able to pre-allocate for that new matrix however I like, right? Was I supposed to clear the matrix object somehow before calling MatSetType? (I didn't do any sort of clear operation.)<br>
<br>
  If the type doesn't change then MatSetType() won't do anything. You can try setting the type to BAIJ and then setting the type back to AIJ. This may/should clear out the matrix.<br></blockquote></span><div>Ah, yes.  If the type is the same as before it does quit early, but changing the type and then back will clear out and rebuild the matrix.   We need </div><div><span style="font-size:13.1999998092651px">something like MatReset() to do the equivalent thing. </span></div><span><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<br>
><br>
> As I said earlier, I'll make a dbg PETSc build, so hopefully that will help shed some light on what's going wrong for me.<br></blockquote></span><div>I think it's always a good idea to have a dbg build of PETSc when you doing things like these. </div><span><font color="#888888"><div><br></div><div>Dmitry. </div></font></span><div><div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<br>
   Don't bother, what I suggested won't work.<br>
<br>
  Barry<br>
<br>
<br>
><br>
> Thanks,<br>
> David<br>
><br>
><br>
><br>
><br>
> On Sun, Feb 22, 2015 at 6:02 PM, Dmitry Karpeyev <<a href="mailto:dkarpeev@gmail.com" target="_blank">dkarpeev@gmail.com</a>> wrote:<br>
> David,<br>
> It might be easier to just rebuild the whole matrix from scratch: you would in effect be doing all that with disassembling and resetting the preallocation.<br>
> MatSetType(mat,MATMPIAIJ)<br>
> or<br>
> PetscObjectGetType((<u></u>PetscObjec<u></u>t)mat,&type);<br>
> MatSetType(mat,type);<br>
> followed by<br>
> MatXAIJSetPreallocation(...);<br>
> should do.<br>
> Dmitry.<br>
><br>
><br>
> On Sun Feb 22 2015 at 4:45:46 PM Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
>  Do not call for SeqAIJ matrix. Do not call before the first time you have preallocated and put entries in the matrix and done the MatAssemblyBegin/End()<br>
><br>
>   If it still crashes you'll need to try the debugger<br>
><br>
>   Barry<br>
><br>
> > On Feb 22, 2015, at 4:09 PM, David Knezevic <<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.com</a>> wrote:<br>
> ><br>
> > Hi Barry,<br>
> ><br>
> > Thanks for your help, much appreciated.<br>
> ><br>
> > I added a prototype for MatDisAssemble_MPIAIJ:<br>
> > PETSC_INTERN PetscErrorCode MatDisAssemble_MPIAIJ(Mat);<br>
> ><br>
> > and I added a call to MatDisAssemble_MPIAIJ before MatMPIAIJSetPreallocation. However, I get a segfault on the call to MatDisAssemble_MPIAIJ. The segfault occurs in both serial and parallel.<br>
> ><br>
> > FYI, I'm using Petsc 3.5.2, and I'm not using a non-debug build (though I could rebuild PETSc in debug mode if you think that would help figure out what's happening here).<br>
> ><br>
> > Thanks,<br>
> > David<br>
> ><br>
> ><br>
> ><br>
> > On Sun, Feb 22, 2015 at 1:13 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
> ><br>
> >   David,<br>
> ><br>
> >    This is an obscure little feature of MatMPIAIJ,   each time you change the sparsity pattern before you call the MatMPIAIJSetPreallocation you need to call  MatDisAssemble_MPIAIJ(Mat mat).    This is a private PETSc function so you need to provide your own prototype for it above the function you use it in.<br>
> ><br>
> >   Let us know if this resolves the problem.<br>
> ><br>
> >    Barry<br>
> ><br>
> > We never really intended that people would call MatMPIAIJSetPreallocation() AFTER they had already used the matrix.<br>
> ><br>
> ><br>
> > > On Feb 22, 2015, at 6:50 AM, David Knezevic <<a href="mailto:david.knezevic@akselos.com" target="_blank">david.knezevic@akselos.com</a>> wrote:<br>
> > ><br>
> > > Hi all,<br>
> > ><br>
> > > I've implemented a solver for a contact problem using SNES. The sparsity pattern of the jacobian matrix needs to change at each nonlinear iteration (because the elements which are in contact can change), so I tried to deal with this by calling MatSeqAIJSetPreallocation and MatMPIAIJSetPreallocation during each iteration in order to update the preallocation.<br>
> > ><br>
> > > This seems to work fine in serial, but with two or more MPI processes I run into the error "nnz cannot be greater than row length", e.g.:<br>
> > > nnz cannot be greater than row length: local row 528 value 12 rowlength 0<br>
> > ><br>
> > > This error is from the call to<br>
> > > MatSeqAIJSetPreallocation(b-><u></u>B<u></u>,o_nz,o_nnz); in MatMPIAIJSetPreallocation_<u></u>MPIA<u></u>IJ.<br>
> > ><br>
> > > Any guidance on what the problem might be would be most appreciated. For example, I was wondering if there is a problem with calling SetPreallocation on a matrix that has already been preallocated?<br>
> > ><br>
> > > Some notes:<br>
> > > - I'm using PETSc via libMesh<br>
> > > - The code that triggers this issue is available as a PR on the libMesh github repo, in case anyone is interested: <a href="https://github.com/libMesh/libmesh/pull/460/" target="_blank">https://github.com/libMesh/<u></u>lib<u></u>mesh/pull/460/</a><br>
> > > - I can try to make a minimal pure-PETSc example that reproduces this error, if that would be helpful.<br>
> > ><br>
> > > Many thanks,<br>
> > > David<br>
> > ><br>
> ><br>
> ><br>
><br>
><br>
<br>
</blockquote></div></div></div></div>
</blockquote></div><br></div></div></blockquote></div></div></div></div>
</blockquote></div><br></div></blockquote></div></div></div></div></div>
</blockquote></div></div></div></blockquote></div></div></div></div>
</blockquote></div><br></div></div></blockquote></div></div></div></div></div>
</blockquote></div><br></div></div>