From 371320eb2b7b9550390dfaca0ad1e4086a4b8ec1 Mon Sep 17 00:00:00 2001 From: upadp Date: Wed, 26 Jun 2019 13:45:44 +0200 Subject: [PATCH 1/6] Minor addition to test --- docs/src/hydrotutorial.md | 6 ++++++ docs/src/hydrotutotial.md | 1 + 2 files changed, 7 insertions(+) create mode 100644 docs/src/hydrotutorial.md create mode 100644 docs/src/hydrotutotial.md diff --git a/docs/src/hydrotutorial.md b/docs/src/hydrotutorial.md new file mode 100644 index 000000000..cadc0ee76 --- /dev/null +++ b/docs/src/hydrotutorial.md @@ -0,0 +1,6 @@ +# Tutorial: Hydrodynamic stability as a fourth order PEP + +## Background +Stability analysis of flows is a very important problem in fluid mechanics. Of particular interest are applications where a transition from laminar to turbulent flow is observed. +NEP +\lamb \ No newline at end of file diff --git a/docs/src/hydrotutotial.md b/docs/src/hydrotutotial.md new file mode 100644 index 000000000..d75bb4a90 --- /dev/null +++ b/docs/src/hydrotutotial.md @@ -0,0 +1 @@ +# Tutorial: Hydrodynamic stability as a fourth order PEP \ No newline at end of file From b7da5bedb02b352dbd5573fd1500ff8aec5f9720 Mon Sep 17 00:00:00 2001 From: upadp Date: Thu, 27 Jun 2019 16:02:14 +0200 Subject: [PATCH 2/6] Minor --- docs/make.jl | 1 + docs/src/hydrotutorial.md | 11 ++++++++--- docs/src/hydrotutotial.md | 1 - docs/src/test.jl | 2 ++ 4 files changed, 11 insertions(+), 4 deletions(-) delete mode 100644 docs/src/hydrotutotial.md create mode 100644 docs/src/test.jl diff --git a/docs/make.jl b/docs/make.jl index 87c74ff7d..6626cf6e1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -20,6 +20,7 @@ makedocs( "Tutorial 4 (Python NEP)" => "tutorial_call_python.md", "Tutorial 5 (MATLAB 1)" => "tutorial_matlab1.md", "Tutorial 6 (FORTRAN 1)" => "tutorial_fortran1.md", + "Tutorial 7 (Hydrodynamic stability PEP)" => "hydrotutorial.md" ] ) diff --git a/docs/src/hydrotutorial.md b/docs/src/hydrotutorial.md index cadc0ee76..93f2f3746 100644 --- a/docs/src/hydrotutorial.md +++ b/docs/src/hydrotutorial.md @@ -1,6 +1,11 @@ # Tutorial: Hydrodynamic stability as a fourth order PEP ## Background -Stability analysis of flows is a very important problem in fluid mechanics. Of particular interest are applications where a transition from laminar to turbulent flow is observed. -NEP -\lamb \ No newline at end of file +Stability analysis of flows is a very important problem in fluid mechanics. Of particular interest are applications where a transition from laminar to turbulent flow is observed and this begins with an instability in the laminar flow. Considering a mean base laminar flow $$\overline{U} = \begin{pmatrix}U(y)& 0& 0\end{pmatrix}$$ and the perturbation $$u' = \begin{pmatrix}u& v& w\end{pmatrix}$$. Linearizing the Navier-Stokes equations around the mean flow gives us the Orr-Sommerfeld and Squire equations + +```math +\dfrac{\partial u}{\partial t}+U\dfrac{\partial u}{\partial x}+vU' = -\dfrac{\partial p}{\partial x}+\frac{1}{Re}\nabla^2u\\ +\dfrac{\partial v}{\partial t}+U\dfrac{\partial v}{\partial x} = -\dfrac{\partial p}{\partial y}+\frac{1}{Re}\nabla^2v\\ +\dfrac{\partial w}{\partial t}+U\dfrac{\partial w}{\partial x} = -\dfrac{\partial p}{\partial z}+\frac{1}{Re}\nabla^2w +``` +Taking \ No newline at end of file diff --git a/docs/src/hydrotutotial.md b/docs/src/hydrotutotial.md deleted file mode 100644 index d75bb4a90..000000000 --- a/docs/src/hydrotutotial.md +++ /dev/null @@ -1 +0,0 @@ -# Tutorial: Hydrodynamic stability as a fourth order PEP \ No newline at end of file diff --git a/docs/src/test.jl b/docs/src/test.jl new file mode 100644 index 000000000..71f1ad5be --- /dev/null +++ b/docs/src/test.jl @@ -0,0 +1,2 @@ +M(λ) +λ \ No newline at end of file From 2b1c1822726e915749c84e10b899e50ed156ac9b Mon Sep 17 00:00:00 2001 From: Parikshit Upadhyaya Date: Wed, 3 Jul 2019 00:01:14 +0200 Subject: [PATCH 3/6] Upto formualtion as a PEP --- docs/src/hydrotutorial.md | 40 +++++++++++++++++++++++++++++++++------ 1 file changed, 34 insertions(+), 6 deletions(-) diff --git a/docs/src/hydrotutorial.md b/docs/src/hydrotutorial.md index 93f2f3746..93b3e0d7c 100644 --- a/docs/src/hydrotutorial.md +++ b/docs/src/hydrotutorial.md @@ -1,11 +1,39 @@ -# Tutorial: Hydrodynamic stability as a fourth order PEP +# Tutorial: Stability of parallel shear flows ## Background -Stability analysis of flows is a very important problem in fluid mechanics. Of particular interest are applications where a transition from laminar to turbulent flow is observed and this begins with an instability in the laminar flow. Considering a mean base laminar flow $$\overline{U} = \begin{pmatrix}U(y)& 0& 0\end{pmatrix}$$ and the perturbation $$u' = \begin{pmatrix}u& v& w\end{pmatrix}$$. Linearizing the Navier-Stokes equations around the mean flow gives us the Orr-Sommerfeld and Squire equations +Stability analysis of flows is a very important problem in fluid mechanics. Of particular interest are applications where a transition from laminar to turbulent flow is observed and this begins with an instability in the laminar flow. Considering a mean base laminar flow $$\overline{U} = \begin{pmatrix}U(y)& 0& 0\end{pmatrix}$$ and the perturbation $$u' = \begin{pmatrix}u& v& w\end{pmatrix}$$. The normal vorticity is defined as $$\eta = \dfrac{\partial u}{\partial z}-\dfrac{\partial w}{\partial x}$$. Linearizing the Navier-Stokes equations around the mean flow and then eliminating pressure gives us the Orr-Sommerfeld and Squire equations, which are a system of fourth order PDEs describing the dynamics of $v$ and $\eta$. ```math -\dfrac{\partial u}{\partial t}+U\dfrac{\partial u}{\partial x}+vU' = -\dfrac{\partial p}{\partial x}+\frac{1}{Re}\nabla^2u\\ -\dfrac{\partial v}{\partial t}+U\dfrac{\partial v}{\partial x} = -\dfrac{\partial p}{\partial y}+\frac{1}{Re}\nabla^2v\\ -\dfrac{\partial w}{\partial t}+U\dfrac{\partial w}{\partial x} = -\dfrac{\partial p}{\partial z}+\frac{1}{Re}\nabla^2w +\left(\Big(\dfrac{\partial }{\partial t}+U\dfrac{\partial }{\partial x}\Big)\nabla^2-U''\dfrac{\partial }{\partial x}-\frac{1}{Re}\nabla^4\right)v = 0,\\ +\left(\dfrac{\partial }{\partial t}+U\dfrac{\partial }{\partial x}-\frac{1}{Re}\nabla^2\right)\eta = -U''\dfrac{\partial v}{\partial z}. ``` -Taking \ No newline at end of file + +## Formulation as a nonlinear eigenvalue problem +The ansatz $v(x,y,z,t) = \tilde{v}(y)\exp(i(\alpha x+\beta z-\omega t))$ and $\eta(x,y,z,t) = \tilde{\eta}(y)\exp(i(\alpha x+\beta z-\omega t))$, implying plane wave perturabtions transforms the system to + +```math +\left((-i\omega+i\alpha U)(\mathcal{D}^2-\alpha^2-\beta^2)-i\alpha U''-\frac{1}{Re}(\mathcal{D}^2-\alpha^2-\beta^2)^2\right)\tilde{v} = 0,\\ +\left((-i\omega+i\alpha U)-\frac{1}{Re}(\mathcal{D}^2-\alpha^2-\beta^2)\right)\eta = -i\beta U'\tilde{v}, +``` +where $\mathcal{D}$ denotes the 1-D differential operator $\frac{\partial }{\partial y}$. In this example, we consider the boundary conditions $\tilde{v} = \mathcal{D}\tilde{v} = \tilde{\eta} = 0$. + +We assume that $$\beta$$ and $$\omega$$ are given and we wish to solve for $$\alpha$$. This is usually done by using a transformation of the form $$\begin{pmatrix}\tilde{v}\\ \tilde{\eta}\end{pmatrix} = \begin{pmatrix}\tilde{V}\\ \tilde{E}\end{pmatrix}\exp(-\alpha y)$$ which reduces the power of $$\alpha$$ from four to two. The problem is then discretized and solved as a quadratic eigenvalue problem. Since NEP-PACK allows us to tackle polynomial eigenvalue problems of any order, we can simply discretize $$\mathcal{D}$$ to $$D$$ and expand, which leads to a polynomial eigenvalue problem of fourth order. + +We define diagonal matrices $$U_0$$, $$U_1$$ and $$U_2$$ as follows. +```math +U_0(i,i) = U(y_i),\\ +U_1(i,i) = U'(y_i),\\ +U_2(i,i) = U''(y_i).\\ +``` +Here $$\{y_i\}_{i=1,\ldots,n}$$ denotes the y-coordinates of the $$n$$ grid points used for discretization. The discretized problem on expansion is + +```math +\left[-\frac{I}{Re}\alpha^4-iU_0\alpha^3+\left(\frac{2D^2}{Re}+\left(i\omega- +\frac{2\beta^2}{Re}\right)I\right)\alpha^2+\left(iU_0(D^2-\beta^2I)-U_2\right)\alpha +\Big(\frac{2\beta^2D^2}{Re}-\frac{D^4}{Re}-\frac{\beta^4I}{Re}+i\omega(\beta^2I-D^2)\Big)\right]\tilde{v} = 0\\ +i\beta U_1\tilde{v}+\left[\frac{I}{Re}\alpha^2 + iU_0 \alpha +\left(\left(\frac{\beta^2}{Re}-i\omega\right)I-\frac{D^2}{Re}\right)\right]\tilde{\eta} +``` + + +## Implementation in NEP-PACK + + From 47256086edfa1167bccb8d36badfe02ff0d7bbad Mon Sep 17 00:00:00 2001 From: Parikshit Upadhyaya Date: Wed, 3 Jul 2019 23:00:22 +0200 Subject: [PATCH 4/6] Matrices for PEP added --- docs/src/hydrotutorial.md | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/docs/src/hydrotutorial.md b/docs/src/hydrotutorial.md index 93b3e0d7c..90ed64dc2 100644 --- a/docs/src/hydrotutorial.md +++ b/docs/src/hydrotutorial.md @@ -30,10 +30,18 @@ Here $$\{y_i\}_{i=1,\ldots,n}$$ denotes the y-coordinates of the $$n$$ grid poin ```math \left[-\frac{I}{Re}\alpha^4-iU_0\alpha^3+\left(\frac{2D^2}{Re}+\left(i\omega- \frac{2\beta^2}{Re}\right)I\right)\alpha^2+\left(iU_0(D^2-\beta^2I)-U_2\right)\alpha +\Big(\frac{2\beta^2D^2}{Re}-\frac{D^4}{Re}-\frac{\beta^4I}{Re}+i\omega(\beta^2I-D^2)\Big)\right]\tilde{v} = 0\\ -i\beta U_1\tilde{v}+\left[\frac{I}{Re}\alpha^2 + iU_0 \alpha +\left(\left(\frac{\beta^2}{Re}-i\omega\right)I-\frac{D^2}{Re}\right)\right]\tilde{\eta} +i\beta U_1\tilde{v}+\left[\frac{I}{Re}\alpha^2 + iU_0 \alpha +\left(\left(\frac{\beta^2}{Re}-i\omega\right)I-\frac{D^2}{Re}\right)\right]\tilde{\eta} = 0 ``` ## Implementation in NEP-PACK - +The discretized problem can be seen as a polynomial eigenvalue problem where $$M(\lambda) = A_0+A_1\lambda+A_2\lambda^2+A_3\lambda^3+A_4\lambda^4$$. The matrices are +```math +A_0 = \begin{pmatrix}\frac{2\beta^2D^2}{Re}-\frac{D^4}{Re}-\frac{\beta^4I}{Re}+i\omega(\beta^2I-D^2)& 0\\ i\beta U_1&\left(\frac{\beta^2}{Re}-i\omega\right)I-\frac{D^2}{Re}\end{pmatrix}\\ +A_1 = \begin{pmatrix}-iU_0& 0\\0& 0\end{pmatrix}\\ +A_2 = \begin{pmatrix}\frac{2D^2}{Re}+\left(i\omega- +\frac{2\beta^2}{Re}\right)I& 0\\0& \frac{I}{Re}\end{pmatrix}\\ +A_3 = \begin{pmatrix}iU_0(D^2-\beta^2I)-U_2& 0\\0& iU_0\end{pmatrix}\\ +A_4 = \begin{pmatrix}-\frac{I}{Re}& 0\\0& 0\end{pmatrix} +``` From 113eebf7490a3380ebdc53cf5d085d12ce23b6fa Mon Sep 17 00:00:00 2001 From: upadp Date: Fri, 9 Aug 2019 12:47:12 +0200 Subject: [PATCH 5/6] Hydrotutorial draft final version --- docs/src/hydro_cheb/cheb4c.jl | 72 ++++++++++++++++++++ docs/src/hydro_cheb/chebdif.jl | 55 +++++++++++++++ docs/src/hydrotutorial.md | 118 +++++++++++++++++++++++++++++---- docs/src/test.jl | 2 - 4 files changed, 233 insertions(+), 14 deletions(-) create mode 100644 docs/src/hydro_cheb/cheb4c.jl create mode 100644 docs/src/hydro_cheb/chebdif.jl delete mode 100644 docs/src/test.jl diff --git a/docs/src/hydro_cheb/cheb4c.jl b/docs/src/hydro_cheb/cheb4c.jl new file mode 100644 index 000000000..d211620c1 --- /dev/null +++ b/docs/src/hydro_cheb/cheb4c.jl @@ -0,0 +1,72 @@ +using LinearAlgebra + +function cheb4c(N) + + # The function x,D4 = cheb4c(N) computes the fourth derivative + # matrix D4 on Chebyshev interior points, incorporating the + # clamped boundary conditions. + # It is adapted from the function with the same name in dmsuite + # developed for Matlab by JAC Weideman and SC Reddy. + # + # Adapted by M. Beneitez at KTH Mechanics. + # beneitez@mech.kth.se + # + # Input: + # N-2: Order of the differentiation matrix. + # + # Output: + # yF: Interior Chebyshev points + # D4: Size [N-2,N-2] Fourth derivative matrix + # + + eye = Matrix{Float64}(I, N-2, N-2); + + n1 = Int(floor(N/2-1)); n2 = Int(ceil(N/2-1)); # Indices used for flipping trick. + + k = collect(1:N-2); # Compute theta vector. + th = k.*pi./(N-1) + + x = sin.(pi.*collect(N-3:-2:3-N)./(2*(N-1))); # Compute Chebyshev points. + + s = [sin.(th[1:n1]);reverse(sin.(th[1:n2]), dims = 1)] + + alpha = s.^4; + beta1 = -4*s.^2 .*x./alpha; + beta2 = 4*(3*x.^2 .-1)./alpha; + beta3 = 24*x./alpha; + beta4 = 24 ./alpha; + B = [beta1'; beta2'; beta3'; beta4']; + + T = repeat(th/2,1,N-2); + DX = 2*sin.(T'+T).*sin.(T'-T); # Trigonometric identity. + DX = [DX[1:n1,:]; -reverse(reverse(DX[1:n2,:], dims = 2), dims=1)]; # Flipping trick. + DX[eye.==true] = ones(N-2,1); # Put 1's on the main diagonal of DX. + + ss = s.^2 .*(-1.0).^k; + S = repeat(ss, 1, N-2); + C = S./S'; + + Z = 1 ./DX; + Z[eye.==true] = zeros(N-2,1); + + X = Z'; + X = X[X .!= 0.0]; + X = reshape(X,N-3,N-2); + + Y = ones(N-3,N-2); + D = Matrix{Float64}(I, N-2, N-2); + DM = zeros(N-2,N-2,4); + + for ell = 1:4 + Y = cumsum([B[ell,:]'; ell.*Y[1:N-3,:].*X], dims = 1); + D = ell.*Z.*(C.*repeat(diag(D), 1, N-2)-D); + D[eye.==true] = Y[N-2,:]; + DM[:,:,ell] = D; + end + + D4 = DM[:,:,4]; + + return x,D4 + + +end diff --git a/docs/src/hydro_cheb/chebdif.jl b/docs/src/hydro_cheb/chebdif.jl new file mode 100644 index 000000000..56156f5b4 --- /dev/null +++ b/docs/src/hydro_cheb/chebdif.jl @@ -0,0 +1,55 @@ +using LinearAlgebra +using ToeplitzMatrices + +function chebdif(N, M) + + # The function x,DM = chebdif(N,M) computes the differentiation + # matrices D1, D2, ..., DM on Chebyshev nodes. + # It is adapted from the function with the same name in dmsuite + # developed for Matlab by JAC Weideman and SC Reddy. + # + # Adapted by M. Beneitez at KTH Mechanics. + # beneitez@mech.kth.se + # + # Input: + # N: Size of the differentiation matrix. + # M: Number of derivates required + # Note: 0 < M <= N-1. + # + # Output: + # DM: DM[1:N,1:N,ell] contains the ell-th derivative matrix, ell=1...M + # + + eye = Matrix{Float64}(I, N, N); # Identity matrix. + + n1 = Int(floor(N/2)); n2 = Int(ceil(N/2)); # Indices used for flipping trick. + + k = collect(0:N-1); # Compute theta vector. + th = k.*pi./(N-1) + + x = sin.(pi.*collect(N-1:-2:1-N)./(2*(N-1))); # Compute Chebyshev points. + + T = repeat(th/2,1,N); + DX = 2*sin.(T'+T).*sin.(T'-T); # Trigonometric identity. + DX = [DX[1:n1,:]; -reverse(reverse(DX[1:n2,:], dims = 2), dims=1)]; # Flipping trick. + DX[eye.==true] = ones(N,1); # Put 1's on the main diagonal of DX. + + C = Matrix(Toeplitz((-1).^vec(k),(-1).^vec(k))); # C is the matrix with + C[1,:] = C[1,:].*2; C[N,:] = C[N,:].*2; # entries c[k]/c[j] + C[:,1] = C[:,1]/2; C[:,N] = C[:,N]/2 + + Z = 1 ./DX; # Z contains entries 1/(x[k]-x[j]) + Z[eye.==true] = zeros(N,1); # with zeros on the diagonal. + + D = eye; # D contains diff(). matrices. + DM = zeros(N,N,M); + + for ell = 1:M + D = ell*Z.*(C.*repeat(diag(D),1,N) - D); # Off-diagonals + D[eye.==true] = -sum(D', dims = 1); # Correct main diagonal of D + DM[:,:,ell] = D; # Store current D in DM + end + + return x,DM + +end diff --git a/docs/src/hydrotutorial.md b/docs/src/hydrotutorial.md index 90ed64dc2..69224ede7 100644 --- a/docs/src/hydrotutorial.md +++ b/docs/src/hydrotutorial.md @@ -15,27 +15,28 @@ The ansatz $v(x,y,z,t) = \tilde{v}(y)\exp(i(\alpha x+\beta z-\omega t))$ and $\e \left((-i\omega+i\alpha U)(\mathcal{D}^2-\alpha^2-\beta^2)-i\alpha U''-\frac{1}{Re}(\mathcal{D}^2-\alpha^2-\beta^2)^2\right)\tilde{v} = 0,\\ \left((-i\omega+i\alpha U)-\frac{1}{Re}(\mathcal{D}^2-\alpha^2-\beta^2)\right)\eta = -i\beta U'\tilde{v}, ``` -where $\mathcal{D}$ denotes the 1-D differential operator $\frac{\partial }{\partial y}$. In this example, we consider the boundary conditions $\tilde{v} = \mathcal{D}\tilde{v} = \tilde{\eta} = 0$. +where $\mathcal{D}$ denotes the 1-D differential operator $\frac{\partial }{\partial y}$. In this example, we consider the boundary conditions $\tilde{v} = \mathcal{D}\tilde{v} = \tilde{\eta} = 0$.We assume that $$\beta$$ and $$\omega$$ are given and we wish to solve for $$\alpha$$. -We assume that $$\beta$$ and $$\omega$$ are given and we wish to solve for $$\alpha$$. This is usually done by using a transformation of the form $$\begin{pmatrix}\tilde{v}\\ \tilde{\eta}\end{pmatrix} = \begin{pmatrix}\tilde{V}\\ \tilde{E}\end{pmatrix}\exp(-\alpha y)$$ which reduces the power of $$\alpha$$ from four to two. The problem is then discretized and solved as a quadratic eigenvalue problem. Since NEP-PACK allows us to tackle polynomial eigenvalue problems of any order, we can simply discretize $$\mathcal{D}$$ to $$D$$ and expand, which leads to a polynomial eigenvalue problem of fourth order. -We define diagonal matrices $$U_0$$, $$U_1$$ and $$U_2$$ as follows. +This is usually done by using a transformation of the form $$\begin{pmatrix}\tilde{v}\\ \tilde{\eta}\end{pmatrix} = \begin{pmatrix}\tilde{V}\\ \tilde{E}\end{pmatrix}\exp(-\alpha y)$$ which reduces the power of $$\alpha$$ from four to two. The problem is then discretized and solved as a quadratic eigenvalue problem. See [Chapter 7, Stability and Transition in Shear Flows, Schmid, Peter J., Henningson, Dan S] for details. Since NEP-PACK allows us to tackle polynomial eigenvalue problems of any order, we can simply discretize $$\mathcal{D}$$ to $$D$$ (using a suitable numerical discretization) and expand, which leads to a polynomial eigenvalue problem of fourth order. + +We define diagonal matrices $$U_0$$, $$U_1$$ and $$U_2$$ as follows. ```math -U_0(i,i) = U(y_i),\\ -U_1(i,i) = U'(y_i),\\ -U_2(i,i) = U''(y_i).\\ +U_0 = diag(U(y_0),U(y_1),\ldots,U(y_n)),\\ +U_1 = diag(U'(y_0),U'(y_1),\ldots,U'(y_n)),\\ +U_2 = diag(U''(y_0),U''(y_1),\ldots,U''(y_n)).\\ ``` -Here $$\{y_i\}_{i=1,\ldots,n}$$ denotes the y-coordinates of the $$n$$ grid points used for discretization. The discretized problem on expansion is - +Here $$\{y_i\}_{i=1,\ldots,n}$$ denotes the y-coordinates of the $$n$$ grid points used for discretization. The discretized problem on expansion is ```math \left[-\frac{I}{Re}\alpha^4-iU_0\alpha^3+\left(\frac{2D^2}{Re}+\left(i\omega- \frac{2\beta^2}{Re}\right)I\right)\alpha^2+\left(iU_0(D^2-\beta^2I)-U_2\right)\alpha +\Big(\frac{2\beta^2D^2}{Re}-\frac{D^4}{Re}-\frac{\beta^4I}{Re}+i\omega(\beta^2I-D^2)\Big)\right]\tilde{v} = 0\\ i\beta U_1\tilde{v}+\left[\frac{I}{Re}\alpha^2 + iU_0 \alpha +\left(\left(\frac{\beta^2}{Re}-i\omega\right)I-\frac{D^2}{Re}\right)\right]\tilde{\eta} = 0 ``` - - -## Implementation in NEP-PACK -The discretized problem can be seen as a polynomial eigenvalue problem where $$M(\lambda) = A_0+A_1\lambda+A_2\lambda^2+A_3\lambda^3+A_4\lambda^4$$. The matrices are +This can be be formulated as a polynomial eigenvalue problem where +```math +M(\lambda) = A_0+A_1\lambda+A_2\lambda^2+A_3\lambda^3+A_4\lambda^4, +``` +and the matrices $$A_i,_{i=1,\ldots,n}$$ are given by ```math A_0 = \begin{pmatrix}\frac{2\beta^2D^2}{Re}-\frac{D^4}{Re}-\frac{\beta^4I}{Re}+i\omega(\beta^2I-D^2)& 0\\ i\beta U_1&\left(\frac{\beta^2}{Re}-i\omega\right)I-\frac{D^2}{Re}\end{pmatrix}\\ A_1 = \begin{pmatrix}-iU_0& 0\\0& 0\end{pmatrix}\\ @@ -44,4 +45,97 @@ A_2 = \begin{pmatrix}\frac{2D^2}{Re}+\left(i\omega- A_3 = \begin{pmatrix}iU_0(D^2-\beta^2I)-U_2& 0\\0& iU_0\end{pmatrix}\\ A_4 = \begin{pmatrix}-\frac{I}{Re}& 0\\0& 0\end{pmatrix} ``` +## Problem setup +We begin by initializing the parameters to the values used to generate the data in Table 7.1 in Schmid-Henningson. +```julia +## Parameters +N = 256; # Number of interior points +Re = 2000; # Reynolds number +ω = 0.3; # Input frequency +β = 0.0; # Spanwise wavenumber +``` +To set up the matrices $$A_i$$, we need the discretized matrices corresponding to the operators $$\mathcal{D}^2$$ and $$\mathcal{D}^4$$. Here, we do this using Chebyshev nodes. +```julia +## Çhebyshev discretization of differential operators +yF,DM = chebdif(N+2, 4); +D2 = DM[2:N+1,2:N+1,2]; #D^2 +yF,D4 = cheb4c(N+2); +eye = Matrix{Float64}(I, N, N); #D^4 +``` +The code for implementation of `chebdif()` and `cheb4c()` is provided in [GITHUB LINK]. We can now set up the coefficient matrices and create a corresponding PEP object. +```julia +#Setuo coefficient matrices +A4 = [-eye/Re zeros(N,N);zeros(N,N) zeros(N,N)]; +A3 = [-1im*diagm(0 => U) zeros(N,N);zeros(N,N) zeros(N,N)]; +A2 = [(1im*ω-2β^2/Re)*eye+2D2/Re zeros(N,N);zeros(N,N) eye/Re]; +A1 = [1im*(diagm(0 => U)*(D2-eye*β^2)-Upp*eye) zeros(N,N);zeros(N,N) 1im*diagm(0=>U)]; +A0 = [2β^2*D2/Re-D4/Re-β^4*eye/Re+1im*ω*(β^2*eye-D2) zeros(N,N);1im*β*diagm(0 => Up) (-1im*ω+β^2/Re)*eye-D2/Re]; + +#Create a PEP object +nep = PEP([A0,A1,A2,A3,A4]); +``` +## Shifting and scaling +Before we proceed to using one of NEP-PACK's methods, we have to consider the issue of ill-conditioned coefficient matrices. We can use the reference eigenvalue $\lambda_0 = 0.97857+0.044394i$ to do a singularity test on $$M(\lambda_0)$$ to see this. +```julia +julia> norm(A0) +8.37194982854379e13 + +julia> norm(A1) +473949.06740743306 + +julia> norm(A2) +303285.4108872535 + +julia> norm(A3) +9.817076958035932 + +julia> norm(A4) +0.008 + +julia> minimum(svdvals(compute_Mder(nep,0.97857+1im*0.044394))) +0.00030839497238639225 +``` +We can get around this issue by scaling the PEP with NEP-PACK's `shift_and_scale()` , and solving the scaled problem $$T(\lambda) = M(100\lambda)$$ instead. +```julia +sc=100; +nep1 = shift_and_scale(nep,scale=sc); +``` +The scaled PEP performs much better in the singularity test. +```julia +julia> minimum(svdvals(compute_Mder(nep2,0.0097857+1im*0.00044394))) +3.8570057996438436e-10 +``` + +## Solution +In this example, we are interested in computing several eigenvalues and our region of interest for the spectrum is in the first quadrant. We use the Tensor Infinite Arnoldi (TIAR) method, as in [REFERENCE TO IAR]. The method is called twice with different shifts `σ`. +```julia +λ1,v1 = tiar(nep2,σ=0.006,v=ones(size(nep,1)),logger=1,neigs=10,maxit=200,tol=1e-14) +λ2,v2 = tiar(nep2,σ=0.005+0.005i,v=ones(size(nep,1)),logger=1,neigs=10,maxit=200,tol=1e-14) +``` +The computed eigenvalues are scaled back to get the eigenvalues of the original problem. +```julia + julia> λ_orig = 100*λtotal +20-element Array{Complex{Float64},1}: + 0.3765784040323032 + 0.09959915134763689im + 0.30865495875240445 + 0.008960297181538185im + 0.4087137042139992 + 0.15906877547743775im + 0.9787481874161135 + 0.0443939782417711im + 0.3430534698620533 + 0.049837687199345705im + -0.2863097014631293 - 0.9011417554715162im + 0.6116671743160434 + 0.14049254864376023im + 0.40933722321954447 + 0.15820580776369225im + 0.43950860634751715 + 0.22808195062035772im + 0.37687009160849716 + 0.09924325053688597im + 0.47944942696208637 + 0.40059913080096726im + 0.4934801111510124 + 0.520295324777746im + 0.496596553706975 + 0.480303769976205im + 0.49307795048742775 + 0.6085434637464817im + 0.5095613980300516 + 0.6639488620533516im + 0.4998069928360967 + 0.3870365082453997im + 0.4861657916302239 + 0.45751028495939855im + 0.4766598864386299 + 0.5557787925858408im + 0.4684227095574837 + 0.4353215518427161im + 0.5014319544500476 + 0.5889845608687214im +``` +All the three eigenvalues in Table 7.1 in Schmid-Henningson are computed this way. diff --git a/docs/src/test.jl b/docs/src/test.jl deleted file mode 100644 index 71f1ad5be..000000000 --- a/docs/src/test.jl +++ /dev/null @@ -1,2 +0,0 @@ -M(λ) -λ \ No newline at end of file From 6bb601c9ecdb298e33adea589ec6a81823ed06aa Mon Sep 17 00:00:00 2001 From: upadp Date: Fri, 9 Aug 2019 15:03:34 +0200 Subject: [PATCH 6/6] Orr Sommerfeld tutorial final version --- docs/make.jl | 2 +- docs/src/hydrotutorial.md | 99 +++++++++++++++++++++++++-------------- 2 files changed, 64 insertions(+), 37 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 6626cf6e1..74635135d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -20,7 +20,7 @@ makedocs( "Tutorial 4 (Python NEP)" => "tutorial_call_python.md", "Tutorial 5 (MATLAB 1)" => "tutorial_matlab1.md", "Tutorial 6 (FORTRAN 1)" => "tutorial_fortran1.md", - "Tutorial 7 (Hydrodynamic stability PEP)" => "hydrotutorial.md" + "Tutorial 7 (Orr-Somerfeld)" => "hydrotutorial.md" ] ) diff --git a/docs/src/hydrotutorial.md b/docs/src/hydrotutorial.md index 69224ede7..624181836 100644 --- a/docs/src/hydrotutorial.md +++ b/docs/src/hydrotutorial.md @@ -1,24 +1,42 @@ # Tutorial: Stability of parallel shear flows ## Background -Stability analysis of flows is a very important problem in fluid mechanics. Of particular interest are applications where a transition from laminar to turbulent flow is observed and this begins with an instability in the laminar flow. Considering a mean base laminar flow $$\overline{U} = \begin{pmatrix}U(y)& 0& 0\end{pmatrix}$$ and the perturbation $$u' = \begin{pmatrix}u& v& w\end{pmatrix}$$. The normal vorticity is defined as $$\eta = \dfrac{\partial u}{\partial z}-\dfrac{\partial w}{\partial x}$$. Linearizing the Navier-Stokes equations around the mean flow and then eliminating pressure gives us the Orr-Sommerfeld and Squire equations, which are a system of fourth order PDEs describing the dynamics of $v$ and $\eta$. - +Stability analysis of flows is a very important problem in fluid mechanics. Linearizing the Navier-Stokes equations around the mean flow and then eliminating pressure gives us the Orr-Sommerfeld and Squire equations, which are a system of fourth order PDEs describing the dynamics: ```math -\left(\Big(\dfrac{\partial }{\partial t}+U\dfrac{\partial }{\partial x}\Big)\nabla^2-U''\dfrac{\partial }{\partial x}-\frac{1}{Re}\nabla^4\right)v = 0,\\ -\left(\dfrac{\partial }{\partial t}+U\dfrac{\partial }{\partial x}-\frac{1}{Re}\nabla^2\right)\eta = -U''\dfrac{\partial v}{\partial z}. +\left(\Big(\dfrac{\partial }{\partial t}+U\dfrac{\partial }{\partial x}\Big)\nabla^2-U''\dfrac{\partial }{\partial x}-\frac{1}{Re}\nabla^4\right)v = 0\\ +\left(\dfrac{\partial }{\partial t}+U\dfrac{\partial }{\partial x}-\frac{1}{Re}\nabla^2\right)\eta = -U''\dfrac{\partial v}{\partial z} ``` +More precisely, these equations stem from modeling using a +mean base laminar flow $$\overline{U} = \begin{pmatrix}U(y)& 0& 0\end{pmatrix}$$ and the perturbation $$u' = \begin{pmatrix}u& v& w\end{pmatrix}$$. The normal vorticity is denoted $$\eta$$. -## Formulation as a nonlinear eigenvalue problem -The ansatz $v(x,y,z,t) = \tilde{v}(y)\exp(i(\alpha x+\beta z-\omega t))$ and $\eta(x,y,z,t) = \tilde{\eta}(y)\exp(i(\alpha x+\beta z-\omega t))$, implying plane wave perturabtions transforms the system to +This is a text-book model for fluid flows. This model is often used to study stability by making a plane wave ansatz and a transformation, and subsequently rewriting the discretized problem as a large standard eigenvalue problem [Chapter 7, Stability and Transition in Shear Flows, Schmid, Peter J., Henningson, Dan S](https://www.springer.com/gp/book/9780387989853#aboutBook). We will here show that the plane wave ansatz directly leads to a NEP, which can be solved with methods in NEP-PACK without any additional transformations. We reproduce computational results in the above reference. + +We thank Miguel Beneitez and Prabal Negi, Department of Mechanics, KTH for the valuable discussions and the discretization code. +## Formulation as a nonlinear eigenvalue problem +To study stability we use the plane wave perturbations ansatz +```math +v(x,y,z,t) = \tilde{v}(y)\exp(i(\alpha x+\beta z-\omega t))\\ +\eta(x,y,z,t) = \tilde{\eta}(y)\exp(i(\alpha x+\beta z-\omega t)), +``` +which transforms the system to ```math \left((-i\omega+i\alpha U)(\mathcal{D}^2-\alpha^2-\beta^2)-i\alpha U''-\frac{1}{Re}(\mathcal{D}^2-\alpha^2-\beta^2)^2\right)\tilde{v} = 0,\\ \left((-i\omega+i\alpha U)-\frac{1}{Re}(\mathcal{D}^2-\alpha^2-\beta^2)\right)\eta = -i\beta U'\tilde{v}, ``` -where $\mathcal{D}$ denotes the 1-D differential operator $\frac{\partial }{\partial y}$. In this example, we consider the boundary conditions $\tilde{v} = \mathcal{D}\tilde{v} = \tilde{\eta} = 0$.We assume that $$\beta$$ and $$\omega$$ are given and we wish to solve for $$\alpha$$. +where $\mathcal{D}$ denotes the 1-D differential operator $\frac{\partial }{\partial y}$. In this example, we consider the boundary conditions $\tilde{v} = \mathcal{D}\tilde{v} = \tilde{\eta} = 0$. We assume that $$\beta$$ and $$\omega$$ are given and we wish to solve for the eigenvalue $$\alpha\in\mathbb{C}$$. -This is usually done by using a transformation of the form $$\begin{pmatrix}\tilde{v}\\ \tilde{\eta}\end{pmatrix} = \begin{pmatrix}\tilde{V}\\ \tilde{E}\end{pmatrix}\exp(-\alpha y)$$ which reduces the power of $$\alpha$$ from four to two. The problem is then discretized and solved as a quadratic eigenvalue problem. See [Chapter 7, Stability and Transition in Shear Flows, Schmid, Peter J., Henningson, Dan S] for details. Since NEP-PACK allows us to tackle polynomial eigenvalue problems of any order, we can simply discretize $$\mathcal{D}$$ to $$D$$ (using a suitable numerical discretization) and expand, which leads to a polynomial eigenvalue problem of fourth order. +This is usually done by using a transformation of the form +```math +\begin{pmatrix}\tilde{v}\\ \tilde{\eta}\end{pmatrix} = \begin{pmatrix}\tilde{V}\\ \tilde{E}\end{pmatrix}\exp(-\alpha y) +``` +which reduces the power of $$\alpha$$ from four to two. The problem is then discretized and solved as a quadratic eigenvalue problem. See Chapter 7 in Schmid-Henningson for details. + +Rather than using the transformation approach described in the +book of +Schmid-Henningson, we +simply discretize $$\mathcal{D}$$ to $$D$$ (using a suitable numerical discretization) which leads to a polynomial eigenvalue problem of fourth order. We define diagonal matrices $$U_0$$, $$U_1$$ and $$U_2$$ as follows. ```math @@ -26,17 +44,18 @@ U_0 = diag(U(y_0),U(y_1),\ldots,U(y_n)),\\ U_1 = diag(U'(y_0),U'(y_1),\ldots,U'(y_n)),\\ U_2 = diag(U''(y_0),U''(y_1),\ldots,U''(y_n)).\\ ``` -Here $$\{y_i\}_{i=1,\ldots,n}$$ denotes the y-coordinates of the $$n$$ grid points used for discretization. The discretized problem on expansion is +Here $$\{y_i\}_{i=1,\ldots,n}$$ denotes the y-coordinates of the $$n$$ grid points used for discretization. The discretized problem after expanding the powers gives ```math \left[-\frac{I}{Re}\alpha^4-iU_0\alpha^3+\left(\frac{2D^2}{Re}+\left(i\omega- -\frac{2\beta^2}{Re}\right)I\right)\alpha^2+\left(iU_0(D^2-\beta^2I)-U_2\right)\alpha +\Big(\frac{2\beta^2D^2}{Re}-\frac{D^4}{Re}-\frac{\beta^4I}{Re}+i\omega(\beta^2I-D^2)\Big)\right]\tilde{v} = 0\\ +\frac{2\beta^2}{Re}\right)I\right)\alpha^2\right.+\\ +\left.\left(iU_0(D^2-\beta^2I)-U_2\right)\alpha +\Big(\frac{2\beta^2D^2}{Re}-\frac{D^4}{Re}-\frac{\beta^4I}{Re}+i\omega(\beta^2I-D^2)\Big)\right]\tilde{v} = 0\\ i\beta U_1\tilde{v}+\left[\frac{I}{Re}\alpha^2 + iU_0 \alpha +\left(\left(\frac{\beta^2}{Re}-i\omega\right)I-\frac{D^2}{Re}\right)\right]\tilde{\eta} = 0 ``` This can be be formulated as a polynomial eigenvalue problem where ```math M(\lambda) = A_0+A_1\lambda+A_2\lambda^2+A_3\lambda^3+A_4\lambda^4, ``` -and the matrices $$A_i,_{i=1,\ldots,n}$$ are given by +and the matrices $$A_0,\ldots,A_4$$ are given by ```math A_0 = \begin{pmatrix}\frac{2\beta^2D^2}{Re}-\frac{D^4}{Re}-\frac{\beta^4I}{Re}+i\omega(\beta^2I-D^2)& 0\\ i\beta U_1&\left(\frac{\beta^2}{Re}-i\omega\right)I-\frac{D^2}{Re}\end{pmatrix}\\ A_1 = \begin{pmatrix}-iU_0& 0\\0& 0\end{pmatrix}\\ @@ -45,7 +64,7 @@ A_2 = \begin{pmatrix}\frac{2D^2}{Re}+\left(i\omega- A_3 = \begin{pmatrix}iU_0(D^2-\beta^2I)-U_2& 0\\0& iU_0\end{pmatrix}\\ A_4 = \begin{pmatrix}-\frac{I}{Re}& 0\\0& 0\end{pmatrix} ``` -## Problem setup +## Problem setup in NEP-PACK We begin by initializing the parameters to the values used to generate the data in Table 7.1 in Schmid-Henningson. ```julia ## Parameters @@ -63,9 +82,9 @@ D2 = DM[2:N+1,2:N+1,2]; #D^2 yF,D4 = cheb4c(N+2); eye = Matrix{Float64}(I, N, N); #D^4 ``` -The code for implementation of `chebdif()` and `cheb4c()` is provided in [GITHUB LINK]. We can now set up the coefficient matrices and create a corresponding PEP object. +The code for implementation of `chebdif()` and `cheb4c()` is provided in [cheb4c.jl](https://nep-pack.github.io/NonlinearEigenproblems.jl/cheb4c.jl) and [chebdif.jl](https://nep-pack.github.io/NonlinearEigenproblems.jl/chebdif.jl). We can now set up the coefficient matrices and create a corresponding PEP object. ```julia -#Setuo coefficient matrices +#Setup coefficient matrices A4 = [-eye/Re zeros(N,N);zeros(N,N) zeros(N,N)]; A3 = [-1im*diagm(0 => U) zeros(N,N);zeros(N,N) zeros(N,N)]; A2 = [(1im*ω-2β^2/Re)*eye+2D2/Re zeros(N,N);zeros(N,N) eye/Re]; @@ -75,47 +94,46 @@ A0 = [2β^2*D2/Re-D4/Re-β^4*eye/Re+1im*ω*(β^2*eye-D2) zeros(N,N);1im*β*diagm #Create a PEP object nep = PEP([A0,A1,A2,A3,A4]); ``` -## Shifting and scaling -Before we proceed to using one of NEP-PACK's methods, we have to consider the issue of ill-conditioned coefficient matrices. We can use the reference eigenvalue $\lambda_0 = 0.97857+0.044394i$ to do a singularity test on $$M(\lambda_0)$$ to see this. -```julia +## Solving with NEP-PACK +Before we proceed to using one of NEP-PACK's methods, we have to consider the issue of poorly scaled coefficient matrices. Hence, direct application of NEP-PACK's solvers on this problem is not adequate. +```julia-repl julia> norm(A0) 8.37194982854379e13 - julia> norm(A1) 473949.06740743306 - julia> norm(A2) 303285.4108872535 - julia> norm(A3) 9.817076958035932 - julia> norm(A4) 0.008 - -julia> minimum(svdvals(compute_Mder(nep,0.97857+1im*0.044394))) -0.00030839497238639225 ``` We can get around this issue by scaling the PEP with NEP-PACK's `shift_and_scale()` , and solving the scaled problem $$T(\lambda) = M(100\lambda)$$ instead. ```julia sc=100; nep1 = shift_and_scale(nep,scale=sc); +mult_scale= norm(nep1.A[end]); +nep2 = PEP(nep1.A ./ mult_scale) ``` -The scaled PEP performs much better in the singularity test. -```julia -julia> minimum(svdvals(compute_Mder(nep2,0.0097857+1im*0.00044394))) -3.8570057996438436e-10 +The scaled PEP has much better scaled coefficient matrices. +```julia-repl +julia> norm.(get_Av(nep2)) +5-element Array{Float64,1}: + 1.0464937285679738e8 + 59.24363342592913 + 3791.0676360906696 + 12.271346197544913 + 1.0 ``` - -## Solution -In this example, we are interested in computing several eigenvalues and our region of interest for the spectrum is in the first quadrant. We use the Tensor Infinite Arnoldi (TIAR) method, as in [REFERENCE TO IAR]. The method is called twice with different shifts `σ`. -```julia -λ1,v1 = tiar(nep2,σ=0.006,v=ones(size(nep,1)),logger=1,neigs=10,maxit=200,tol=1e-14) -λ2,v2 = tiar(nep2,σ=0.005+0.005i,v=ones(size(nep,1)),logger=1,neigs=10,maxit=200,tol=1e-14) +In this example, we are interested in computing several eigenvalues and our region of interest for the spectrum is in the first quadrant. We use the Tensor Infinite Arnoldi (TIAR) method implemented in [`tiar()`](methods.md#NonlinearEigenproblems.NEPSolver.tiar). The method is called twice with different shifts `σ`. +```julia-repl +λ1,v1 = tiar(nep2,σ=0.006,v=ones(size(nep,1)),neigs=10,maxit=200,tol=1e-14) +λ2,v2 = tiar(nep2,σ=0.005+0.005i,v=ones(size(nep,1)),neigs=10,maxit=200,tol=1e-14) +λtotal = [λ1;λ2]; ``` The computed eigenvalues are scaled back to get the eigenvalues of the original problem. ```julia - julia> λ_orig = 100*λtotal + julia> λ_orig = sc*λtotal 20-element Array{Complex{Float64},1}: 0.3765784040323032 + 0.09959915134763689im 0.30865495875240445 + 0.008960297181538185im @@ -138,4 +156,13 @@ The computed eigenvalues are scaled back to get the eigenvalues of the original 0.4684227095574837 + 0.4353215518427161im 0.5014319544500476 + 0.5889845608687214im ``` -All the three eigenvalues in Table 7.1 in Schmid-Henningson are computed this way. +This reproduces all the three eigenvalues in Table 7.1 in Schmid-Henningson. We plot the computed eigenvalues. +```@raw html +
+ +``` +This is in agreement with Figure 7.2 from Schmid-Henningson for the eigenvalues in the first quadrant. +```@raw html +
+ +``` \ No newline at end of file