Best way to compute h-ellipticity?

[bueler]

(This is a support question not a bug or other real issue.)

Suppose I define an operator but not a smoother (splitting):

from lfa_lab import *                                                              
fine = Grid(2, [1.0/32, 1.0/32])                                                   
L = gallery.poisson_2d(fine)

I want the h-ellipticity as defined by Trottenberg et al seciton 4.7, i.e. E_h(L) = min(|lshigh|) / max(|ls|) where ls is the symbol of L and lshigh is the symbol over the high-frequency modes.

I can imagine various awkward methods for computing E_h(L) but I wonder what you would suggest?

PS Just found your nice library!

[hrittich]

This is, indeed, a good question. When computing the symbol of a linear operator, LFA Lab actually computes a block diagonal matrix, representing a sampling of the (discrete time) Fourier transform of the operator, i.e., each diagonal block stores the interactions of the harmonic frequencies. In the simplest case of a (constant) stencil operator, this matrix is just diagonal. You can get the matrix by calling .symbol().matrix() on an operator. Then the .block(i) method gives you the i-th block. The issue now is to figure out, which frequency belongs to a specific entry of a block, while the blocks and entries are ordered to make the implementation easy, which is not the most intuitive ordering.

Long story short, I probably should add one or more conveniece methods, to make the computation of h-ellipticity more straight forward.

[hrittich]

I have added a function to compute the h-ellipticity of an operator (see here). Turns out that there is a (more or less) simple way to use a hp_filter to extract the desired information. For more information feel free to look at the code.

[bueler]

Works for me! Thanks.