Non-homogeneity in boundary conditions for EVPs

[Bindesh Tripathi]

Dear all,

While solving an eigenvalue problem in Dedalus, is it mandatory for all the eigenfunctions to vanish at the boundaries? Can I not have a boundary condition like below (a non-homogeneous equation)?

problem.add_bc("left(psi)   = -Lz" )

problem.add_bc("right(psi) =  Lz")

When I write the BCs in the above format, Dedalus suggests that the RHS is not homogeneous and exceeds the tolerance of 1.0e-10. When I write them in the following fashion:

problem.add_bc("left(psi)  + Lz = 0" )

problem.add_bc("right(psi) - Lz = 0")

an error is encountered -- RecursionError: maximum recursion depth exceeded while calling a Python object. 

I reckon that this is a pretty common error. Is there no way around to deal with non-homogeneous boundary condition for EVPs in Dedalus?

Context: I am attempting to do a 2D eigenvalue problem and hence wrote a script that generates tens/hundreds of equations coupled in the Fourier direction and the Chebyshev direction, which I then pass to Dedalus. It appears that somehow I would need to employ non-homogeneous boundary conditions (otherwise if I change my background state to make it satisfy homogeneous boundary conditions, the eigenvalue equation/problem no longer remains homogeneous). 

I would much appreciate if any suggestions masters have here dealing with non-homogeneity in BCs for EVPs.

Thanks a lot for your time and thoughts,

Bindesh

UW-Madison

[Geoffrey Vasil]

Bindesh, 

The answer to your question is most yes. “Eigenfunctions" need to vanish on the boundary. But I put “Eigenfunctions” in quotes, because it may not be exactly your eigenfunction, but maybe its derivative, or something like that. 

Here’s how to see it.

Putting aside differential equations for the time being, all eigenvalue problems reduce to effectively matrix problems. The most general form the can take is 

L(v) = s M(v),

where L, M are both *linear* operators and s is an eigenvalue. What is that exactly? It needs to satisfy the linearity property, 

L(a v + b w) = a L(v) + b L(w).

For *all* coefficient a, b. 

Same for M. The problem is that (regardless what you learned in high school) the function 

L(v)  + z 

is *nonlinear*. What to I mean? 

L(a v + b w) + z = a (L(v) + z) + b (L(w) + z) + (1-a-b) z.

This only works for the *particular coefficients* a+b = 1. But it needs to work *for all* a, b, so the constant terms means the function is not linear. You could even say it’s nonlinear. Even though that sounds crazy. A constant function is nonlinear. Its sub-linear. 

So what about boundary condition? This is linear:

(aψ + bφ)(x=0) = a ψ(x=0) + b φ(x=0)

You could also use the derivative. 

(aψ + bφ)'(x=0) = a ψ'(x=0) + b φ'(x=0)

This works because function evaluation and differentiation are both linear (in the strict sense). 

How does this fit into a generalised eigenvalue problem?

The 

L(v) = s M(v)

Part of the vector v could represent ψ. And part of the L could represent ψ(x=0) or left(ψ). And for that part of the L matrix, the corresponding M matrix could have the operation “multiply by zero”, which is also a linear operation. 

The M matrix could also not multiply by zero. In that case you could have something like 

ψ(x=0) = s ψ'(x=0) 

And the eigenvalue can be part of the boundary condition. This actually happens frequently in applications. 

But what you can’t have is 

ψ(x=0) - s ψ'(x=0) = 1 

This is non-linear. 

Here’s my favourite rule of thumb for nearly all math: 

Work out what happens in the discrete finite cases, and the continuous infinite thing will basically be the same. 

So imagine 

3 x + y + 2 z = s ( 2x + y + z)

2 x - 2y + 2z = s (  x + 0 y - 7 z) 

x + y + z       = 9 

So far so good with the firs two equations. But the last one really throws things off. If you get rid of the pesky 9, then 

3 x + y + 2 z = s ( 2x + y + z)

2 x - 2y + 2z = s (  x + 0 y - 7 z) 

x + y + z       = s ( 0 x + 0 y + 0 z) 

And this *is* a generalised eigenvalue problem. 

L(x,y,z) - s M(x,y,z)

The characteristic polynomials is det ( L - s M) = - 4 -11 s + 7 s* s

Notice that this only has degree 2. That’s ok. It means it only has two roots; s ~ -0.3, 1.8.  The last (0,0,0) row means the other eigenvalue is off at infinity. But it’s really not there at all. 

But hopefully all this makes it clear that *some* linear combination of the eigenfunction and its derivative does need to vanish at the boundary. 

Hope this helps!

Cheers, 

Geoff 

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/7f9166e8-fba2-4cd9-96d8-fb923d2c1bb6n%40googlegroups.com.

[BINDESH TRIPATHI]

Dear Geoff,

 

Thank you very much for guiding me in such great details – much appreciated! You are so excellent in demonstrating matters with cleanest paradigmatic examples. I adore your ability to educate finely!!!

 

Your explanation makes perfect sense – all we have got to do for the EVP calculations is to express the eigenfunctions (and their derivatives) in linear combination (no matter whether at the boundaries or in the bulk of the domain) so that each linear combination is separately zero. The task afterward is a straightforward matrix eigenvalue computation. This is excellent! I appreciate you highlighting this point.

 

Following your suggestion, I am now able to write a different set of equations that allows me (I believe) to deal with non-vanishing boundary conditions of “some function”.

 

problem.add_equation(“dt(psi) – ...*dz(psi) = 0”, tau=0)

 

#tau=0 ‘coz I will introduce two new equations and apply BCs on them as I do not wish to have the third equation floating around with only two BCs in hand. Perhaps I could write this equation as problem.add_equation (“dt(psi) – ...*psi_z = 0”) as well to avoid tau=0 ‘coz it’s simply an algebraic equation now.

 

problem.add_equation(“psi_z  – dz(psi)   = 0”)

problem.add_equation(“psi_zz – dz(psi_z) = 0”)

 

problem.add_bc(“left(psi_zz)  = 0”)

problem.add_bc(“right(psi_zz) = 0”)

 

 

Please note that I was earlier attempting to use the BCs such that psi goes to -Lz and +Lz at the left and right boundaries. This new set of equations with BCs applied to the higher-order derivatives (instead of the original function itself) appears to allow psi to become ±Lz at the boundaries. Please correct me if I am wrong. My line of thinking is:

 

  left(psi_zz) = 0

⇒ lim_{z ➝ -Lz} psi_zz = 0

⇒ lim_{z ➝ -Lz} psi_z  = constant

⇒ lim_{z ➝ -Lz} psi    = constant*z

⇒ left(psi)    = -constant*Lz

 

Similarly, right(psi) = constant*Lz.

 

However, I am unsure how to inform Dedalus that the constant should be 1. Does Dedalus inherently assume this to be the case?

 

I would highly treasure any comments if Geoff or other excellent folks here have to offer. If these above steps of imposing BCs are correct, then I trust I have been able to leverage Dedalus to do 2D EVP as I write hundreds of equations (not manually though!) for the variables psi at each kx, which, of course, are coupled in Fourier directions in a 2D EVP problem (I retain dz() terms though). It’s basically a higher dimensional version of 1D EVP where the eigenfunction now is: X = [psi_00(z), psi_01(z), psi_02(z), …] and the EVP becomes M dt(X) = L X.

 

Really appreciative of the Dedalus community, I am,

Bindesh.

--
You received this message because you are subscribed to a topic in the Google Groups "Dedalus Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dedalus-users/PeQLSCUArRU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/6F009D97-B062-489C-A175-071C501F3522%40gmail.com.

[christophe...@gmail.com]

You can use Python's string interpolation to add constants to the BC which aren't important enough to be parameters in your problem. Below I apply projective boundary conditions for a homoclinic orbit by pre-computing 3x3 matrices Lm and Lp, and interpolating them into the BCs using the f"...{value}..." syntax:

problem.add_bc(f"{Lm[0,0]}*left(u)  + {Lm[0,1]}*left(v)  + {Lm[0,2]}*left(uz)  = 0")
problem.add_bc(f"{Lp[0,0]}*right(u) + {Lp[0,1]}*right(v) + {Lp[0,2]}*right(uz) = 0")
problem.add_bc(f"{Lp[1,0]}*right(u) + {Lp[1,1]}*right(v) + {Lp[1,2]}*right(uz) = 0")

Maybe that will be useful for your problem?

[BINDESH TRIPATHI]

Many thanks, Christopher, for chipping-in the string interpolation (f-strings) in Python.

 

However, my issue with the BCs is more subtle than this. As you can see below, I am applying the BC on psi_zz and not on psi (so I can’t use f-strings here).

  left(psi_zz) = 0
⇒ lim_{z ➝ -Lz} psi_zz = 0
⇒ lim_{z ➝ -Lz} psi_z  = constant
⇒ lim_{z ➝ -Lz} psi    = constant*z
⇒ left(psi)    = -constant*Lz

I was hoping that by writing the first equation above in Dedalus, I am implying the last one holds true (which is what I want). However, the last one has a constant. I have no direct handle over the value of constant here. I want it to be 1. Please note that I can’t write left(psi) = -constant*Lz directly in Dedalus in EVP solver as it has no problem variable that multiplies the constant term (thus called non-homogeneous BC).

 

I was able to run the EVP calculation with left(psi_zz) = 0. However, I have not confident what it is producing is correct (because of the ambiguity of constant). Any further help or suggestions would be highly appreciated!

 

Thanks a lot,

Bindesh

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/6a88a730-0e86-4073-9fd3-8bcc17f23b28n%40googlegroups.com.

[BINDESH TRIPATHI]

Addendum: There should have been two constants as like below

 

left(psi_zz) = 0 
⇒ lim_{z ➝ -Lz} psi_zz = 0

⇒ lim_{z ➝ -Lz} psi_z  = constant1 
⇒ lim_{z ➝ -Lz} psi    = constant1*z + constant2
⇒ left(psi)    = -constant1*Lz + constant2

 

Thank you,

Bindesh

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CH2PR06MB6631322CFD77FB9D978608C5C72F9%40CH2PR06MB6631.namprd06.prod.outlook.com.

[BINDESH TRIPATHI]

Dear all,

 

After deeper contemplation over the weekend, I realized that there is a better way to tackle this problem of non-homogeneity in BCs or non-homogeneous term in eigenvalue equations (there is no reason why nature can’t have a constant term in the eigenvalue equation; as Geoff mentioned earlier, the problem variable of such an equation is NOT the “eigenfunction” and hence needs to be massaged to make it completely linear).

 

I thought it might be a good idea to record the resolution to this issue here (if somebody else happens to come across a similar problem, if not myself :) ).

 

I recast my state variable -- magnetic flux, psi, as:

psi = background (x-averaged in 2D simulations) + fluctuations (fluctuations go to zero at the boundaries; good for EVP) 

 

This yields in 2D magnetic induction equation, dt(psi) = {phi, psi} + 1/Rm \nabla^2{psi}, where phi(x,z) is the stream function:

 

dt(psi_k) – sum_{n=-kmax}^{kmax} {A_{nk} * psi_k} - sum_{n=-kmax}^{kmax} {B_{nk} * dz(psi_k)} = 0 ...(1)

 

where n runs from -kmax to +kmax via zero, but psi_k in dt(psi_k) is defined only for the fluctuating wavenumbers as the background is not supposed to evolve in time here. So, k in dt(psi_k) in the LHS runs from -kmax, -kmax+1,  …, -2, -1, 1, 2, kmax-1, kmax (because k=0 is the background). This is how I encounter the non-homogeneous term (or constant term independent of the eigenfunction), i.e., psi_{k=0} term floating around inside the sum_n. In the above equation, psi_k(z) refers to Fourier transform of psi(x,z) at wavenumber kx and I have dropped the resistive term to compute non-dissipative 2D eigenmodes (i.e., modes that have spatial structure in both kx- and z-directions).

 

With this prelude, let me share how I overcome this challenge of non-homogeneous term in EVP. Equation (1) above is of the form:

dt(X) = L.X + nonhomogeneous-term(function of z-coordinate) … (2)

 

where

X = [psi_{-kmax}, psi_{-kmax+1}, …, psi_{-2}, psi{-1}, psi_{1}, psi{2}, …, psi_{kmax-1}, psi(kmax)].

Each of the variables, psi_{…}, here are functions of z. Note that psi_{0} does not evolve over time and hence is not considered. Apart from psi_{0}, all other variables go to zero at the boundaries as they belong to the fluctuating spectrum.

 

The non-homogeneous term in equation (2) arises because Fourier-transforming the induction equation couples different Fourier modes, one of which is the background itself.

 

Equation (2) can be rewritten as, with N(z) being the non-homogeneous term (a similar, but more paradigmatic form of the problem can be found here):

dt(X(z,t)) = L. (X(z,t)) + N(z) … (3)

 

Now, consider a decomposition X(z,t) = S(z,t) + H(z), where S(z,t) and H(z) would be chosen later. (H(z) is like time-averaged component of the non-zero kx spectrum of psi.)

dt(S(z,t)) = L. (S(z,t)) + L. H(z) + N(z) …(4)

 

This equation now can be passed to an EVP solver in Dedalus if

L. H(z) + N(z) = 0 …(5)

 

which lends us:

dt(S(z,t)) = L. (S(z,t)) …(6)

 

Equation (6) is perfectly linear, and the boundary conditions on S(z,t) are zeros as it comes from the fluctuating spectrum (non-zero kx of psi). Equation (5) can be solved for H(z) as a boundary value problem because the linear operator L has d/dz term in it (but we know what exactly N(z) is from the x-averaged background state, so not a problem!).

 

Thus, this non-homogeneous EVP is neatly overcome by defining a time-averaged component of the psi for non-zero kx, which is represented by H(z) that corresponds to the term (I believe) H_0(z) exp(lambda*t) where lambda=0. (The full expansion of psi would be psi = H_0 exp(0*t) + sum_{r} H_r exp(r*t)). That’s it!

 

Thanks a lot for this wonderful community where lots of brainstorming and mutual sharing of ideas prosper,

Bindesh.