Compute diagonal with Matrix-Free on adaptive meshes

[Daniel Jodlbauer]

Dear all!

To verify my MatrixFree implementation, I compared its application to the classical matrix-vector multiplication (call it matrix A).
This is done by computing the matrix M of the operator MF by applying it to all unit-vectors.

However, when I compute the diagonal in the same way as LaplaceOperator does it (copied it), I get values different from the assembled diagonal once I have hanging nodes.
This error only occurs in the compute_diagonal function, i.e surprisingly M == A, but A(i,i) != compute_diagonal(i) (if hanging nodes are present)

My test-case starts with a 2x2 grid and refines one cell; apart from hanging node constraints no other constraints are active.

Some more details:

Constraints: (these look suspicious as well, shouldn't there be 4 constrained dofs instead of 2?)
    5 4:  0.5
    5 8:  0.5
    7 6:  0.5
    7 8:  0.5

Error occurs at (ignoring 5 and 7, since they are constrained)

dof | M, A | diag
 4  | 1.5  |  2.0
 6  | 1.5  |  2.0
 8  | 2.8  |  4.0

I will provide a mwe tomorrow, but maybe someone else has an idea already.


Thanks

Daniel

[Wolfgang Bangerth]

Daniel,

> To verify my MatrixFree implementation, I compared its application to

> the classical matrix-vector multiplication (call it matrix *A*).
> This is done by computing the matrix *M* of the operator *MF* by

> applying it to all unit-vectors.
>
> However, when I compute the diagonal in the same way as LaplaceOperator
> does it (copied it), I get values different from the assembled diagonal
> once I have hanging nodes.

The rows and columns corresponding to hanging nodes are empty with the
exception of the diagonal entry -- for which we use a value that has the
correct order of magnitude, but that is otherwise unspecified. In other
words, what diagonal value you have there is unimportant as long as it
is nonzero because these degrees of freedom don't couple with all of the
other degrees of freedom, and as long as you overwrite the values of the
computed solution through ConstraintMatrix::distribute().

It is not surprising to me that you get different diagonal entries when
you use two different methods. The question is whether you get a
different *solution* vector after distributing to hanging nodes.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Daniel Jodlbauer]

But thats exactly my point, the error occurs in the dofs which constrain a hanging node, not the hanging nodes (dofs 5 and 7) itself. I agree that the constrained dofs 5 and 7 can have arbitrary values.

I will check whether the solution is going to be different in any case. I was just afraid that wrong values on the diagonal could cause problems in the MGSmoother afterwards.

[Martin Kronbichler]

Dear Daniel,

the problem is, as far as I can tell, the fact that once you assemble into a matrix and once into a vector. The paper

K. Kormann: A Time-Space Adaptive Method for the Schrödinger Equation, Commun. Comput. Phys. 2016, doi: 10.4208/cicp.101214.021015a

describes this effect in section 5.3. Can you check there if this is what you see?

Best,
Martin

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Daniel Jodlbauer]

The effect described in the paper looks indeed similar, thanks for the hint.