Re: [firedrake] operations on matrices
On 5 January 2015 at 09:48, Cotter, Colin J <colin.cotter@imperial.ac.uk> wrote:
Does it help if the matrices are columnwise-blocks? I could assemble the entire matrix in a "column kernel" and then perform surgery on it?
The issue with doing things via a column kernel is that you would lose all of the UFL infrastructure, so you would have to code your own assembly in C.
--cjc
On 05/01/15 09:39, David Ham wrote:
On 5 January 2015 at 09:37, Cotter, Colin J <colin.cotter@imperial.ac.uk <mailto:colin.cotter@imperial.ac.uk>> wrote:
Oh yes, good isolation of the problem.
If alpha also depends on values of A, do we have a problem there too?
Yes. Same problem. There is no current way for a parallel loop to read the entries of a matrix.
-cjc
On 05/01/15 09:33, David Ham wrote: > Hi Colin, > > There is no way for a parallel loop to read from a matrix. However
the
> operation you describe appears to be: > > assemble A > assemble alpha > > scale entries of A by the corresponding entries of alpha. > > The last step is clearly the problem. I wonder if this could be
achieved
> by some PETSc operation on the matrices. > > On 5 January 2015 at 09:02, Cotter, Colin J <
colin.cotter@imperial.ac.uk <mailto:colin.cotter@imperial.ac.uk>
> <mailto:colin.cotter@imperial.ac.uk <mailto:colin.cotter@imperial.ac.uk>>> wrote: > > Dear all, > Happy New Year! > > Perhaps I made the mistake of making some complex explanation before > asking my question. > > What is the best way to make adjustments to matrix entries as part of a > loop over elements? > > cheers > --cjc > > On 22/12/14 11:13, Cotter, Colin J wrote: > > Dear Firedrakers, > > I've been recently revisiting the "algebraic flux correction" > schemes > > of Dmitri Kuzmin, with the aim of getting a conservative+bounded > > advection scheme for temperature in our NWP setup. These schemes > involve > > the following steps: > > > > 1) Forming the consistent mass matrix (which is column-diagonal) > M_C for > > the temperature space. > > 2) Constructing the following matrix with the same sparsity as M_C: > > > > A_{ij} = (M_C)_{ij}(T_i-T_j) > > > > where T_i is the value of temperature at node i. > > > > 3) "Limiting" the matrix by replacing > > > > A_{ij} -> A_{ij}\alpha_{ij} > > > > where \alpha_{ij} depends on various field values at nodes i and > j (only > > needs to be evaluated when nodes i and j share an element). > > > > 4) Evaluating Ax where x is the vector containing 1s, and adding x to > > the RHS of mass-matrix projection equation before solving. > > > > My question is: how to implement this in an efficient and > parallel-safe > > way in the Firedrake/PyOP2 framework? In particular, step
(3)
> involves > > looping over elements, and correcting matrix entries. Also, I'm > not sure > > of the best way to assemble A. > > > > all the best > > --Colin > > > _______________________________________________ > firedrake mailing list > firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> <mailto:firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk>> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > > -- > Dr David Ham > Departments of Mathematics and Computing > Imperial College London > > http://www.imperial.ac.uk/people/david.ham > > > _______________________________________________ > firedrake mailing list > firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr David Ham Departments of Mathematics and Computing Imperial College London
http://www.imperial.ac.uk/people/david.ham
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr David Ham Departments of Mathematics and Computing Imperial College London http://www.imperial.ac.uk/people/david.ham
Ignoring the issue of acting on global matrices for a minute, what would be the way to do the following: 1) assemble a local matrix from a UFL expression 2) modify local matrix entries a_ij making use of field values f_i and f_j, 3) insert the modified local matrix into a global matrix This involves code generation for 1) but hand-tooling for (2-3). all the best --cjc On 05/01/15 10:02, David Ham wrote:
On 5 January 2015 at 09:48, Cotter, Colin J <colin.cotter@imperial.ac.uk <mailto:colin.cotter@imperial.ac.uk>> wrote:
Does it help if the matrices are columnwise-blocks? I could assemble the entire matrix in a "column kernel" and then perform surgery on it?
The issue with doing things via a column kernel is that you would lose all of the UFL infrastructure, so you would have to code your own assembly in C.
--cjc
On 05/01/15 09:39, David Ham wrote: > > > On 5 January 2015 at 09:37, Cotter, Colin J <colin.cotter@imperial.ac.uk <mailto:colin.cotter@imperial.ac.uk> > <mailto:colin.cotter@imperial.ac.uk <mailto:colin.cotter@imperial.ac.uk>>> wrote: > > Oh yes, good isolation of the problem. > > If alpha also depends on values of A, do we have a problem there too? > > > Yes. Same problem. There is no current way for a parallel loop to read > the entries of a matrix. > > -cjc > > On 05/01/15 09:33, David Ham wrote: > > Hi Colin, > > > > There is no way for a parallel loop to read from a matrix. However the > > operation you describe appears to be: > > > > assemble A > > assemble alpha > > > > scale entries of A by the corresponding entries of alpha. > > > > The last step is clearly the problem. I wonder if this could be achieved > > by some PETSc operation on the matrices. > > > > On 5 January 2015 at 09:02, Cotter, Colin J <colin.cotter@imperial.ac.uk <mailto:colin.cotter@imperial.ac.uk> <mailto:colin.cotter@imperial.ac.uk <mailto:colin.cotter@imperial.ac.uk>> > > <mailto:colin.cotter@imperial.ac.uk <mailto:colin.cotter@imperial.ac.uk> > <mailto:colin.cotter@imperial.ac.uk <mailto:colin.cotter@imperial.ac.uk>>>> wrote: > > > > Dear all, > > Happy New Year! > > > > Perhaps I made the mistake of making some complex explanation > before > > asking my question. > > > > What is the best way to make adjustments to matrix entries as > part of a > > loop over elements? > > > > cheers > > --cjc > > > > On 22/12/14 11:13, Cotter, Colin J wrote: > > > Dear Firedrakers, > > > I've been recently revisiting the "algebraic flux > correction" > > schemes > > > of Dmitri Kuzmin, with the aim of getting a > conservative+bounded > > > advection scheme for temperature in our NWP setup. These > schemes > > involve > > > the following steps: > > > > > > 1) Forming the consistent mass matrix (which is > column-diagonal) > > M_C for > > > the temperature space. > > > 2) Constructing the following matrix with the same > sparsity as M_C: > > > > > > A_{ij} = (M_C)_{ij}(T_i-T_j) > > > > > > where T_i is the value of temperature at node i. > > > > > > 3) "Limiting" the matrix by replacing > > > > > > A_{ij} -> A_{ij}\alpha_{ij} > > > > > > where \alpha_{ij} depends on various field values at nodes > i and > > j (only > > > needs to be evaluated when nodes i and j share an element). > > > > > > 4) Evaluating Ax where x is the vector containing 1s, and > adding x to > > > the RHS of mass-matrix projection equation before solving. > > > > > > My question is: how to implement this in an efficient and > > parallel-safe > > > way in the Firedrake/PyOP2 framework? In particular, step (3) > > involves > > > looping over elements, and correcting matrix entries. > Also, I'm > > not sure > > > of the best way to assemble A. > > > > > > all the best > > > --Colin > > > > > > _______________________________________________ > > firedrake mailing list > > firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> <mailto:firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk>> > <mailto:firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> <mailto:firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk>>> > >https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > > > > > > > -- > > Dr David Ham > > Departments of Mathematics and Computing > > Imperial College London > > > >http://www.imperial.ac.uk/people/david.ham > > > > > > _______________________________________________ > > firedrake mailing list > > firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> <mailto:firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk>> > > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > > _______________________________________________ > firedrake mailing list > firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> <mailto:firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk>> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake > > > > > -- > Dr David Ham > Departments of Mathematics and Computing > Imperial College London > > http://www.imperial.ac.uk/people/david.ham > > > _______________________________________________ > firedrake mailing list > firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> > https://mailman.ic.ac.uk/mailman/listinfo/firedrake >
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr David Ham Departments of Mathematics and Computing Imperial College London
http://www.imperial.ac.uk/people/david.ham
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 07/01/15 08:46, Colin Cotter wrote:
Ignoring the issue of acting on global matrices for a minute, what would be the way to do the following:
1) assemble a local matrix from a UFL expression 2) modify local matrix entries a_ij making use of field values f_i and f_j, 3) insert the modified local matrix into a global matrix
This involves code generation for 1) but hand-tooling for (2-3).
So (1) produces a kernel that assembles a single element tensor. (3) is just normal element tensor insertion. So the difficulty is inserting (2) between (1) and (3) to modify the element tensor as a result of f_i and f_j. I'm going to guess here that you're unable to express the action of (2) in terms of UFL expressions. So one hacky way to do this is as follows (but, it's hacky, there are probably better ways): 1. Generate the kernel for (1) 2. Generate the kernel that takes the output of (1), an nxm element tensor and reads field values to modify the element tensor. 3. Construct a "super-kernel" that takes all the argument from (1) and the additional field arguments from (2) which calls first 1., then 2. to generate the correctly modified element tensor. 4. Use the kernel from 3. as your assembly kernel in a par_loop (rather than the kernel from 1.) Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.10 (GNU/Linux) Comment: Using GnuPG with undefined - http://www.enigmail.net/ iQEcBAEBAgAGBQJUrP4dAAoJECOc1kQ8PEYviXgH/0fPeNVA9Z8QMEmqBUlUjUH/ 2fGlyICoYmxmMPJ2d0x2BzSCtiGvo3nMfjejQC7+lPyKlSjdMnkv/K5rLKCD4UZQ nsHU3Vtwv10c9doa/hjXsIohqPu5vmUvmjRr4O+A9/+HknfReXdg7VByTDiQIpTg glKDQRAGWdzGfroom8b5Jd4mwXr7GkHNHKFTU33PhBi4tYrJHDSGh0AY4ZcbfzvG Ql6o1+RAYKXR3njPhQFvwZLYyZ9OFy1xEOmauSb/DxPNDHf8UZ7gLdIVSa/HTwd/ z3WqXRrapCCLXhgDezFXyUXsEQcJhPHPHZ/3+P1HQ6fBkiu3p7qtcTySLBpBG60= =icRz -----END PGP SIGNATURE-----
On 07/01/15 09:36, Lawrence Mitchell wrote:
So one hacky way to do this is as follows (but, it's hacky, there are probably better ways):
1. Generate the kernel for (1)
2. Generate the kernel that takes the output of (1), an nxm element tensor and reads field values to modify the element tensor.
3. Construct a "super-kernel" that takes all the argument from (1) and the additional field arguments from (2) which calls first 1., then 2. to generate the correctly modified element tensor.
4. Use the kernel from 3. as your assembly kernel in a par_loop (rather than the kernel from 1.)
Thanks. So, how do I do that in practice, given that I want to write the kernel for 2 in a c string? cheers --cjc
participants (3)
-
Colin Cotter
-
David Ham
-
Lawrence Mitchell