Tutorial 02 Beyn's Global Eigenvalue Solver

In Tutorial 01 you learnt the basics of the Helmholtz solver. This tutorial will make you more familiar with the global eigenvalue solver. The solver implements an algorithm invented by W.-J. Beyn in [1]. The implementation closely follows [2].

Model set-up.

The model is the same Rijke tube configuration as in Tutorial 01:

using WavesAndEigenvalues.Helmholtz
mesh=Mesh("Rijke_mm.msh",scale=0.001) #load mesh
dscrp=Dict() #initialize model descriptor
dscrp["Interior"]=(:interior, ()) #define resonant cavity
dscrp["Outlet"]=(:admittance, (:Y,1E15)) #specify outlet BC
γ=1.4 #ratio of specific heats
ρ=1.225 #density at the reference location upstream to the flame in kg/m^3
Tu=300.0    #K unburnt gas temperature
Tb=1200.0    #K burnt gas temperature
P0=101325.0 # ambient pressure in Pa
A=pi*0.025^2 # cross sectional area of the tube
Q02U0=P0*(Tb/Tu-1)*A*γ/(γ-1) #the ratio of mean heat release to mean velocity Q02U0
x_ref=[0.0; 0.0; -0.00101] #reference point
n_ref=[0.0; 0.0; 1.00] #directional unit vector of reference velocity
n=0.01 #interaction index
τ=0.001 #time delay
dscrp["Flame"]=(:flame,(γ,ρ,Q02U0,x_ref,n_ref,:n,:τ,n,τ)) #flame dynamics
R=287.05 # J/(kg*K) specific gas constant (air)
speedofsound(x,y,z) = z<0. ? sqrt(γ*R*Tu) : sqrt(γ*R*Tb)
c=generate_field(mesh,speedofsound)
L=discretize(mesh,dscrp,c)

1006×1006-dimensional operator family: 

ω^2*M+K+n*exp(-iωτ)*Q+ω*Y*C

Parameters
----------
n	0.01 + 0.0im
λ	Inf + 0.0im
ω	0.0 + 0.0im
τ	0.001 + 0.0im
Y	1.0e15 + 0.0im

The global eigenvalue solver has a range of parameters. Mandatory, are only the linear operator family L you like to solve and the contour Γ.

The contour is specified as a list of complex numbers. Each of which defining a vertex of the polygon you want to scan for eigenvalues. In Tutorial 01 we specified the contour as rectangle by

Γ=[150.0+5.0im, 150.0-5.0im, 1000.0-5.0im, 1000.0+5.0im].*2*pi

4-element Array{Complex{Float64},1}:
 942.4777960769379 + 31.41592653589793im
 942.4777960769379 - 31.41592653589793im
 6283.185307179586 - 31.41592653589793im
 6283.185307179586 + 31.41592653589793im

The vertices are traversed in the order you specify them. However, the direction does not matter, meaning that equivalently you can specify Γ as:

Γ=[1000.0+5.0im, 1000.0-5.0im, 150.0-5.0im, 150.0+5.0im]

4-element Array{Complex{Float64},1}:
 1000.0 + 5.0im
 1000.0 - 5.0im
  150.0 - 5.0im
  150.0 + 5.0im

Note that you should not specify the first point as the last point. This will be automatically done for you by the solver.

Again, Γ is a polygon, so, e.g., specifying only three points will make the search region a triangle. And to get curved regions you just approximate them by a lot of vertices. For instance, this polygon is a good approximation of a circle.

radius=10
center_point=100
Γ_circle=radius.*[cos(α)+sin(α)*1.0im for α=0:2*pi/1000:(2*pi-2*pi/1000)].+center_point

1000-element Array{Complex{Float64},1}:
              110.0 + 0.0im
 109.99980260856137 + 0.06283143965558952im
 109.99921044203816 + 0.1256603988335261im
  109.9982235238081 + 0.18848439715408175im
   109.996841892833 + 0.2513009544333748im
 109.99506560365731 + 0.3141075907812829im
  109.9928947264059 + 0.3769018266993454im
 109.99032934678125 + 0.43968118317864907im
 109.98736956606018 + 0.5024431817976955im
 109.98401550108974 + 0.5651853448202453im
                    ⋮
 109.98401550108974 - 0.5651853448202481im
 109.98736956606018 - 0.5024431817976934im
 109.99032934678125 - 0.43968118317865074im
  109.9928947264059 - 0.37690182669934214im
 109.99506560365731 - 0.3141075907812836im
   109.996841892833 - 0.25130095443337047im
  109.9982235238081 - 0.18848439715408138im
 109.99921044203816 - 0.1256603988335207im
 109.99980260856137 - 0.06283143965558806im

Note that your contour does not need to be convex!

Number of eigenvalues

Beyn's algorithm needs an initial guess on how many eigenvalues l you expect inside your contour. Per default this is l=5. If this guess is less than the actual eigenvalues. The algorithms will miserably fail to compute any eigenvalue inside your contour correctly. You may, therefore, choose a large number for l to be on the save site. However, this will dramatically increase your computation time.

There are several strategies to deal with this problem and each of which might be well suited in a certain situation:

Split your domain in several smaller domains. This reduces the potential number of eigenvalues in each of the subdomains.
Just repeatedly rerun the solver with increasing values of l until the found eigenvalues do not change anymore.

Quadrature points

As Beyn's algorithm is based on contour integration a numerical quadrature method is performed under the hood. More precisely, the code performs a Gauss-Legendre integration on each of the edges of the Polygon Γ. The number of quadrature points for each of these integrations is N=32 per default. But you can specify it as an optional parameter when calling the solver. For instance as the circular contour is already discretized by 1000 points, it is not necessary to utilize N=16 quadrature points on each of the edges. Let's say N=4 should be fine. Then, you would call the solver like

Ω,P=beyn(L,Γ_circle,N=4,output=true)

(Complex{Float64}[], Array{Complex{Float64}}(undef,1006,0))

Test singular values

One step in Beyn's algorithm requires a singular value decomposition. Non-zero singular values and singular vectors associated with these are meant to be disregarded. Unfortunately, on a computer zero could also mean a very small number, and small here is problem-dependent. The code will disregard any singular value σ<tol where tol == 0.0 Per default. This means, that per default only singular values that are exactly 0 will be disregarded. It is very unlikely that this will be the case for any singular value. You may specify your own threshold by the optional key-word argument tol

Ω,P=beyn(L,Γ, tol=1E-10)

(Complex{Float64}[], Array{Complex{Float64}}(undef,1006,0))

but this is rather an option for experts. Even the authors of the code do not know a systematic way of well defining this threshold for given configuration.

Test the position of the eigenvalues

Because there could be a lot of spurious eigenvalues as there is no correct implementation of the singular value test. A very simple test to check whether your eigenvalues are correct is to see whether they are inside your contour Γ. This test is performed per default, but you can disable it using the keyword pos_test.

Ω,P=beyn(L,Γ, pos_test=false)

(Complex{Float64}[-227.59378160988854 + 0.013602819904629755im, 0.00036401480533119037 - 0.03890434385869311im, 177.83430676140006 + 5.429335912057693im, 558.885267571265 - 4.204510046668872im, 587.347194875147 + 6.160822425376081im], Complex{Float64}[0.035793838308427556 + 0.004487759414966951im 0.00822005351563242 + 0.00043895537317368357im … -0.026925492420329343 + 0.0010053296076950897im -0.0251758505732799 - 0.0007621835774363298im; 0.03269795943233084 + 0.003965455720839242im 0.008220095392531867 + 0.00043943607433099847im … -0.024286849171794474 + 0.00029399621664918344im -0.02738402091647882 - 0.0004657210439259702im; … ; 0.011972594819804066 + 0.0028930305605216504im 0.008220492659332803 + 0.0004345453228687188im … -0.034342167400046905 + 0.0022984425081145304im -0.023860223763101315 - 0.001284108555237522im; -0.04093846230311263 + 0.00037037002762738434im 0.008221385223486997 + 0.00042553514347187007im … -0.041045478106863396 - 0.003536716142963243im -0.06336953109667608 + 0.0036724651222049464im])

Note, that if you disable the position test and do not specify a threshold for the singular values, you will always be returned with l eigenvalues by Beyn's algorithm. (And sometimes even eigenvalues that are outside of your contour are fairly correct... )

Toggle the progress bar

Especially when searching for many eigenvalues, in large domains the algorithm may take some time to finish. You can, therefore, optionally display a progress bar together with an estimate of how long it will take for the routine to finish. This is done using the optional keyword output.

Ω,P=beyn(L,Γ, output=true)

(Complex{Float64}[558.885267571265 - 4.204510046668872im], Complex{Float64}[-0.026925492420329343 + 0.0010053296076950897im; -0.024286849171794474 + 0.00029399621664918344im; … ; -0.034342167400046905 + 0.0022984425081145304im; -0.041045478106863396 - 0.003536716142963243im])

Summary

Beyn's algorithm is a powerful tool for finding all eigenvalues in a specified domain. However, its computational costs are quite high and it therefore is best suited for solving small to medium-sized problems. A good strategy is to combine it with a local iterative solver such that solutions from Beyn's algorithm are fed into an iterative solver for further improvement. The next tutorial will further explain iterative solvers.

References

[1] W.-J. Beyn, An integral method for solving nonlinear eigenvalue problems, 2012, Lin. Alg. Appl., doi:10.1016/j.laa.2011.03.030

[2] P.E. Buschmann, G.A. Mensah, J.P. Moeck, Solution of Thermoacoustic Eigenvalue Problems with a Non-Iterative Method, J. Eng. Gas Turbines Power, Mar 2020, 142(3): 031022 (11 pages) doi:10.1115/1.4045076

This page was generated using Literate.jl.