Solver fails to solve the problem with viscosity nu = 0.001 [m^2 / s]. Simple 3D Flow in Square Pipe

Question

Hi,

I am trying to solve a basic 3D Navier-Stokes equations using FEniCS.

There are many examples in the internet of solving Navier-Stokes equations using FEniCS.
But almost all in 2D or in 3D with high viscosity (like 1 or 1/8).

I easily solved 3D Flow in Square Pipe with high viscosity nu = (1/8) [m^2 / s].

But when I tried to solve the solution with viscosity nu = 0.001 [m^2 / s], the solver crashed.
I tried to increase mesh size, decrease time step and change initial velocity/pressure, but it wasn't help.
I also tried to solve the problem with different boundary conditions:
Inlet velocity and outlet pressure or Inlet/Outlet pressure.
On the walls I set NoSlipBoundary condition of course.
The solution fails in both cases.

My Fenics/dolfin code for 3D Flow in Square Pipe I write below.

Box_chorin_nu1000_viscosity_u_test.py
(The case with Inlet velocity and outlet pressure boundary conditions)

from dolfin import *

set_log_level(WARNING)

# No-slip boundary
class NoslipBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (near(x[0], 0) or near(x[1], 0) or near(x[0], 1) or near(x[1], 1))

# Inflow boundary
class InflowBoundary(SubDomain):
    def inside(self, x, on_boundary):        
        return near(x[2], 0) and on_boundary

# Outflow boundary
class OutflowBoundary(SubDomain):
    def inside(self, x, on_boundary):        
        return near(x[2], 10) and on_boundary

#---------------------------------------------------------------------

def Box_chorin():  

  mesh = BoxMesh(0, 0, 0, 1, 1, 10, 12, 12, 60)
  h = mesh.hmin()

  #! problem specific
  f = Constant((0., 0., 0.))         # force
  nu = Constant(1./1000.)           # kinematic viscosity  
  #nu = Constant(1./8.) 
  dt = 0.2*h/0.3                 # time step CFL with 1 = max. velocity
  k = Constant(dt)               # time step 
  T = 20.0                       # total simulation time
  u0 = Constant((0., 0., 0.))        # initial velocity
  p0 = Constant(0.)             # initial pressure

  #! solver specific
  V = VectorFunctionSpace(mesh, 'CG', 2)
  Q = FunctionSpace(mesh, 'CG', 1)

  u = TrialFunction(V)
  v = TestFunction(V)
  p = TrialFunction(Q)
  q = TestFunction(Q)

  u0 = interpolate(u0, V)  
  p0 = interpolate(p0, Q)
  us = Function(V)
  u1 = Function(V)
  p1 = Function(Q)

  # tentative velocity
  F0 = (1./k)*inner(u - u0, v)*dx + inner(dot(grad(u0), u0), v)*dx\
       + nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
  a0, L0 = system(F0)

  # projection
  F1 = inner(grad(p), grad(q))*dx + (1./k)*q*div(us)*dx
  a1, L1 = system(F1)

  # finalize
  F2 = (1./k)*inner(u - us, v)*dx + inner(grad(p1), v)*dx 
  a2, L2 = system(F2)

  # boundary conditions
  b_v = DirichletBC(V, Constant((0.0, 0.0, 0.0)), NoslipBoundary())
  b_v1 = DirichletBC(V, Constant((0.0, 0.0, 0.1)), InflowBoundary())  
  b_p1 = DirichletBC(Q, Constant(0.),  OutflowBoundary())
  bcs_v = [b_v, b_v1]
  bcs_p = [b_p1]

  A0 = assemble(a0)
  A1 = assemble(a1)
  A2 = assemble(a2)

  solver02 = KrylovSolver('gmres', 'ilu')
  solver1 = KrylovSolver('cg', 'petsc_amg')

  ufile = File("results_Box_nu1000_viscosity_u_test/Box_chorin_nu1000_u_ns-fenics-master_u.pvd")
  pfile = File("results_Box_nu1000_viscosity_u_test/Box_chorin_nu1000_u_ns-fenics-master_p.pvd")

  t = 0
  while t < T:
    t += dt

    b = assemble(L0)
    [bc.apply(A0, b) for bc in bcs_v]
    solver02.solve(A0, us.vector(), b)

    b = assemble(L1)
    [bc.apply(A1, b) for bc in bcs_p]
    solver1.solve(A1, p1.vector(), b)

    b = assemble(L2)
    [bc.apply(A2, b) for bc in bcs_v]
    solver02.solve(A2, u1.vector(), b)

    ufile << u1
    pfile << p1

    u0.assign(u1)
    print t

  return u0, p1

#----------------------------------------------------------------------

if __name__ == '__main__':
  u, p = Box_chorin()

Box_chorin_nu1000_viscosity_p_test.py
(The case with Inlet/Outlet pressure boundary conditions)

  # boundary conditions
  b_v = DirichletBC(V, Constant((0.0, 0.0, 0.0)), NoslipBoundary())
  b_p0 = DirichletBC(Q, Constant(0.032), InflowBoundary())
  b_p1 = DirichletBC(Q, Constant(0.),  OutflowBoundary())
  bcs_v = [b_v]
  bcs_p = [b_p0, b_p1]

Also I solved the problem in another CFD package, and it solved the problem in both cases (viscosity nu = (1/8) [m^2 / s] and viscosity nu = 0.001 [m^2 / s]) without any errors.

In Fenics I use Chorin’s method for solving Incompressible Navier-Stokes equations, like in the Fenics tutorials here:
http://fenicsproject.org/documentation/dolfin/1.3.0/python/demo/documented/navier-stokes/python/documentation.html

For solving problem I use this Fenics application (ns-fenics-master), which is based on nsbench Fenics application:
https://github.com/MiroK/ns-fenics

Other test and demo examples work without any errors so I don't even suggest what may cause the solver fail.

Could you hint, help to notice or suggest why the solver fails to solve the problem with viscosity nu = 0.001 [m^2 / s]?

Even ordinary water has nu = 1^(-6) or 1e-6 [m^2 / s] for temperature T = 20 C and it is far from nu = (1/8) [m^2 / s].

Any help appreciated!

Thanks in advance!

Regards, Maksim

Answer 1

This program runs fine as is on my development version of dolfin. No crashing, nice and stable.

Btw, you should not use "cg" solver on pressure when you apply Dirichlet boundary conditions. The matrix becomes unsymmetric unless you assemble with assemble_system. Furthermore, you could solve this problem in 2D, or you can use periodic boundary conditions in the z-direction (use a constant pressure gradient as forcing).

Comments

Thank you for suggestion!

I am using Fenics dolfin version 1.3 under Ubuntu OS.

I change solver for pressure from 'cg' to 'gmres':

#solver1 = KrylovSolver('cg', 'petsc_amg')
solver1 = KrylovSolver('gmres', 'petsc_amg')

but the problem also fails on T = 45.5877257682.
The solution almost convergent and then solver fails.
The result before the fail similarly to a result is observed in other CFD package, besides the crash of course.

The error:

Traceback (most recent call last):
  File "Box_chorin_nu1000_viscosity_u_test.py", line 116, in <module>
    u, p = Box_chorin()
  File "Box_chorin_nu1000_viscosity_u_test.py", line 93, in Box_chorin
    solver02.solve(A0, us.vector(), b)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = inf).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 1.3.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

I am also increase solution time from T = 20.0 to T = 200.0 to figure out the error.

Also one can add "iter" to decrease amount of output results.

iter = 0
t = 0
while t < T:
  t += dt
  iter += 1

  if ((iter % 20) == 0):
      print iter 
      ufile << u1
      pfile << p1

---

Are you after a steady solution? Chorin is a very poor steady state solver. Furthermore, your convection term is fully explicit and thus not very stable. If you're after a steady solution I recommend using a fully coupled instead of a fractional step solver that is only intended for unsteady flow.

---

"Are you after a steady solution?"
Yes, I try to get a steady solution for this problem.

"Chorin is a very poor steady state solver. Furthermore, your convection term is fully explicit and thus not very stable."
Yes, it is.

"If you're after a steady solution I recommend using a fully coupled instead of a fractional step solver that is only intended for unsteady flow."
Could you hint where I can find an example of a fully coupled Naviers-Stokes solver (like solvers in nsbench project)?

Or it is like in the Fenics dolfin tutorials below?
(from http://fenicsproject.org/documentation/dolfin/1.3.0/python/demo/index.html)
18. Stokes equations
19. Stokes equations with Mini elements
20. Stokes equations with stabilized first order elements
21. Stokes equations with Taylor-Hood elements

Thanks in advance!

---

These tutorials are a good start. Just add the nonlinear convection and solve using Newton iterations. See a nonlinear demo.