Hi !
I’ve solved most of my issues from my previous post but now I’m stuck with weird eigenvectors. Let me explain what I’m trying to achieve.
I’m interested in looking at the stability of my system for a given Reynolds number. In other words, let \dfrac{\partial w}{\partial t} = \mathcal{F}(w) the representation of our Navier-Stokes equations.
We can decompose the state vector into a steady base state and an unsteady perturbation such that w(x,t)=\bar{w}(x)+\epsilon w'(x,t). Inserting this into the previous equation
\qquad \qquad \qquad \qquad \qquad \dfrac{\partial \bar{w}}{\partial t} + \dfrac{\partial w'}{\partial t} = \mathcal{F}(\bar{w}+w')
The base state is steady so the time derivative can be dropped. In our case, \mathcal{F} is a non-linear operator so we need to linearize it around the base state.
\qquad \qquad \qquad \qquad \qquad \mathcal{F}(\bar{w} + w')=\mathcal{F}(\bar{w})+\left. \dfrac{\partial \mathcal{F}(w)}{\partial w} \right\rvert_{\bar{w}}w' + \mathcal{O}(w'^2)
With \mathcal{F}(\bar{w}) the solution of the steady Navier-Stokes equations (i.e. \mathcal{F}(\bar{w})=0). By considering the previous equation, neglecting the higher order terms and writing \mathcal{J}(\bar{w})= \left. \dfrac{\partial \mathcal{F}(w)}{\partial w} \right\rvert_{\bar{w}} we obtain
\qquad \qquad \qquad \qquad \qquad \dfrac{\partial w'}{\partial t} = \mathcal{J}(\bar{w})w'
Where \mathcal{J}(\bar{w}) is the Jacobian matrix. The perturbations are solved by characterising them as normal modes w'=\hat{w}e^{\sigma t}. Putting this back in the last equation
\qquad \qquad \qquad \qquad \qquad \sigma \hat{w}= \mathcal{J}(\bar{w})\hat{w}
With \sigma and \hat{w} respectively the eigenvalues and eigenvectors of the system. For recall, if one eigenvalue has a positive real part (>0) then the system is absolutely unstable.
Now the implementation (MWE in the end), I’m considering as well known stable test case (i.e. all eigenvalues should be negative until Re = 47)
Computation of the baseflow OK and verified, assembly of the Jacobian UNSURE, solving the eigenproblem OK but weird eigenvalues/vectors when plotting.
For the Jacobian assembly I have two options:
- Replace my BC with homogeneous Dirichlet in order to have a well defined weak form, problem is I have different positive eigenvalues and weird eigenvectors
- Ignore BC, no high eigenvalues comming from the BC but still positive (should be all negative) and still weird eigenvectors
Most of the people studying Navier-Stokes stability are consedering linearized Navier-Stokes equations in order to assemble an approximation of the Jacobian. I would like to take advantage of FEniCS and the derivative function.
So far, I know that the baseflow calculation is working properly. I have some doubt on the derivative function and the way that I’m assembling the matrix since the eigenvectors are really weird (almost zero everywhere).
Meshfile needed to run the code:
https://drive.google.com/file/d/1gQWqMFP3Mk_5I9Ym6P20UrI3gvEW1m6l/view?usp=sharing
MWE:
from dolfin import *
import matplotlib.pyplot as plt
comm = MPI.comm_world
solfile = XDMFFile(comm, "BF_RealEigenvector.pvd")
solfile.parameters["rewrite_function_mesh"] = False
solfile.parameters["functions_share_mesh"] = True
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["std_out_all_processes"] = False
#Apply or not BC to the jacobian matrix
WithBC = False
nev = 50
Re = 40
nu = 1 / Re
mesh_file = "cylinder_coarse.h5"
#Read the mesh and the facet, 1 = left+top+bot, 2 = cylinder
mesh = Mesh()
try:
with HDF5File(comm, mesh_file, "r") as h5file:
h5file.read(mesh, "mesh", False)
facet = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
h5file.read(facet, "facet")
except FileNotFoundError as fnf_error:
print(fnf_error)
#Taylor-Hood element [P2, P2]-P1
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, V * Q)
#All the test/trial functions needed to solve the problem
v, q = TestFunctions(W)
w = Function(W)
u, p = split(w)
u_export, p_export = w.split()
u_export.rename("Velocity", "0")
p_export.rename("Pressure", "1")
#BC
inflow = Constant((1.0, 0.0))
noslip = Constant((0.0, 0.0))
inflowBC = DirichletBC(W.sub(0), inflow, facet, 1)
noslipBC = DirichletBC(W.sub(0), noslip, facet, 2)
bcs = [inflowBC, noslipBC]
#Weak form of Steady Navier-Stokes
F = (Constant(nu) * inner(grad(u), grad(v))\
+ dot(dot(grad(u), u), v)\
- p * div(v)\
- q * div(u)) * dx
#Creating the solver object
J = derivative(F, w)
problem = NonlinearVariationalProblem(F, w, bcs, J)
solver = NonlinearVariationalSolver(problem)
#Solving the problem
with Timer() as t:
begin("BaseFlow computation at Re = %d" %Re)
solver.solve()
end()
info("Computation duration = %f s" %t.elapsed()[0])
solfile.write(u_export, 0)
solfile.write(p_export, 0)
#So far no problem, the baseflow looks fine
#Now studying the stability of the system around the baseflow
if WithBC:
#New BC, homogeneous Dirichlet everywhere
inflowBC = DirichletBC(W.sub(0), noslip, facet, 1)
noslipBC = DirichletBC(W.sub(0), noslip, facet, 2)
bcs = [inflowBC, noslipBC]
else:
bcs = []
A = PETScMatrix()
#Assembling the jacobian of the weak form evaluated at
#the baseflow previously calculated and stored in w
J = derivative(F, w)
dummy = inner(Constant((1.0, 1.0)), v) * dx
assemble_system(J, dummy, bcs, A_tensor = A)
#Solving the standard eigenvalue problem, looking for the
#largest real part
eigensolver = SLEPcEigenSolver(A)
PETScOptions.set("eps_view")
eigensolver.parameters["spectrum"] = "largest real"
eigensolver.solve(nev)
#Exporting the real part of the eigenvectors and plotting
#eigenvalues
Real = []
Complex = []
for i in range(nev):
r, c, rv, cv = eigensolver.get_eigenpair(i)
w.vector()[:] = rv
solfile.write(u_export, i + 1)
solfile.write(p_export, i + 1)
Real.append(r)
Complex.append(c)
solfile.close()
plt.plot(Real, Complex, "ro")
plt.xlabel("Real part")
plt.ylabel("Complex part")
plt.show()