Re: [firedrake] Nullspace & VectorSpaceBasis
Four comments: 1) If I understand correctly, PETSc insists that the diagonal entries of the sparse matrix are filled in. Usually, the LHS form will cause this to happen; however, that obviously isn't the case for your form! 2) It should be possible to build the nullspace you want, but I'm clueless on the actual syntax. This wouldn't get around (1) without further fixes from our end. 3) Could the equations given not be included in the weak form for the rest of the problem? 4) If not, there might be a hacky "add-a-weighted-mass-matrix-to-the-LHS" fix you could use?.... Andrew On 11 December 2014 at 12:46, Anna Kalogirou <a.kalogirou@leeds.ac.uk> wrote:
Dear all,
I am trying to solve a simple problem (linear water wave equations), where I have the Laplace equation in a square domain with Neumann boundary conditions on the three sides and on the top boundary two equations are satisfied.
I have the following code: -------------------------------------- V = FunctionSpace(mesh,"CG",1)
eta0 = Function(V).interpolate(Expression("cos(2*pi*x[0])")) phi0 = Function(V) eta1 = Function(V) phi1 = Function(V) eta = TrialFunction(V) phi = TrialFunction(V) v = TestFunction(V)
aeta = v*eta*ds(4) Leta = v*eta0*ds(4) + dt*inner(grad(v),grad(phi0))*dx
eta_problem = LinearVariationalProblem(aeta,Leta,eta1) eta_solver = LinearVariationalSolver(eta_problem)
aphi = v*phi*ds(4) Lphi = (v*phi0 - g*dt*v*eta1)*ds(4)
phi_problem = LinearVariationalProblem(aphi,Lphi,phi1) phi_solver = LinearVariationalSolver(phi_problem) --------------------------------------
The error I get tells me that Object is in wrong state, and Matrix is missing diagonal entry 39.
I thought this had to do with the problem being singular, so I tried to build the nullspace but I don't know how to define the VectorSpaceBasis. Any help would be appreciated.
Thank you very much.
Regards, Anna.
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Thanks Andrew. About your comment 1), how can I fix this? The two equations valid on one of the boundaries can actually be combined into one by paying the price of a higher time derivative. This is not a problem for this simple system, so if I use aphi = v*phi*ds(4) + g*dt**2*inner(grad(v),grad(phi))*dx Lphi = v*(2*phi1-phi0)*ds(4) phi_problem = LinearVariationalProblem(aphi,Lphi,phi2) phi_solver = LinearVariationalSolver(phi_problem) is equivalent to what I had before, and this time the code works and gives a good result (compared to the exact solution). However, I need to understand why the other case doesn't work because that will later be generalised to become nonlinear. Anna. On 11/12/14 13:04, Andrew McRae wrote:
Four comments:
1) If I understand correctly, PETSc insists that the diagonal entries of the sparse matrix are filled in. Usually, the LHS form will cause this to happen; however, that obviously isn't the case for your form!
2) It should be possible to build the nullspace you want, but I'm clueless on the actual syntax. This wouldn't get around (1) without further fixes from our end.
3) Could the equations given not be included in the weak form for the rest of the problem?
4) If not, there might be a hacky "add-a-weighted-mass-matrix-to-the-LHS" fix you could use?....
Andrew
On 11 December 2014 at 12:46, Anna Kalogirou <a.kalogirou@leeds.ac.uk <mailto:a.kalogirou@leeds.ac.uk>> wrote:
Dear all,
I am trying to solve a simple problem (linear water wave equations), where I have the Laplace equation in a square domain with Neumann boundary conditions on the three sides and on the top boundary two equations are satisfied.
I have the following code: -------------------------------------- V = FunctionSpace(mesh,"CG",1)
eta0 = Function(V).interpolate(Expression("cos(2*pi*x[0])")) phi0 = Function(V) eta1 = Function(V) phi1 = Function(V) eta = TrialFunction(V) phi = TrialFunction(V) v = TestFunction(V)
aeta = v*eta*ds(4) Leta = v*eta0*ds(4) + dt*inner(grad(v),grad(phi0))*dx
eta_problem = LinearVariationalProblem(aeta,Leta,eta1) eta_solver = LinearVariationalSolver(eta_problem)
aphi = v*phi*ds(4) Lphi = (v*phi0 - g*dt*v*eta1)*ds(4)
phi_problem = LinearVariationalProblem(aphi,Lphi,phi1) phi_solver = LinearVariationalSolver(phi_problem) --------------------------------------
The error I get tells me that Object is in wrong state, and Matrix is missing diagonal entry 39.
I thought this had to do with the problem being singular, so I tried to build the nullspace but I don't know how to define the VectorSpaceBasis. Any help would be appreciated.
Thank you very much.
Regards, Anna.
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds
On 11 Dec 2014, at 19:17, Anna Kalogirou <a.kalogirou@leeds.ac.uk> wrote:
Thanks Andrew. About your comment 1), how can I fix this?
This required a change in PyOP2 (https://github.com/OP2/PyOP2/pull/417) to always create matrices with (at least) zeros in all diagonal entries. Further comments in line below.
The two equations valid on one of the boundaries can actually be combined into one by paying the price of a higher time derivative. This is not a problem for this simple system, so if I use
aphi = v*phi*ds(4) + g*dt**2*inner(grad(v),grad(phi))*dx Lphi = v*(2*phi1-phi0)*ds(4)
phi_problem = LinearVariationalProblem(aphi,Lphi,phi2) phi_solver = LinearVariationalSolver(phi_problem)
is equivalent to what I had before, and this time the code works and gives a good result (compared to the exact solution).
However, I need to understand why the other case doesn't work because that will later be generalised to become nonlinear.
Hopefully with that change, the problem will go away again.
On 11/12/14 13:04, Andrew McRae wrote:
Four comments:
1) If I understand correctly, PETSc insists that the diagonal entries of the sparse matrix are filled in. Usually, the LHS form will cause this to happen; however, that obviously isn't the case for your form!
2) It should be possible to build the nullspace you want, but I'm clueless on the actual syntax. This wouldn't get around (1) without further fixes from our end.
3) Could the equations given not be included in the weak form for the rest of the problem?
4) If not, there might be a hacky "add-a-weighted-mass-matrix-to-the-LHS" fix you could use?....
Andrew
On 11 December 2014 at 12:46, Anna Kalogirou <a.kalogirou@leeds.ac.uk> wrote: Dear all,
I am trying to solve a simple problem (linear water wave equations), where I have the Laplace equation in a square domain with Neumann boundary conditions on the three sides and on the top boundary two equations are satisfied.
I have the following code: -------------------------------------- V = FunctionSpace(mesh,"CG",1)
eta0 = Function(V).interpolate(Expression("cos(2*pi*x[0])")) phi0 = Function(V) eta1 = Function(V) phi1 = Function(V) eta = TrialFunction(V) phi = TrialFunction(V) v = TestFunction(V)
aeta = v*eta*ds(4) Leta = v*eta0*ds(4) + dt*inner(grad(v),grad(phi0))*dx
eta_problem = LinearVariationalProblem(aeta,Leta,eta1) eta_solver = LinearVariationalSolver(eta_problem)
aphi = v*phi*ds(4) Lphi = (v*phi0 - g*dt*v*eta1)*ds(4)
phi_problem = LinearVariationalProblem(aphi,Lphi,phi1) phi_solver = LinearVariationalSolver(phi_problem) --------------------------------------
The error I get tells me that Object is in wrong state, and Matrix is missing diagonal entry 39.
I thought this had to do with the problem being singular, so I tried to build the nullspace but I don't know how to define the VectorSpaceBasis. Any help would be appreciated.
So this operator certainly has a nullspace, which you can build, but it's quite large. Since any Function orthogonal to the boundary dofs is in the space. Can you describe how the documentation for the vector space basis is insufficient? To build it, you must provide a list of orthonormal Functions spanning the null space of the operator: basis = [] for i in range(nbasis): f = Function(V) f.interpolate(...) basis.append(f) nullspace = VectorSpaceBasis(basis) solve(...., nullspace=nullspace) This is the left null space of the operator (used to orthogonalise the solution). You need to arrange that your right hand side is consistent (i.e. orthogonal to the right null space of the operator [if the operator is symmetric this is the same as the left null space]). Lawrence
Dear all, I am experiencing again the same error I reported in December, when trying to solve an equation only on the boundary of a domain (the code is the same as the one in my email below dated 11 December 12:46). The problem was fixed back then but after all these months and after various updates, I get it again. Basically PETSc is complaining about Object being in wrong state and Matrix missing a diagonal entry. The error arises when trying to solve this VP: a = v*phi*ds(4) L = (v*phi0 - g*dt*v*eta0)*ds(4) Note that I get two different errors on two different computers (on the other computer I get a Zero pivot error in LU factorization, on row 0), but the one I am reporting here is based on the computer with the latest petsc/pyop2/firedrake update. Thank you! Best. Anna. On 12/12/14 12:32, Lawrence Mitchell wrote:
On 11 Dec 2014, at 19:17, Anna Kalogirou <a.kalogirou@leeds.ac.uk> wrote:
Thanks Andrew. About your comment 1), how can I fix this? This required a change in PyOP2 (https://github.com/OP2/PyOP2/pull/417) to always create matrices with (at least) zeros in all diagonal entries.
Further comments in line below.
The two equations valid on one of the boundaries can actually be combined into one by paying the price of a higher time derivative. This is not a problem for this simple system, so if I use
aphi = v*phi*ds(4) + g*dt**2*inner(grad(v),grad(phi))*dx Lphi = v*(2*phi1-phi0)*ds(4)
phi_problem = LinearVariationalProblem(aphi,Lphi,phi2) phi_solver = LinearVariationalSolver(phi_problem)
is equivalent to what I had before, and this time the code works and gives a good result (compared to the exact solution).
However, I need to understand why the other case doesn't work because that will later be generalised to become nonlinear. Hopefully with that change, the problem will go away again.
On 11/12/14 13:04, Andrew McRae wrote:
Four comments:
1) If I understand correctly, PETSc insists that the diagonal entries of the sparse matrix are filled in. Usually, the LHS form will cause this to happen; however, that obviously isn't the case for your form!
2) It should be possible to build the nullspace you want, but I'm clueless on the actual syntax. This wouldn't get around (1) without further fixes from our end.
3) Could the equations given not be included in the weak form for the rest of the problem?
4) If not, there might be a hacky "add-a-weighted-mass-matrix-to-the-LHS" fix you could use?....
Andrew
On 11 December 2014 at 12:46, Anna Kalogirou <a.kalogirou@leeds.ac.uk> wrote: Dear all,
I am trying to solve a simple problem (linear water wave equations), where I have the Laplace equation in a square domain with Neumann boundary conditions on the three sides and on the top boundary two equations are satisfied.
I have the following code: -------------------------------------- V = FunctionSpace(mesh,"CG",1)
eta0 = Function(V).interpolate(Expression("cos(2*pi*x[0])")) phi0 = Function(V) eta1 = Function(V) phi1 = Function(V) eta = TrialFunction(V) phi = TrialFunction(V) v = TestFunction(V)
aeta = v*eta*ds(4) Leta = v*eta0*ds(4) + dt*inner(grad(v),grad(phi0))*dx
eta_problem = LinearVariationalProblem(aeta,Leta,eta1) eta_solver = LinearVariationalSolver(eta_problem)
aphi = v*phi*ds(4) Lphi = (v*phi0 - g*dt*v*eta1)*ds(4)
phi_problem = LinearVariationalProblem(aphi,Lphi,phi1) phi_solver = LinearVariationalSolver(phi_problem) --------------------------------------
The error I get tells me that Object is in wrong state, and Matrix is missing diagonal entry 39.
I thought this had to do with the problem being singular, so I tried to build the nullspace but I don't know how to define the VectorSpaceBasis. Any help would be appreciated. So this operator certainly has a nullspace, which you can build, but it's quite large. Since any Function orthogonal to the boundary dofs is in the space.
Can you describe how the documentation for the vector space basis is insufficient?
To build it, you must provide a list of orthonormal Functions spanning the null space of the operator:
basis = [] for i in range(nbasis): f = Function(V) f.interpolate(...) basis.append(f)
nullspace = VectorSpaceBasis(basis)
solve(...., nullspace=nullspace)
This is the left null space of the operator (used to orthogonalise the solution). You need to arrange that your right hand side is consistent (i.e. orthogonal to the right null space of the operator [if the operator is symmetric this is the same as the left null space]).
Lawrence
-- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Hi Anna, On 29/04/15 14:50, Anna Kalogirou wrote:
Dear all,
I am experiencing again the same error I reported in December, when trying to solve an equation only on the boundary of a domain (the code is the same as the one in my email below dated 11 December 12:46). The problem was fixed back then but after all these months and after various updates, I get it again.
Basically PETSc is complaining about Object being in wrong state and Matrix missing a diagonal entry. The error arises when trying to solve this VP:
a = v*phi*ds(4) L = (v*phi0 - g*dt*v*eta0)*ds(4)
Note that I get two different errors on two different computers (on the other computer I get a Zero pivot error in LU factorization, on row 0), but the one I am reporting here is based on the computer with the latest petsc/pyop2/firedrake update.
Thanks, the missing diagonal entry was indeed a regression on our part. I have fixed this in PyOP2 so if you just update PyOP2 again it should go away. Additionally I actually added tests that the diagonals get created this time, so hopefully the same issue will not pop up again the next time one of us touches this part of the code. Now to address the zero pivot error. This is due to a change in PETSc: previously (certainly back in December) the default ILU preconditioner would shift the diagonal to try and avoid zero pivots. This was changed such that now the default is no shifting. With your operator, given that a large number of the entries /are/ zero, you therefore hit this zero pivot error. To rectify this, pass the solver parameter option: solver_parameters={'pc_factor_shift_type': 'nonzero'} If you're running with MPI, you should also pass the pair: 'sub_pc_factor_shift_type': 'nonzero' Since the default PC in parallel is process-block ILU. Cheers, Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 iQEcBAEBAgAGBQJVQQ0eAAoJECOc1kQ8PEYvgEEH+wY26EMvhg2RHR+hwBTVBruz DpJOGdyJaBIIgMkDCP4/F0Q7aNyAmX5VBjD4EOlLh3zUJu0DIUT3AkGpYV9cYQ16 nRSed7eKsTb03HaiNwcxrNxJ2cGIiQeBWdJ8TadX3HYEGlBacTI7oMAzAjp+G0V2 U4rOlq+6ilTdiPJMvfADMLxU70tV6Ni4QgsQNdT0yWygLA/cefFuRlM1D/jE4EPq Wr8NTp8e4bu6Z1q42sFiK0kIGahTu42hdYTnUHJx7g2mbYh3qBHEaILnb8f5yXR4 OPRaLhCjiQAjYMsjoYexXpBgCX/zlVibYsCuJ+f1EbXuyVIQ/FJG64XhW/Dzrk4= =aHDU -----END PGP SIGNATURE-----
Thanks Lawrence. Unfortunately a lot of tests fail after the PyOP2 installation... On 29/04/15 17:56, Lawrence Mitchell wrote:
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
Hi Anna,
On 29/04/15 14:50, Anna Kalogirou wrote:
Dear all,
I am experiencing again the same error I reported in December, when trying to solve an equation only on the boundary of a domain (the code is the same as the one in my email below dated 11 December 12:46). The problem was fixed back then but after all these months and after various updates, I get it again.
Basically PETSc is complaining about Object being in wrong state and Matrix missing a diagonal entry. The error arises when trying to solve this VP:
a = v*phi*ds(4) L = (v*phi0 - g*dt*v*eta0)*ds(4)
Note that I get two different errors on two different computers (on the other computer I get a Zero pivot error in LU factorization, on row 0), but the one I am reporting here is based on the computer with the latest petsc/pyop2/firedrake update. Thanks, the missing diagonal entry was indeed a regression on our part. I have fixed this in PyOP2 so if you just update PyOP2 again it should go away. Additionally I actually added tests that the diagonals get created this time, so hopefully the same issue will not pop up again the next time one of us touches this part of the code.
Now to address the zero pivot error. This is due to a change in PETSc: previously (certainly back in December) the default ILU preconditioner would shift the diagonal to try and avoid zero pivots. This was changed such that now the default is no shifting.
With your operator, given that a large number of the entries /are/ zero, you therefore hit this zero pivot error.
To rectify this, pass the solver parameter option:
solver_parameters={'pc_factor_shift_type': 'nonzero'}
If you're running with MPI, you should also pass the pair:
'sub_pc_factor_shift_type': 'nonzero'
Since the default PC in parallel is process-block ILU.
Cheers,
Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1
iQEcBAEBAgAGBQJVQQ0eAAoJECOc1kQ8PEYvgEEH+wY26EMvhg2RHR+hwBTVBruz DpJOGdyJaBIIgMkDCP4/F0Q7aNyAmX5VBjD4EOlLh3zUJu0DIUT3AkGpYV9cYQ16 nRSed7eKsTb03HaiNwcxrNxJ2cGIiQeBWdJ8TadX3HYEGlBacTI7oMAzAjp+G0V2 U4rOlq+6ilTdiPJMvfADMLxU70tV6Ni4QgsQNdT0yWygLA/cefFuRlM1D/jE4EPq Wr8NTp8e4bu6Z1q42sFiK0kIGahTu42hdYTnUHJx7g2mbYh3qBHEaILnb8f5yXR4 OPRaLhCjiQAjYMsjoYexXpBgCX/zlVibYsCuJ+f1EbXuyVIQ/FJG64XhW/Dzrk4= =aHDU -----END PGP SIGNATURE-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr Anna Kalogirou Research Fellow School of Mathematics University of Leeds http://www1.maths.leeds.ac.uk/~matak/
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 29/04/15 18:10, Anna Kalogirou wrote:
Thanks Lawrence.
Unfortunately a lot of tests fail after the PyOP2 installation...
Are these the pyop2 tests, firedrake tests, or your own? Note that you'll need to have updated the full set of dependencies, maybe also best to run firedrake/scripts/firedrake-clean Cheers, Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 iQEcBAEBAgAGBQJVQRH6AAoJECOc1kQ8PEYvn04H+wdzl2EGaUoE8pdbaVAGIn+o 6et39CeiwnlFJihr9cw/JXuItZbRvMCX50pv5MGOkrYjVdYHqyroWnazPUC55ods 3wWKWHIJQEHpmHx9dLplxpyHIal/oRmTV6MOzQW+X0z4kgPeDnvRin3UhaAUTruU pBBHI/zVTvvcj+raENbDy9fCxo8Z7PqC0N2r9NoFfF45xSUCB5PCqabvtUO4tjbU qnEd5FicioprdN8CL6PgEVe+I//Ya7q2/vGnXDhZwajOujTj3Bj+jF+oGqz+bOk2 LYRxxeM5P/+78NUgs+sG1xq/pcu37o2iGlQpj5fXR/7EosLUAxBpQXbVmnzwk7s= =xIKC -----END PGP SIGNATURE-----
participants (3)
-
Andrew McRae
-
Anna Kalogirou
-
Lawrence Mitchell