On Mon, Oct 31, 2016 at 8:44 AM, Jed Brown <j...@jedbrown.org> wrote:
> Jeremy Theler <jer...@seamplex.com> writes: > > > Hi again > > > > I have been wokring on these issues. Long story short: it is about the > > ordering of the unknown fields in the vector. > > > > Long story: > > The physics is linear elastic problem, you can see it does work with LU > > over a simple cube (warp the displacements to see it does represent an > > elastic problem, E=200e3, nu=0.3): > > > > https://caeplex.com/demo/results.php?id=5817146bdb561 > > > > > > Say my three displacements (unknowns) are u,v,w. I can define the > > unknown vector as (is this called node-based ordering?) > > > > [u1 v1 w1 u2 v2 w2 ... un vn wn]^T > > > > Another option is (is this called unknown-based ordering?) > > > > [u1 u2 ... un v1 v2 ... vn w1 w2 ... wn]^T > > > > > > With lu/preonly the results are the same, although the stiffnes matrixes > > for each case are attached as PNGs. And of course, the near-nullspace > > vectors are different. So PCSetCoordinates() should work with one > > ordering and not with another one, an issue I did not take into > > consideration. > > > > After understanding Matt's point about the near nullspace (and reading > > some interesting comments from Jed on scicomp stackexchange) I did built > > my own vectors (I had to take a look at MatNullSpaceCreateRigidBody() > > because I found out by running the code the nullspace should be an > > orthonormal basis, it should say so in the docs). > > Where? > > "vecs - the vectors that span the null space (excluding the constant > vector); these vectors must be orthonormal." > > https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/ > MatNullSpaceCreate.html > > And if you run in debug mode (default), as you always should until you > are confident that your code is correct, MatNullSpaceCreate tests that > your vectors are orthonormal. > > > Now, there are some results I do not understand. I tried these six > > combinations: > > > > order near-nullspace iterations norm > > ----- -------------- ---------- ---- > > unknown explicit 10 1.6e-6 > > unknown PCSetCoordinates 15 1.7e-7 > > unknown none 15 2.4e-7 > > node explicit fails with error -11 > > node PCSetCoordinates fails with error -11 > > node none 13 3.8e-7 > > Did you set a block size for the "node-based" orderings? Are you sure > the above is labeled correctly? Anyway, PCSetCoordinates uses > "node-based" ordering. Implementation performance will generally be > better with node-based ordering -- it has better memory streaming and > cache behavior. > > The AIJ matrix format will also automatically do an "inode" optimization > to reduce memory bandwidth and enable block smoothing (default > configuration uses SOR smoothing). You can use -mat_no_inode to try > turning that off. > > > Error -11 is > > PETSc's linear solver did not converge with reason > > 'DIVERGED_PCSETUP_FAILED' (-11) > > Isn't there an actual error message? > > > Any explanation (for dumbs)? > > Another thing to take into account: I am setting the dirichlet BCs with > > MatZeroRows(), but I am not updating the columns to keep symmetry. Can > > this pose a problem for GAMG? > > Usually minor, but it is better to maintain symmetry. > If the boundary values are not zero, no way to maintain symmetry unless we reduce the extra part of the matrix. Not updating the columns is better in this situation. Fande,
