Eigenvalues for Orr-Sommerfeld Problems

[Chris Camobreco]

Hello,

Following the online documentation, I setup an eigenvalue problem for Poiseuille flow at Re=2000, alpha = 1, beta = 0, U(y) = 1-y**2. This returned the correct leading eigenvalue to that documented in Schmid and Henningson 2001 (to the 7 sig figs listed therein). I also checked that the eigenvalue from dedalus matched my own MATLAB code based on Trefethen 2000 at the same parameters (which matched to 11 sig figs). It thus appeared that I was setting up the eigenvalue problem correctly in dedalus.

However, the problem I am interested in is not pressure driven Poiseuille flow, but is a wall-driven quasi-two-dimensional Shercliff flow. The equations solved are almost identical to the NSE's, excepting the inclusion of a linear friction term

 -H/Re\vect{u}; H is the friction parameter. This modifies the base flow to U(y) = cosh(H**0.5*y)/cosh(H**0.5), with impermeable lateral duct walls moving at unit streamwise velocity. At H=10, Re=71210.9, alpha =0.979651, beta = 0, my MATLAB code consistently returns a growth rate of -5.940844843682833e-04 with anywhere from 200-800 Chebyshev points across the duct.

Sadly, for this setup, dedalus does not seem to be accurately predicting the growth rate. I have tested primitive pressure-velocity (1), velocity-vorticity (2), Orr-Sommerfeld (3) and adjoint Orr-Sommerfeld (4) formulations, and the eigenvalues returned are the following:

These all have 600 points across the domain

(1 - single) pressure-velocity:         -5.944846823044626534e-04  

(1 - compound) pressure-velocity: -5.950779968328887568e-04

(2 - single) velocity-vorticity:           -5.940739103038325053e-04

(2 - compound) velocity-vorticity:   -5.940018491754990495e-04

(3 - single) forward OS:                     -5.931086385353027176e-04

(3 - compound) forward OS:            -5.930365748691260623e-04

(4 - single) adjoint OS:                      -5.935194119948383447e-04

(4 - compound) adjoint OS:              -5.941127054618644283e-04

MATLAB (elsewhere validated):       -5.940844843682833e-04 

So something is up. Actually, a few things seem off. First, with 600 points, the results should easily be resolution independent, so switching from a single to a compound basis should not make a difference (I would have thought). Furthermore, targeting eigenvalues near zero, I still require a relatively deep search (100-150 eigenvalues) to find the correct leading mode. It also appears to be requiring far more resolution to become resolution independent in either the single or compound cases than my MATLAB equivalent.

Python scripts to compute eigenvalues for the various formulations follow:

(3) and (4) - there are two places where I comment/uncomment; one to adjust the y-basis and the other for whether I am solving the forward or adjoint OS equation

import numpy as np

import h5py

import time

from dedalus import public as de

from dedalus.extras import flow_tools

import logging

logger = logging.getLogger(__name__)

Reynolds = 71210.9

Har = 10

alpha = 0.979651

nys = 150

yb1 = de.Chebyshev('y1',nys, interval=(-1, -0.5))

yb2 = de.Chebyshev('y2',nys, interval=(-0.5, 0))

yb3 = de.Chebyshev('y3',nys, interval=(0, 0.5))

yb4 = de.Chebyshev('y4',nys, interval=(0.5, 1))

# UNCOMMENT FOR COMPOUND BASIS

# y_basis = de.Compound('y',(yb1,yb2,yb3,yb4))

ny = 600

# UNCOMMENT FOR SINGLE BASIS

y_basis = de.Chebyshev('y',ny, interval=(-1, 1))

domain = de.Domain([y_basis], grid_dtype=np.complex128)

problem = de.EVP(domain, variables=['v','vy','vyy','vyyy'],eigenvalue='omega')

problem.parameters['Re'] = Reynolds

problem.parameters['kx'] = alpha

problem.parameters['H'] = Har

problem.parameters['sqH'] = Har**0.5

problem.substitutions['UB'] = 'cosh(sqH*y)/cosh(sqH)'

problem.substitutions['UBy'] = 'sqH*sinh(sqH*y)/cosh(sqH)'

problem.substitutions['UByy'] = 'H*cosh(sqH*y)/cosh(sqH)'

# COMMENT/UNCOMMENT FIRST/SECOND EQUATIONS TO ADJUST FROM FORWARD TO ADJOINT SOLVE

problem.add_equation("-1j*omega*(vyy - kx**2*v) - 1/Re*(dy(vyyy) - 2*kx**2*vyy + kx**4*vy) - 1j*kx*UByy*v + UB*1j*kx*(vyy - kx**2*v) + (H/Re)*(vyy - kx**2*v) = 0")

# problem.add_equation("-1j*omega*(vyy - kx**2*v) - 1/Re*(dy(vyyy) - 2*kx**2*vyy + kx**4*vy) - 2*1j*kx*UBy*vy - UB*1j*kx*(vyy - kx**2*v) + (H/Re)*(vyy - kx**2*v) = 0")

problem.add_equation("vy - dy(v) = 0")

problem.add_equation("vyy - dy(vy) = 0")

problem.add_equation("vyyy - dy(vyy) = 0")

problem.add_bc("left(v) = 0")

problem.add_bc("right(v) = 0")

problem.add_bc("left(vy) = 0")

problem.add_bc("right(vy) = 0")

solver =  problem.build_solver()

# Solve for eigenvalues with sparse search near zero, rebuilding NCCs

solver.solve_sparse(solver.pencils[0], N=150, target=0,tol=0)

# solver.solve_dense(solver.pencils[0])

logger.info('Max growth rate: %e' %np.max(solver.eigenvalues.imag))

revals = solver.eigenvalues.real

ievals = solver.eigenvalues.imag

revecs = solver.eigenvectors.real

ievecs = solver.eigenvectors.imag

np.savetxt('revals.dat',revals)

np.savetxt('ievals.dat',ievals)

np.savetxt('revecs.dat',revecs)

np.savetxt('ievecs.dat',ievecs)

# SEPARATE TEST TO ENSURE THE BASE FLOW IS BEING READ IN CORRECTLY

# ENERGY MATCHES MATLAB

y = domain.grid(0)

f = domain.new_field(name='f')

f['g'] = (np.cosh((Har**0.5)*y)/np.cosh((Har**0.5)))**2

# f['g'] = (1-(y**2))**2

EB = y_basis.Integrate(f)

EBN = EB.evaluate()

print(EBN['g'][0])

(2) PARTIAL CODE - THE RELEVANT SECTION TO SWAP TO VELOCITY-VORTICITY

problem = de.EVP(domain, variables=['o','v','oy','vy'],eigenvalue='omega')

problem.parameters['Re'] = Reynolds

problem.parameters['kx'] = alpha

problem.parameters['H'] = Har

problem.parameters['sqH'] = Har**0.5

problem.add_equation("-1j*omega*o  - 1/Re*(-kx*kx*o + dy(oy)) + (H/Re)*o + (cosh(sqH*y)/cosh(sqH))*1j*kx*o - v*H*(cosh(sqH*y)/cosh(sqH)) = 0")

problem.add_equation("o - 1j*kx*v - dy(vy)/(1j*kx) = 0")

problem.add_equation("oy - dy(o) = 0")

problem.add_equation("vy - dy(v) = 0")

problem.add_bc("left(vy) = 0")

problem.add_bc("right(vy) = 0")

problem.add_bc("left(v) = 0")

problem.add_bc("right(v) = 0")

(1) PARTIAL CODE - THE RELEVANT SECTION TO SWAP TO PRESSURE-VELOCITY

problem = de.EVP(domain, variables=['v','p','vy','py'],eigenvalue='omega')

problem.parameters['Re'] = Reynolds

problem.parameters['kx'] = alpha

problem.parameters['H'] = Har

problem.parameters['sqH'] = Har**0.5

problem.substitutions['UB'] = 'cosh(sqH*y)/cosh(sqH)'

problem.substitutions['UBy'] = 'sqH*sinh(sqH*y)/cosh(sqH)'

problem.add_equation("-1j*omega*v + dy(p) - 1/Re*(-kx*kx*v + dy(vy)) + UB*1j*kx*v = -(H/Re)*v")

# THE PRESSURE EQUATION IS FROM THE DIVERGENCE OF THE NSE's

problem.add_equation("-kx*kx*p + dy(py) + 2*UBy*1j*kx*v = 0")

problem.add_equation("vy - dy(v) = 0")

problem.add_equation("py - dy(p) = 0")

problem.add_bc("left(v) = 0")

problem.add_bc("right(v) = 0")

problem.add_bc("left(vy) = 0")

problem.add_bc("right(vy) = 0")

Any help identifying why there is such variance between setups would be much appreciated. I ask not only for the sake of this linear computation, but also as I plan to use linear initial conditions (from either dedalus, or a MATLAB code, interpolated to dedalus's Chebyshev basis) for nonlinear simulations. However, initial nonlinear tests appear to be requiring more resolution than other tests with an in-house solver, which makes me wonder if I am doing something wrong (and this may be the simplest example replicating the issue).

Thanks,

Chris

[Geoffrey Vasil]

Hey Chris, 

I’ve just glanced at this quickly. So I probably need to have a closer look. But a couple questions did pop up.

I notice you have non-constant coefficient functions that involve cosh(y) ect? If so, my recollection of the classic Orr-Sommerfeld equation is that it only has linear background flow. 

It’s not that there’s anything wrong with using a different background flow, but if you are using a cosh-like profile, then that might be related to the issue. 

Dedalus expands non-constant coefficients in a Chebyshev series and converts the results to operators. If this isn’t done accurately enough, then it will throw the results off more than just about anything else. 

So, I perhaps my suggestion is play with the NCC_cutoff parameter. Other than that, I’ll have a close look at your script a bit later and see if anything else becomes clear. 

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/f8c3bf05-8b76-4c75-95bf-e9c463f21c02n%40googlegroups.com.

[Keaton Burns]

Hi Chris,

Besides the ncc_cutoff that Geoff mentioned, I think you need so many modes because the full eigenvalue of the mode you’re looking for has an O(1) real part.  If you set target=1 (maybe a reasonable prior given the cosh background flow) and ncc_cutoff=1e-12, it looks like Dedalus finds the right mode and matches the growth rate to 8 digits with a resolution of ny=128 and N=10 modes in the sparse search. This is for the pressure-velocity and velocity-vorticity formulations at least — I’m still seeing larger errors in the OS formulation though.

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/2F22CFB2-568E-4327-A8BB-8D2175AE1578%40gmail.com.

[Chris Camobreco]

Hi,

Thanks for the rapid responses. Adjusting the ncc_cutoff and the target parameters, I am seeing the results you described for the lower order setup (8SF agreement). 

With the OS formulation, I typed those out this morning just to cover all my bases before asking about possible errors, and evidently was a little hasty. I have typed "kx**4*vy", when it should have been "kx**4*v". Correcting this my forward problem now gives -5.940844891243494616e-04, and so is also fixed. Sorry.

Assuming ncc_cutoff affects all computations, e.g. creating initial fields and integrals of perturbation=nonlinear - base components when using analysis task handlers, then I expect my nonlinear cases should fall in line as well.

Thanks again,

Chris

[Geoffrey Vasil]

Glad it seems that things seems to be working. 

There’s one thing unrelated to the numerical issues that might help make things less painful in the future. 

You can make substitutions for your multiplication be horizontal wavenumber. 

Something like, 

problem.substitutions[‘dy(u)’] = “1j*ky*u”

You could then do something like 

problem.substitutions[‘dyy(u)’] = “dy(dy(u))”

The nice think is that it’s easier to avoid (or at least notice) minus sign errors. 

And depending on how you use the substitutions, you can sometimes use the same (or similar) equations to setup a time-stepping problem without having to get rid of the wavenumber directly. 

Anyhow, it’s not a big deal for the problem you are working with. But it can come in handy sometimes. 

Cheers, 

Geoff 

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/030ac22e-1ca6-4907-95f1-2532732f9510n%40googlegroups.com.