Re: [firedrake] operations on matrices
On 5 January 2015 at 09:37, Cotter, Colin J <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>> 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> 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
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? --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
On 5 Jan 2015, at 09:48, Colin Cotter <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?
You can do whatever you like to the "element tensor" block before insertion. However, there's no way of getting the element tensor /out/ of a CSR matrix to do surgery on it in a parallel loop. Lawrence
On 5 Jan 2015, at 09:39, David Ham <David.Ham@imperial.ac.uk> wrote:
On 5 January 2015 at 09:37, Cotter, Colin J <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.
You could do this by creating a new output matrix, looping over the rows of the input matrix: out = A.duplicate() start, end = A.getOwnershipRange() for row in range(start, end): acol, aval = A.getRow(row) alphacol, alphaval = alpha.getRow(row) assert all(acol == alphacol), "Column pattern must match" val = aval * alphaval out.setValues(row, acol, val, addv=PETSc.InsertMode.INSERT_VALUES) out.assemble() Completely untested. Lawrence
participants (3)
-
Colin Cotter
-
David Ham
-
Lawrence Mitchell