This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Eigenvalue for navier stokes equation

0 votes

Hi !

I work on a project of instability analysis corcernant a cylinder in 2D but I am having a weird prolème on the mesh . When it is refined as FreeFem ++ it does not work on me FENICS ( I do not get the values ​​), while when I reduced so that the mesh is fine. Please look at the attached code.
Thank you in advance.

from dolfin import *
from mshr import *

# Constants related to the geometry

  xmin = -10.; xmax = 40.
  ymin = -10.0; ymax = 10.
  xcenter = 0.0; ycenter = 0.0; radius = 0.5
  # create domain mesh
  domain = Rectangle(Point(xmin,ymin), Point(xmax,ymax)) - Circle(Point(0.0,0.0),0.5)
  domain.set_subdomain(1, Rectangle(Point(xmin,-3.), Point(xmax,3.)))
  domain.set_subdomain(2, Rectangle(Point(-2.,-1.5), Point(7., 1.5)))
  mesh = generate_mesh(domain, 50)
   # mesh = refine(mesh)
   # refine mesh TWICE around the cylinder
   for refinements in [0,1]:
        cell_markers = CellFunction("bool", mesh)
        cell_markers.set_all(False)
        for cell in cells(mesh):
             p = cell.midpoint()
             if (xmin < p[0] < xmax) and \
                (-2.75 < p[1] < 2.75 ) :
                  cell_markers[cell] = True
        mesh = refine(mesh, cell_markers)

# refine mesh tree times around the cylinder
for refinements in [0,1]:
     cell_markers = CellFunction("bool", mesh)
     cell_markers.set_all(False)
     for cell in cells(mesh):
          p = cell.midpoint()
          if (-2. < p[0] < 7.) and \
             (-1.25 < p[1] < 1.25 ) :
              cell_markers[cell] = True
      mesh = refine(mesh, cell_markers)

  File("mesh.xml.gz") << mesh
  # Plot the mesh
   plot(mesh)
  interactive()

Plus the eigenvalue code ....

from dolfin import *


# Constants related to the geometry
 bmarg = 1.e-3 + DOLFIN_EPS
 xmin = -10.0; xmax = 40.0
 ymin = -10.0; ymax = 10.0
 xcenter = 0.0; ycenter = 0.0; radius = 0.5

  mesh = Mesh("mesh.xml.gz")
  plot(mesh, title = "Mesh domain")

  #boundary conditions using a mesh function
   boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

  # Inflow boundary
  class InflowBoundary(SubDomain):
        def inside(self, x, on_boundary):
               return on_boundary and x[0] < xmin + bmarg

  # No-slip boundary
   class NoslipBoundary(SubDomain):
         def inside(self, x, on_boundary):
               dx = x[0] - xcenter
               dy = x[1] - ycenter
                r = sqrt(dx*dx + dy*dy)
              return on_boundary and r < radius + bmarg

  # Define function spaces
  V = VectorFunctionSpace(mesh, "CG", 2)
  Q = FunctionSpace(mesh, "CG", 1)
  W = V*Q


  # no-slip velocity b.c
  noslipBoundary = NoslipBoundary()
   g0 = Constant( (0., 0.) )
   bc0 = DirichletBC(W.sub(0), g0, noslipBoundary)

  # inlet velocity b.c
  inflowBoundary = InflowBoundary()
  g1 = Constant( (1., 0.) )
  bc1 = DirichletBC(W.sub(0), g1, inflowBoundary)

         # collect b.c.
  bcs = [bc0, bc1]
  # define function
  u = Function(V)
  p = Function(Q)
  Re = 50.0
  print("Reynolds:", Re)
  # Define variational problem for initial guess (stokes solution)
  v, q = TestFunctions(W)
  u, p = TrialFunctions(W)

  def sigma(v):
      return sym(grad(v))

  f = Constant( (0.0, 0.0) )
  F = (inner(sigma(u), grad(v)) - p*div(v) + q*div(u) )*dx + \
         inner(f, v)*dx
  a = lhs(F);  L = rhs(F)
  # Assemble system
  A, b = assemble_system(a, L, bcs)
  up = Function(W)
  solve(A, up.vector(), b, 'lu')
   u_a, p_a = up.split()
  plot(p_a, title = "initial guess pressure")
  plot(u_a, title = "initial guess velocity")
  # Save solution to file
  file = File("initial_guess.pvd") 
  file << up
  # Boundary conditions
  # inlet velocity b.c
  g11 = Constant( (0., 0.) )
  bc11 = DirichletBC(W.sub(0), g11, inflowBoundary)

  # collect b.c.
  bcs1 = [bc0, bc11]

  # define function
  du = Function(V)
  dp = Function(Q)

  # Define variational problem for convergent solution up1
  du, dp = TrialFunctions(W)
  F1 = (inner(grad(du)*u_a, v) + inner(grad(u_a)*du, v) )*dx + \
          (1/Re)*inner(grad(du), grad(v))*dx + \
          (- dp*div(v) + q*div(du))*dx + \
          (inner(grad(u_a)*u_a, v))*dx + (1/Re)*inner(grad(u_a), grad(v))*dx + \
          (- p_a*div(v) + q*div(u_a))*dx
  a1 = lhs(F1);  L1 =  rhs(F1)
         # Newton iteration
 err = 1.0
 tol = 2.e-6
 iter = 0
 maxiter = 10
 up1 = Function(W)
 u1, p1 = up1.split()
 while err > tol and iter < maxiter:
         A1, b1 = assemble_system(a1, L1, bcs1)
         solve(A1, up1.vector(), b1, 'lu') 
         err = norm(up1, 'L2')/25200
         print " iteration=%g, res=%g", iter, err
         iter += 1
         up.vector()[:] += up1.vector()
   # save
   File("base_flow.pvd") << up
   File("velocity.pvd") << u_a
   File("pressure.pvd") << p_a

 U0, P0 = up.split()  # base flow
 plot(P0, title = "base flow pressure")
 plot(U0, title = "base flow velocity")

 ############ normal modes problems u = U0 + du1exp(iwt), #######
u2, p2 = TrialFunctions(W)

a2 = (inner(grad(u2)*U0, v) + inner(grad(U0)*u2, v) + (1/Re)*inner(grad(u2), grad(v))*dx-\
        ( p2*div(v) + q*div(u2) )*dx 
m = - inner(u2, v)*dx 
A0 = PETScMatrix()
assemble(a2, tensor=A0)
M = PETScMatrix()
assemble(m, tensor=M)
for bc in bcs:
    (bc.apply(A0) and bc.apply(M))

 #create eigensolver
 eigensolver = SLEPcEigenSolver(A0, M)
 # eigensolver.parameters['spectrum'] = 'smallest magnitude'
 eigensolver.parameters["spectral_transform"] = "shift-and-invert"
 eigensolver.parameters["spectral_shift"] = 4
 eigensolver.parameters['verbose'] = True

 print 'solving'
 num = 6
 eigensolver.solve(num)
 print 'solving termine'
  for i in range(num):
      r, c, rx, cx = eigensolver.get_eigenpair(i)
      print("Eigenvalue:", i, r, c)

 # Initialize function and assign eigenvector
 Utild = Function(W)
 Utild.vector()[:] = rx
 utild, ptild = Utild.split()
     # Plot eigenfunction
  plot(utild)
  plot(ptild)
  interactive()
asked Dec 14, 2014 by fall FEniCS Novice (230 points)

Hello. It is no fun to answer questions of type "I don't get the values I want, here is the code." Please specify:

  1. What is the problem you want to solve?
  2. What results do you expect?
  3. When and where does your code fail?

1 Answer

0 votes

Hi
yes... the aim of this project is to get the eigenvalues for a given Reynolds number but my problem is with this Mesh, I use above it does not walk but if I change the dimension of the mesh like this it will be fine .

from dolfin import *
from mshr import *

# Constants related to the geometry
bmarg = 1.e-3 + DOLFIN_EPS
xmin = -8; xmax = 13
ymin = -4.0; ymax = 4
xcenter = 0.0; ycenter = 0.0; radius = 0.5
# create domain mesh
domain = Rectangle(xmin,ymin,xmax,ymax) - Circle(0.0,0.0,0.5)
domain.set_subdomain(1, Rectangle(-4, -1.5, 4,1.5)
domain.set_subdomain(2, Rectangle(-1.5, -1.0, 1.5, 1))
mesh =  generate_mesh(domain, 10) etc ...

if I use this Mesh for the second code which calculate the eigenvalues I will get them but you notice that this mesh is not well refined that why I want to use the first one for more accuracy. I think the main problem is the mesh and my goal is to get the eigenvalues in the parts of the second code.

Thanks for your answer and if it' s not clear enough let me know I wil try to make details about this

answered Dec 15, 2014 by fall FEniCS Novice (230 points)

Did you solve this problem ? I am also trying eigenvalue computation but for me the error says

[0]PETSC ERROR: Detected zero pivot in LU factorization: see
http://www.mcs.anl.gov/petsc/documentation/faq.html#ZeroPivot!
[0]PETSC ERROR: Zero pivot row 60138 value 0 tolerance 2.22045e-14!

Since mass matrix is singular, how does the eigensolver work in this case ?

Try applying a shift.

could you explain me how to do it please? if it 's different to the shift and invert I used

For cylinder inside a channel, I have an example code for eigenvalue computation using scipy here https://code.google.com/p/cfdlab/source/browse/trunk/fenics/2d/ns_cylinder/

The slepc eigensolver does not always work for me, and gives zero pivot error.

For cylinder in external flow problem, I am not able to get the eigenfunctions yet. Did you make any progress.I have one code here https://code.google.com/p/cfdlab/source/browse/trunk/fenics/2d/cylinder_control/

hi
I try this but I have a new error ( NotImplementedError: shifted eigenproblem not supported yet) if you have any idea. thanks

Maybe your scipy version is too old. Try upgrading scipy. I am using default scipy on OSX Mavericks and it works on that.

yes it's true i think the same too but i will test it next week in lab. I would like to know if you have other codes in instability in 2D or 3D. I like to carry on on after my project. thanks

...