Dear firedrakers, I’ve been trying to disentangle the setup costs in my various solvers again. Recall that I solve a velocity/pressure system and then precondition it with a Schur-complement preconditioner, which is either AMG or geometric multigrid for the pressure system. However, to solve the velocity/pressure system I need to assemble the u/p operator at some point and this cost will be the same in both solvers. In my matrix-free code I assemble matrices for the velocity/pressure system and then apply them as follows in my geometric multigrid: mat = assemble(EXPRESSION) # line 1 mat_handle = M.handle # line 2 with r.dat.vec as v: with u.dat.vec_ro as x: mat_handle.mult(x,v) I find that the really expensive bit is line 2. Is this the correct way of doing it? The setup seems to be much more expensive than the setup of the corresponding LinearVariationalSolver which I construct for the PETSc fieldsplit solver. In the PETSc solver (with AMG preconditioner) on the other hand I also assemble an operator for the mixed velocity/pressure system in the corresponding setup routine (in this routine I assemble the velocity/pressure operator and form a LinearVariationalSolver from it, and tell it to use the fieldsplit preconditoner). If I measure the time for this, it is much less than the setup time in the corresponding geometric multigrid operator. However, I then notice that the first AMG iteration is much more expensive than the next ones, so my suspicion is that this contains the setup costs for both the velocity/pressure system and for the AMG. To make any sense of my results, I really need to disentangle this somehow, and figure out what is AMG setup time and what is the cost of setting up the mixed operator. I might simply do something wrong in the setup of my mixed system. Thanks a lot, Eike
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Hi Eike, On 21/04/15 10:36, Eike Mueller wrote:
Dear firedrakers,
I’ve been trying to disentangle the setup costs in my various solvers again. Recall that I solve a velocity/pressure system and then precondition it with a Schur-complement preconditioner, which is either AMG or geometric multigrid for the pressure system. However, to solve the velocity/pressure system I need to assemble the u/p operator at some point and this cost will be the same in both solvers.
In my matrix-free code I assemble matrices for the velocity/pressure system and then apply them as follows in my geometric multigrid:
mat = assemble(EXPRESSION) # line 1 mat_handle = M.handle # line 2
with r.dat.vec as v: with u.dat.vec_ro as x: mat_handle.mult(x,v)
I find that the really expensive bit is line 2. Is this the correct way of doing it? The setup seems to be much more expensive than the setup of the corresponding LinearVariationalSolver which I construct for the PETSc fieldsplit solver.
In the PETSc solver (with AMG preconditioner) on the other hand I also assemble an operator for the mixed velocity/pressure system in the corresponding setup routine (in this routine I assemble the velocity/pressure operator and form a LinearVariationalSolver from it, and tell it to use the fieldsplit preconditoner). If I measure the time for this, it is much less than the setup time in the corresponding geometric multigrid operator. However, I then notice that the first AMG iteration is much more expensive than the next ones, so my suspicion is that this contains the setup costs for both the velocity/pressure system and for the AMG.
To make any sense of my results, I really need to disentangle this somehow, and figure out what is AMG setup time and what is the cost of setting up the mixed operator. I might simply do something wrong in the setup of my mixed system.
It is currently the case, although it will soon change as part of my attempts to reduce the cost of cheap iterations over the mesh, that Matrices are not actually created until you ask for the matrix handle. So there are two different levels of things going on: mat = assemble(expression) mat_handle = mat.M.handle This gives you a firedrake matrix that is a promise to produce an assembled operator. When you ask for the handle, two things happen: 1) The matrix is allocated and zeroed 2) The assembly is performed If you want to measure the time to allocate and zero the operator: mat_handle = mat._M.handle This gives you a pointer to the mat handle but doesn't perform the assembly. mat.M This forces the assembly. A (hopefully soon to land) merge will change the initialisation of PyOP2 matrices such that they are always allocated and zeroed immediately (rather than delaying until you ask for the matrix handle). Comparing to the LinearVariationalSolver: The assembly of the operator doesn't actually occur until the first call to solve, since it is at this point that the solver asks for the jacobian which forces the assembly. IOW, when you build the LVS, you again are only really doing symbolic level operations. To measure the cost of operator setup you can profile the form_jacobian method inside firedrake/solving_utils.py Hope that helps! Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 iQEcBAEBAgAGBQJVNhxpAAoJECOc1kQ8PEYvF2gH/3K66sbCaIyeYUNQBiWVURtw 7rAqk8mTRNbJuxryWQGT+T58mO21o7MP3/X0DHrFkKckPeK4SD1Kf1ul2yBm1AfQ lVBLmqWy2qDiaievL3TCGC5CqCKmkAFWs13iPnTtwIgKuH4rFYEFOlWLfLhCA7S7 5getCdTuJ1hzxDgRlhvbu4SfP5vq0hOWj/kUDvJsdqkDNR++ANhzQIJtyZiGbPg5 51wnETNZLXizuBRnFwEQz4U2aTuxTUQOcuQY4wcoDPSGr0RxhqbULTl70Bm9SIQL FtGSU2hLG4DU3qbKqxRWl21Nt01U2Ayr5Zpjp2HviGqp9eqGAgx7lKIxrm7yj98= =vqym -----END PGP SIGNATURE-----
Hi Lawrence, thanks! When I build a LinearSolver, can I solve the same problem (but for different RHSs) to two different tolerances? Currently I construct the LinearSolver twice, and that is very wasteful since I assemble the matrix twice). Thanks, Eike -- Dr Eike Hermann Mueller Lecturer in Scientific Computing Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom +44 1225 38 5557 e.mueller@bath.ac.uk<mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ On 21 Apr 2015, at 10:46, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk<mailto:lawrence.mitchell@imperial.ac.uk>> wrote: -----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Hi Eike, On 21/04/15 10:36, Eike Mueller wrote: Dear firedrakers, I’ve been trying to disentangle the setup costs in my various solvers again. Recall that I solve a velocity/pressure system and then precondition it with a Schur-complement preconditioner, which is either AMG or geometric multigrid for the pressure system. However, to solve the velocity/pressure system I need to assemble the u/p operator at some point and this cost will be the same in both solvers. In my matrix-free code I assemble matrices for the velocity/pressure system and then apply them as follows in my geometric multigrid: mat = assemble(EXPRESSION) # line 1 mat_handle = M.handle # line 2 with r.dat.vec as v: with u.dat.vec_ro as x: mat_handle.mult(x,v) I find that the really expensive bit is line 2. Is this the correct way of doing it? The setup seems to be much more expensive than the setup of the corresponding LinearVariationalSolver which I construct for the PETSc fieldsplit solver. In the PETSc solver (with AMG preconditioner) on the other hand I also assemble an operator for the mixed velocity/pressure system in the corresponding setup routine (in this routine I assemble the velocity/pressure operator and form a LinearVariationalSolver from it, and tell it to use the fieldsplit preconditoner). If I measure the time for this, it is much less than the setup time in the corresponding geometric multigrid operator. However, I then notice that the first AMG iteration is much more expensive than the next ones, so my suspicion is that this contains the setup costs for both the velocity/pressure system and for the AMG. To make any sense of my results, I really need to disentangle this somehow, and figure out what is AMG setup time and what is the cost of setting up the mixed operator. I might simply do something wrong in the setup of my mixed system. It is currently the case, although it will soon change as part of my attempts to reduce the cost of cheap iterations over the mesh, that Matrices are not actually created until you ask for the matrix handle. So there are two different levels of things going on: mat = assemble(expression) mat_handle = mat.M.handle This gives you a firedrake matrix that is a promise to produce an assembled operator. When you ask for the handle, two things happen: 1) The matrix is allocated and zeroed 2) The assembly is performed If you want to measure the time to allocate and zero the operator: mat_handle = mat._M.handle This gives you a pointer to the mat handle but doesn't perform the assembly. mat.M This forces the assembly. A (hopefully soon to land) merge will change the initialisation of PyOP2 matrices such that they are always allocated and zeroed immediately (rather than delaying until you ask for the matrix handle). Comparing to the LinearVariationalSolver: The assembly of the operator doesn't actually occur until the first call to solve, since it is at this point that the solver asks for the jacobian which forces the assembly. IOW, when you build the LVS, you again are only really doing symbolic level operations. To measure the cost of operator setup you can profile the form_jacobian method inside firedrake/solving_utils.py Hope that helps! Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 iQEcBAEBAgAGBQJVNhxpAAoJECOc1kQ8PEYvF2gH/3K66sbCaIyeYUNQBiWVURtw 7rAqk8mTRNbJuxryWQGT+T58mO21o7MP3/X0DHrFkKckPeK4SD1Kf1ul2yBm1AfQ lVBLmqWy2qDiaievL3TCGC5CqCKmkAFWs13iPnTtwIgKuH4rFYEFOlWLfLhCA7S7 5getCdTuJ1hzxDgRlhvbu4SfP5vq0hOWj/kUDvJsdqkDNR++ANhzQIJtyZiGbPg5 51wnETNZLXizuBRnFwEQz4U2aTuxTUQOcuQY4wcoDPSGr0RxhqbULTl70Bm9SIQL FtGSU2hLG4DU3qbKqxRWl21Nt01U2Ayr5Zpjp2HviGqp9eqGAgx7lKIxrm7yj98= =vqym -----END PGP SIGNATURE----- _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk<mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 21/04/15 11:20, Eike Mueller wrote:
Hi Lawrence,
thanks! When I build a LinearSolver, can I solve the same problem (but for different RHSs) to two different tolerances? Currently I construct the LinearSolver twice, and that is very wasteful since I assemble the matrix twice).
Yes: foo = LinearSolver(..., solver_parameters={'ksp_rtol': 1e-10}) foo.solve(...) foo.parameters = {'ksp_rtol': 1e-4} Note that setting this dict completely overwrites all the solver options. Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 iQEcBAEBAgAGBQJVNiULAAoJECOc1kQ8PEYvD0cH/jQTIJp1xRouxpWv02CydaLQ 7uKPkQ7RA2E2HfdSqEwJDz1cv7uVfi9mevlz8sY8sKfzb3QTodyzFY6/oFI+IFb6 W7Ae/d2LmhFRLv0aeGVNh7KneNHepZNVXMvzjxKmBNy7OQpnuifMVCVlE+KLYw/u M/+um5X5RRintrbq1NGNx/Lle6FSfds9r7cO2rgB/nNpI0FUHvCwuvOx21R1kNHe ym+DcJdCP3O/V/Nc5IzyHrhCbf5W40MGyNSW0wgw0gKAninYiGQwkyd49AnKemuf pijIZB3xjO1AcIt8X0fS/C5glotdpISc2QxqDGouc38MgOib4ov4gR6KGi2h77k= =qd87 -----END PGP SIGNATURE-----
Hi Lawrence, it looks like I’m setting up the same operator multiple times in my code, so this is very wasteful and explains why my setup costs are so high. Basically I need four assembled operators which describe the 2x2 matrix of the velocity-pressure system in several parts of my code. What happens if I construct those centrally, but then use them in different parts of the code? For example, I might want to apply the UU,UP,PU and PP-operators when I apply the mixed system matrix in the outer loop of the iterative solver (and I would do this in my MixedOperator class). For this I would extract the handles with .M.handle as below. But somewhere else I need to solve for the velocity mass matrix (I would do this in the class Mutilde), so I would throw the same UU operator into a LinearSolver class. Would the operator be built multiple times in this case? Or can I call mat.M once, and this will then be the only time where the operator gets assembled? If that is the case, then I need some class which holds the UU,UP,PU and PP operators and then feed them to the other classes appropriately. Thanks, Eike
On 21 Apr 2015, at 11:23, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
On 21/04/15 11:20, Eike Mueller wrote:
Hi Lawrence,
thanks! When I build a LinearSolver, can I solve the same problem (but for different RHSs) to two different tolerances? Currently I construct the LinearSolver twice, and that is very wasteful since I assemble the matrix twice).
Yes:
foo = LinearSolver(..., solver_parameters={'ksp_rtol': 1e-10})
foo.solve(...)
foo.parameters = {'ksp_rtol': 1e-4}
Note that setting this dict completely overwrites all the solver options.
Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1
iQEcBAEBAgAGBQJVNiULAAoJECOc1kQ8PEYvD0cH/jQTIJp1xRouxpWv02CydaLQ 7uKPkQ7RA2E2HfdSqEwJDz1cv7uVfi9mevlz8sY8sKfzb3QTodyzFY6/oFI+IFb6 W7Ae/d2LmhFRLv0aeGVNh7KneNHepZNVXMvzjxKmBNy7OQpnuifMVCVlE+KLYw/u M/+um5X5RRintrbq1NGNx/Lle6FSfds9r7cO2rgB/nNpI0FUHvCwuvOx21R1kNHe ym+DcJdCP3O/V/Nc5IzyHrhCbf5W40MGyNSW0wgw0gKAninYiGQwkyd49AnKemuf pijIZB3xjO1AcIt8X0fS/C5glotdpISc2QxqDGouc38MgOib4ov4gR6KGi2h77k= =qd87 -----END PGP SIGNATURE-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 21/04/15 11:42, Eike Mueller wrote:
Hi Lawrence,
it looks like I’m setting up the same operator multiple times in my code, so this is very wasteful and explains why my setup costs are so high.
Basically I need four assembled operators which describe the 2x2 matrix of the velocity-pressure system in several parts of my code. What happens if I construct those centrally, but then use them in different parts of the code? For example, I might want to apply the UU,UP,PU and PP-operators when I apply the mixed system matrix in the outer loop of the iterative solver (and I would do this in my MixedOperator class). For this I would extract the handles with .M.handle as below. But somewhere else I need to solve for the velocity mass matrix (I would do this in the class Mutilde), so I would throw the same UU operator into a LinearSolver class. Would the operator be built multiple times in this case? Or can I call mat.M once, and this will then be the only time where the operator gets assembled? If that is the case, then I need some class which holds the UU,UP,PU and PP operators and then feed them to the other classes appropriately.
So you build UU, UP, PU, and PP separately? Given they never change, I would assemble them up front and then just use them as and where. If you always do (whereever you need it): UU = assemble(...) Then you will pay the cost multiple times. If you build everything up front then calling mat.M a second time doesn't reassemble (since now the promise to assemble is fulfilled) and the operator is ready. Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 iQEcBAEBAgAGBQJVNitHAAoJECOc1kQ8PEYvrB8H/0SW4RYh7IEIsv2UjU9O6vyW DEVlLdiHMJTmbgJu9kJ7N9/6LabdPqBOVr15AB/IpRbm3GD0BahOevCaFbI3Ck+r 7hyhgZyBAsZXzggI+cjXKPOTJtvmXrvCzM4uoK0EHbK496gKwoa4o8xc6JrUKt+5 7J6594g6PH6nudIzd51Gm6YRpF0ldPU6HvJ5Hg9EUuSV1R5mtlPvhM8BU9cxpB4N Cq+BR0EnSvruL238nTvSHMyRitrqPe01eHFPby426xAzPFrcdy3WERLOO6GQXN1k IVuGl1cFV5L8ylsocxOlv/6odNEJX7+06dvPqLizu+haSoqbMtb8t3k2Uqq1LOw= =g9Oo -----END PGP SIGNATURE-----
Hi Lawrence,
So you build UU, UP, PU, and PP separately?
Yes, but as you say, I should only build them once, not multiple times, which is incredibly wasteful.
Given they never change, I would assemble them up front and then just use them as and where. If you always do (whereever you need it):
UU = assemble(...)
Then you will pay the cost multiple times.
If you build everything up front then calling mat.M a second time doesn't reassemble (since now the promise to assemble is fulfilled) and the operator is ready.
Ok, good, so that sounds like what I’m going to do will work and I only pay the cost of assembling the mixed operator once. Thanks, Eike
Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1
iQEcBAEBAgAGBQJVNitHAAoJECOc1kQ8PEYvrB8H/0SW4RYh7IEIsv2UjU9O6vyW DEVlLdiHMJTmbgJu9kJ7N9/6LabdPqBOVr15AB/IpRbm3GD0BahOevCaFbI3Ck+r 7hyhgZyBAsZXzggI+cjXKPOTJtvmXrvCzM4uoK0EHbK496gKwoa4o8xc6JrUKt+5 7J6594g6PH6nudIzd51Gm6YRpF0ldPU6HvJ5Hg9EUuSV1R5mtlPvhM8BU9cxpB4N Cq+BR0EnSvruL238nTvSHMyRitrqPe01eHFPby426xAzPFrcdy3WERLOO6GQXN1k IVuGl1cFV5L8ylsocxOlv/6odNEJX7+06dvPqLizu+haSoqbMtb8t3k2Uqq1LOw= =g9Oo -----END PGP SIGNATURE-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
Eike Mueller
-
Lawrence Mitchell