a flow around a sphere

Question

Hi everybody

I try to model the base flow of a flow around a sphere in 3D with Fenics and for doing that I have an unbelievable problem which you 'll see below.

here is the Mesh first.

  import sys, math

     if len(sys.argv) < 2:
         sys.exit("Must specify Reynolds number; restart is optional")

# Whether to use previous solution as initial condition
 restart = "no"
 if len(sys.argv) == 3:
    restart = sys.argv[2]

# Set parameter values
Re   = float(sys.argv[1])
D    = 0.1
Uinf = 1.0
nu   = D * Uinf / Re

 from dolfin import *
 import numpy as np
 import scipy.linalg as la


 # Constants related to the geometry
  bmarg = 1.e-3 + DOLFIN_EPS
  xmin = -0.2; xmax = 0.5
  ymin = 0.0; ymax = 0.0
  zmin = 0.0; zmax = 0.0
  xcenter = 0.0; ycenter = 0.0; zcenter = 0.0; radius = 0.05

   mesh = Mesh("mesh3d.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
        dz = x[2] - zcenter
        r = sqrt(dx*dx + dy*dy + dz*dz)
       return on_boundary and r < radius + bmarg

 # Define function spaces
 V = VectorFunctionSpace(mesh, "CG", 3)
 Q = FunctionSpace(mesh, "CG", 1)
 W = V*Q
 # no-slip velocity b.c
 noslipBoundary = NoslipBoundary()
 g0 = Constant( (0.0, 0.0, 0.0) )
 bc0 = DirichletBC(W.sub(0), g0, noslipBoundary)

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


 # collect b.c.
 bc = [bc0, bc1,]
 # define function
 # Define test functions
 (v,q) = TestFunctions(W)

 # Define trial functions
  w     = Function(W)
 (u,p) = (as_vector((w[0], w[1], w[1])), w[3])

  # Reynolds number 
   print("Reynolds:", Re)


   parameters["form_compiler"]["quadrature_degree"] = 3
   parameters["form_compiler"]["optimize"] = True
   parameters["form_compiler"]["representation"] = 'quadrature'

   # Weak form
   F =   inner(grad(u)*u, v)*dx        \
           + (1./Re)*inner(grad(u), grad(v))*dx \
           -  p*div(v)*dx                   \
            - q*div(u)*dx

  # Derivative of weak form
   dw = TrialFunction(W)
   dF = derivative(F, w, dw)

   problem = NonlinearVariationalProblem(F, w, bc, dF)
   solver  = NonlinearVariationalSolver(problem)
   # Set linear solver parameters
   itsolver = solver.parameters["newton_solver"]
   itsolver["absolute_tolerance"] = 1.0e-10
   itsolver["relative_tolerance"] = 1.0e-10
   # To see various solver options, uncomment following line
   #info(solver.parameters, True); quit()
   # If you want to initialize solution from previous computation
   if restart == "yes":
      print "Setting initial condition from file ..."
   File("steady.xml") >> w.vector()

   # Solve the problem
    solver.solve()

   # Save steady solution
    File("steady.xml") << w.vector()

   # Save vtk for visualization
   (u,p) = w.split()
   File("velocity.pvd") << u
   File("pressure.pvd") << p


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

   interactive()

and this is the error I get
fall@fall-ENLK13BZ:~/Downloads/a3d$ python newton_3D.py 210.0
('Reynolds:', 210.0)
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 2.577e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-10)

UMFPACK V5.4.0 (May 20, 2009): ERROR: out of memory

Traceback (most recent call last):
File "newton_3D.py", line 107, in
solver.solve()
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 successfully call PETSc function 'KSPSolve'.
*** Reason: PETSc error code is: 76.
*** Where: This error was encountered inside /build/buildd/dolfin-1.4.0+dfsg/dolfin/la/PETScLUSolver.cpp.
*** Process: unknown

*** DOLFIN version: 1.4.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

I'll really appreciate for a help

thanks in advance

a flow around a sphere

1 Answer