October 2014 - petsc-dev - mailman3.cels.anl.gov

Re: [petsc-dev] Sanity Check for were initial conditions are set in example ex45
by Barry Smith 31 Oct '14

31 Oct '14

The initial guess function is called only once in KSPSetUp() (which is automatically called by KSPSolve() if you don’t call it). if (ksp->dmActive) { DMKSP kdm; ierr = DMGetDMKSP(ksp->dm,&kdm);CHKERRQ(ierr); if (kdm->ops->computeinitialguess && ksp->setupstage != KSP_SETUP_NEWRHS) { /* only computes initial guess the first time through */ ierr = (*kdm->ops->computeinitialguess)(ksp,ksp->vec_sol,kdm->initialguessctx);CHKERRQ(ierr); ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); It is not called inside of computing the right hand side. Barry > On Oct 31, 2014, at 1:29 PM, Alletto, John M <john.m.alletto(a)lmco.com> wrote: > > All, > > Please provide me with a sanity check of where the initial conditions get set in the Laplacian program example 45 > > Currently I assign the initial voltages to a 3-D array called barray in ComputeRHS, these very between 0 and 100; > > I initially believed initialization was occurring in ComputeRHS, but because my final output goes from 0 to 1, I am questioning my understanding of where the initial conditions go. > > Can you set me straight? > > > > Thanks > John

1 0

Sanity Check for were initial conditions are set in example ex45
by Alletto, John M 31 Oct '14

31 Oct '14

All, Please provide me with a sanity check of where the initial conditions get set in the Laplacian program example 45 Currently I assign the initial voltages to a 3-D array called barray in ComputeRHS, these very between 0 and 100; I initially believed initialization was occurring in ComputeRHS, but because my final output goes from 0 to 1, I am questioning my understanding of where the initial conditions go. Can you set me straight? Thanks John

1 0

Re: [petsc-dev] [petsc-users] gamg failure with petsc-dev
by Stephan Kramer 30 Oct '14

30 Oct '14

On 29/10/14 20:06, Jed Brown wrote: > Stephan Kramer <s.kramer(a)imperial.ac.uk> writes: > >> On 18/09/14 15:41, Mark Adams wrote: >>> Thanks for updating me Stephan, >>> >>> This is a common problem. I'm surprised we have not seen it more. >>> >>> Note, the problem with the solver dying when MatSetBlockSize is set >>> (correctly). Is some sort of bug. You have Stokes and use FieldSplit, >>> so GAMG is working on a 2D elasticity like problem. Is that correct? >>> And the equations are ordered node first: [v0_x, v0_y, v1_x, v1_y, >>> v2_x, ....] correct? >> >> Correct - although we're not using FieldSplit because of historic >> reasons (the code predates fieldsplit) - but yes, it's the velocity >> block K of a Stokes system (K G; G^T 0) with variable viscosity. >> >> I finally figured out why the convergence was so much worse if we do set >> the block size: It turns out to do with the scaling of the number I put >> on the diagonal for the eliminated DOFs associated with strong bcs. This >> problem has free slip bcs and I eliminate the DOF associated with the >> normal component by zeroing out the row and column and putting an >> arbitrary number on the diagonal (I chose 1.0). > > Are the slip surfaces always aligned with the coordinate axes or are you > rotating the vertex coordinate system to line up? If the latter, is > that rotation accounted for when defining the near-null space? The boundaries are non-aligned and the boundary conditions are rotated indeed. Both our physical coordinate space and the velocity components that we solve for are x,y aligned. After the assembly of the matrix (with "natural" bcs only) the dofs at the boundary nodes are rotated with a PTAP to align with the normal and tangential directions. After that the strong bc is applied with a call to MatZeroRowsColumns(). The current problem I'm looking at is in a cylindrical (annulus) domain with free slip on all sides. Thus we have a real null space containing the rotational mode, and a near null space containing all three of the usual modes. We have indeed applied the same rotation to all these null space modes > >> Normally the number doesn't really matter as the row is completely >> decoupled from the rest of the system. However in this case, when >> using gamg with blocks, it affects the strong coupling criterion for >> *both* velocity components associated with the boundary node. This >> criterion is scaled with respect to the 1-norm of the diagonal block, >> thus if the diagonal for the normal component is chosen too big it >> also declares all connections to the dof associated with the >> tangential component weak. In this way all the boundary nodes (the >> pairs of velocity components at these nodes) became isolated clusters >> and were disregarded at the coarse level. After choosing a smaller >> value for the diagonals the solve now converges fine in roughly the >> same number of iterations as when not setting the block size. >> >> The tricky thing here is the choice of the dummy diagonal number becomes >> very sensitive - with a typical gamg threshold of 0.01 it can't be more >> than an order out. Currently I use MatZeroRowsColumns() with a constant >> value for the diagonal. Because I have a variable viscosity it seems I >> can no longer do that. Now I can just go in and look at what the >> diagonal is for the other component at each boundary node and use that >> value instead (or at least a value that's smaller) - but it all seems a >> bit hackish to me. Perhaps we're following the wrong approach altogether >> and should be dealing with the bcs in an other way? Any thoughts would >> be more than welcome... > > I think we should improve the strength of connection measure. For > example, the (2 bs)×(2 bs) matrix > > A B > C D > > is strongly coupled if > > A^{-1} (A - B D^{-1} C) > > is far from the identity. It would normally be measured as > > sqrt(|| A^{-1} B D^{-1} C ||) > > In the case of scalars, this reduces to the standard measure > > sqrt(|bc / ad|). > > > Note that this still doesn't work in more exotic cases like unaligned > anisotropy with certain discretizations, but it should fix the problem > you're seeing. > That sounds like a nice solution as indeed it falls back to the scalar case when one of the two dofs is lifted. Cheers Stephan

1 0

Re: [petsc-dev] [petsc-users] gamg failure with petsc-dev
by Jed Brown 29 Oct '14

29 Oct '14

Stephan Kramer <s.kramer(a)imperial.ac.uk> writes: > On 18/09/14 15:41, Mark Adams wrote: >> Thanks for updating me Stephan, >> >> This is a common problem. I'm surprised we have not seen it more. >> >> Note, the problem with the solver dying when MatSetBlockSize is set >> (correctly). Is some sort of bug. You have Stokes and use FieldSplit, >> so GAMG is working on a 2D elasticity like problem. Is that correct? >> And the equations are ordered node first: [v0_x, v0_y, v1_x, v1_y, >> v2_x, ....] correct? > > Correct - although we're not using FieldSplit because of historic > reasons (the code predates fieldsplit) - but yes, it's the velocity > block K of a Stokes system (K G; G^T 0) with variable viscosity. > > I finally figured out why the convergence was so much worse if we do set > the block size: It turns out to do with the scaling of the number I put > on the diagonal for the eliminated DOFs associated with strong bcs. This > problem has free slip bcs and I eliminate the DOF associated with the > normal component by zeroing out the row and column and putting an > arbitrary number on the diagonal (I chose 1.0). Are the slip surfaces always aligned with the coordinate axes or are you rotating the vertex coordinate system to line up? If the latter, is that rotation accounted for when defining the near-null space? > Normally the number doesn't really matter as the row is completely > decoupled from the rest of the system. However in this case, when > using gamg with blocks, it affects the strong coupling criterion for > *both* velocity components associated with the boundary node. This > criterion is scaled with respect to the 1-norm of the diagonal block, > thus if the diagonal for the normal component is chosen too big it > also declares all connections to the dof associated with the > tangential component weak. In this way all the boundary nodes (the > pairs of velocity components at these nodes) became isolated clusters > and were disregarded at the coarse level. After choosing a smaller > value for the diagonals the solve now converges fine in roughly the > same number of iterations as when not setting the block size. > > The tricky thing here is the choice of the dummy diagonal number becomes > very sensitive - with a typical gamg threshold of 0.01 it can't be more > than an order out. Currently I use MatZeroRowsColumns() with a constant > value for the diagonal. Because I have a variable viscosity it seems I > can no longer do that. Now I can just go in and look at what the > diagonal is for the other component at each boundary node and use that > value instead (or at least a value that's smaller) - but it all seems a > bit hackish to me. Perhaps we're following the wrong approach altogether > and should be dealing with the bcs in an other way? Any thoughts would > be more than welcome... I think we should improve the strength of connection measure. For example, the (2 bs)×(2 bs) matrix A B C D is strongly coupled if A^{-1} (A - B D^{-1} C) is far from the identity. It would normally be measured as sqrt(|| A^{-1} B D^{-1} C ||) In the case of scalars, this reduces to the standard measure sqrt(|bc / ad|). Note that this still doesn't work in more exotic cases like unaligned anisotropy with certain discretizations, but it should fix the problem you're seeing.

1 0

Re: [petsc-dev] Backend independent VecGetArray for GPUs
by Karl Rupp 29 Oct '14

29 Oct '14

Hi, great, thanks for the feedback. I'll let you know once this is master. Best regards, Karli On 10/29/2014 08:16 PM, Lisandro Dalcin wrote: > On 27 October 2014 07:55, Karl Rupp <rupp(a)iue.tuwien.ac.at> wrote: >> >> thanks for the feedback! I'll clean up the feature branch and propagate the >> functionality to next tomorrow. >> > > Karl, Ashwin confirmed all this is working, so I'll wait for you to > merge your feature branch to master and then update petsc4py, please > let me know once done. > >

1 0

Re: [petsc-dev] Backend independent VecGetArray for GPUs
by Lisandro Dalcin 29 Oct '14

29 Oct '14

On 27 October 2014 07:55, Karl Rupp <rupp(a)iue.tuwien.ac.at> wrote: > > thanks for the feedback! I'll clean up the feature branch and propagate the > functionality to next tomorrow. > Karl, Ashwin confirmed all this is working, so I'll wait for you to merge your feature branch to master and then update petsc4py, please let me know once done. -- Lisandro Dalcin ============ Research Scientist Computer, Electrical and Mathematical Sciences & Engineering (CEMSE) Numerical Porous Media Center (NumPor) King Abdullah University of Science and Technology (KAUST) http://numpor.kaust.edu.sa/ 4700 King Abdullah University of Science and Technology al-Khawarizmi Bldg (Bldg 1), Office # 4332 Thuwal 23955-6900, Kingdom of Saudi Arabia http://www.kaust.edu.sa Office Phone: +966 12 808-0459

1 0

Re: [petsc-dev] [petsc-users] gamg failure with petsc-dev
by Stephan Kramer 29 Oct '14

29 Oct '14

On 18/09/14 15:41, Mark Adams wrote: > Thanks for updating me Stephan, > > This is a common problem. I'm surprised we have not seen it more. > > Note, the problem with the solver dying when MatSetBlockSize is set > (correctly). Is some sort of bug. You have Stokes and use FieldSplit, > so GAMG is working on a 2D elasticity like problem. Is that correct? > And the equations are ordered node first: [v0_x, v0_y, v1_x, v1_y, > v2_x, ....] correct? Correct - although we're not using FieldSplit because of historic reasons (the code predates fieldsplit) - but yes, it's the velocity block K of a Stokes system (K G; G^T 0) with variable viscosity. I finally figured out why the convergence was so much worse if we do set the block size: It turns out to do with the scaling of the number I put on the diagonal for the eliminated DOFs associated with strong bcs. This problem has free slip bcs and I eliminate the DOF associated with the normal component by zeroing out the row and column and putting an arbitrary number on the diagonal (I chose 1.0). Normally the number doesn't really matter as the row is completely decoupled from the rest of the system. However in this case, when using gamg with blocks, it affects the strong coupling criterion for *both* velocity components associated with the boundary node. This criterion is scaled with respect to the 1-norm of the diagonal block, thus if the diagonal for the normal component is chosen too big it also declares all connections to the dof associated with the tangential component weak. In this way all the boundary nodes (the pairs of velocity components at these nodes) became isolated clusters and were disregarded at the coarse level. After choosing a smaller value for the diagonals the solve now converges fine in roughly the same number of iterations as when not setting the block size. The tricky thing here is the choice of the dummy diagonal number becomes very sensitive - with a typical gamg threshold of 0.01 it can't be more than an order out. Currently I use MatZeroRowsColumns() with a constant value for the diagonal. Because I have a variable viscosity it seems I can no longer do that. Now I can just go in and look at what the diagonal is for the other component at each boundary node and use that value instead (or at least a value that's smaller) - but it all seems a bit hackish to me. Perhaps we're following the wrong approach altogether and should be dealing with the bcs in an other way? Any thoughts would be more than welcome... Cheers Stephan

1 0

Re: [petsc-dev] [petsc-users] unreliable AMG in PETSc
by Mark Adams 29 Oct '14

29 Oct '14

> > > Mark, > > How can we insure that that the smoother is SPD? Which are known to be. > In his case it was Richardson + SSOR on a bunch of processes. > > The problem is stability. Richardson + SSOR is symmetric and stable for an SPD matrix (I think) but we do not have a global SSOR, so we really have Jacobi/SSOR. The iteration matrix for the smoothers has spectral radius like || Ident - omega D^{-1}A|| for damped (omega) Jacobi. If A is an M-matrix 0 < || D^{-1}A ||_2 < 2 (barely) and so spectral radius of this iteration is less that 1.0 and it is stable (slow but stable). So our Richardson + SSOR should not blow up on an M-matrix (provably) and works pretty well in practice, but even a 9-point stencil (2D) is not a M-matrix. And I think it is hopeless to try to detect indefiniteness, instability in the "random" (often buggy) matrices that we get tossed to us. Just too expensive to do anything reliable. Symmetry is easy to test for, the largest eigenvalue is considered easy but can be tricky, as we have seen. I'm not sure if there are any shortcuts if you just want to know if a symmetric matrix does not have negative eigenvalues? (I'm not optimistic). This is all you need but even the relatively easy high eigenvalue can be hard to make (reasonably) robust and even it has not proof. Mark

1 0

Re: [petsc-dev] GAMG documentation very thin?
by Mark Adams 28 Oct '14

28 Oct '14

Yes, I will give the doc some love. Mark On Mon, Oct 27, 2014 at 5:55 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > The manual page for PCGAMG and the users manual seem very very thin > are how to use PGGAMG, no mention of various useful options and no seeming > mention of near null spaces etc that, as I understand it, care important > for getting good convergence. Remember this is a library, not a framework? > > Barry > > > Library, something you get and use. Framework, something you write a > proposal with the authors of and then they may have some junior person help > you use it until the grant runs low or they are working on another proposal. > > >

1 0

Re: [petsc-dev] [petsc-users] unreliable AMG in PETSc
by Jed Brown 27 Oct '14

27 Oct '14

Barry Smith <bsmith(a)mcs.anl.gov> writes: >> Actually, the problem is more subtle because if one wants an SPD >> multigrid cycle, they only need that the post-smoother is the transpose >> of the pre-smoother. > > Hmm. Is this true? This guy’s matrix causes grief with > indefiniteness but since it uses SSOR smoother I think that > post-smooth should have been the transpose of the pre smooth > automatically? Yes, S_pre = S_post = S is symmetric so S_pre^T = S_post trivially. But S appears to be indefinite. It's simpler to understand indefiniteness for symmetric matrices (like SSOR), but we shouldn't get overzealous and put in checks for whether S is SPD.

1 0