2D Stationary Navier-Stokes equation

1

Hello,

I am using FEniCS for the first time, for a mathematics project at university. I have to modelised an elliptic equation with FEniCS, and I chose the bidimensional flow in a conduct, with Navier-Stokes equation. I wrote this code, but when solving, an error is raised because of an already indexed object …

My code :

from dolfin import *
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

nu=1E-5
tol=1E-14
X=10
Y=2

from mshr import *
nx=100
ny=20
mesh2 = RectangleMesh(Point(0.0, -1.0), Point(X, Y-1.0), nx, ny, "left")

Hh2 = VectorFunctionSpace(mesh2,"P",2)
Q = FunctionSpace(mesh2, "P", 1)
W = Hh2 or Q


class Bord(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[1] - 1) < tol) or (abs(x[1] + 1) < tol))

class Entree(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (abs(x[O] - 0))

class Sortie(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (abs(x[0] - 10) < tol)

b_haut_bas = Bord()
bhbD = DirichletBC(Hh2, (0, 0), b_haut_bas)

b_in = Entree()
b_ex = Sortie()
inD = DirichletBC(Q, Constant(10), b_in)
exD = DirichletBC(Q, Constant(8), b_ex)

bcu = [bhbD,inD,exD]

class BordFlux(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[0] - 0) < tol) or (abs(x[0] - 10) < tol))
    
bfN = BordFlux()
boundaries = MeshFunction("size_t", mesh2, mesh2.topology().dim()-1, 0)
bfN.mark(boundaries, 0)
ds = Measure("ds", domain=mesh2, subdomain_data=boundaries)

v, q = TestFunctions(W)
up = TrialFunction(W)
u, p = split(up)

a1 = nu*inner(grad(u), grad(v))*dx - inner(p, div(v))*dx - inner(q, div(u))*dx
l1 = 0


up = Function(W)
u, p = split(up)
solve(a1 == l1, up, bcu)

Then the error raised is "UFLException: Attempting to index with (Ellipsis, Index(17)), but object is already indexed: ":

Calling FFC just-in-time (JIT) compiler, this may take some time.
Attempting to index with (Ellipsis, Index(17)), but object is already indexed: <Indexed id=140295614367928>
---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
<ipython-input-100-78ff56ec39b2> in <module>
      6 u, p = split(up)
      7 
----> 8 solve(a1 == l1, up, bcu)

/usr/lib/python3/dist-packages/dolfin/fem/solving.py in solve(*args, **kwargs)
    218     # tolerance)
    219     elif isinstance(args[0], ufl.classes.Equation):
--> 220         _solve_varproblem(*args, **kwargs)
    221 
    222     # Default case, just call the wrapped C++ solve function

/usr/lib/python3/dist-packages/dolfin/fem/solving.py in _solve_varproblem(*args, **kwargs)
    259         # Create problem
    260         problem = NonlinearVariationalProblem(eq.lhs, u, bcs, J,
--> 261                                               form_compiler_parameters=form_compiler_parameters)
    262 
    263         # Create solver and call solve

/usr/lib/python3/dist-packages/dolfin/fem/problem.py in __init__(self, F, u, bcs, J, form_compiler_parameters)
     88 
     89         # Wrap forms
---> 90         F = Form(F, form_compiler_parameters=form_compiler_parameters)
     91         if J is not None:
     92             J = Form(J, form_compiler_parameters=form_compiler_parameters)

/usr/lib/python3/dist-packages/dolfin/fem/form.py in __init__(self, form, **kwargs)
     42 
     43         ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
---> 44                            mpi_comm=mesh.mpi_comm())
     45         ufc_form = cpp.fem.make_ufc_form(ufc_form[0])
     46 

/usr/lib/python3/dist-packages/dolfin/jit/jit.py in mpi_jit(*args, **kwargs)
     45         # Just call JIT compiler when running in serial
     46         if MPI.size(mpi_comm) == 1:
---> 47             return local_jit(*args, **kwargs)
     48 
     49         # Default status (0 == ok, 1 == fail)

/usr/lib/python3/dist-packages/dolfin/jit/jit.py in ffc_jit(ufl_form, form_compiler_parameters)
     95     p.update(dict(parameters["form_compiler"]))
     96     p.update(form_compiler_parameters or {})
---> 97     return ffc.jit(ufl_form, parameters=p)
     98 
     99 

/usr/lib/python3/dist-packages/ffc/jitcompiler.py in jit(ufl_object, parameters, indirect)
    215 
    216     # Inspect cache and generate+build if necessary
--> 217     module = jit_build(ufl_object, module_name, parameters)
    218 
    219     # Raise exception on failure to build or import module

/usr/lib/python3/dist-packages/ffc/jitcompiler.py in jit_build(ufl_object, module_name, parameters)
    131                                     name=module_name,
    132                                     params=params,
--> 133                                     generate=jit_generate)
    134     return module
    135 

/usr/lib/python3/dist-packages/dijitso/jit.py in jit(jitable, name, params, generate, send, receive, wait)
    163         elif generate:
    164             # 1) Generate source code
--> 165             header, source, dependencies = generate(jitable, name, signature, params["generator"])
    166             # Ensure we got unicode from generate
    167             header = as_unicode(header)

/usr/lib/python3/dist-packages/ffc/jitcompiler.py in jit_generate(ufl_object, module_name, signature, parameters)
     64 
     65     code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
---> 66             prefix=module_name, parameters=parameters, jit=True)
     67 
     68     # Jit compile dependent objects separately,

/usr/lib/python3/dist-packages/ffc/compiler.py in compile_form(forms, object_names, prefix, parameters, jit)
    141     """This function generates UFC code for a given UFL form or list of UFL forms."""
    142     return compile_ufl_objects(forms, "form", object_names,
--> 143                                prefix, parameters, jit)
    144 
    145 

/usr/lib/python3/dist-packages/ffc/compiler.py in compile_ufl_objects(ufl_objects, kind, object_names, prefix, parameters, jit)
    183     # Stage 1: analysis
    184     cpu_time = time()
--> 185     analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
    186     _print_timing(1, time() - cpu_time)
    187 

/usr/lib/python3/dist-packages/ffc/analysis.py in analyze_ufl_objects(ufl_objects, kind, parameters)
     88         # Analyze forms
     89         form_datas = tuple(_analyze_form(form, parameters)
---> 90                            for form in forms)
     91 
     92         # Extract unique elements accross all forms

/usr/lib/python3/dist-packages/ffc/analysis.py in <genexpr>(.0)
     88         # Analyze forms
     89         form_datas = tuple(_analyze_form(form, parameters)
---> 90                            for form in forms)
     91 
     92         # Extract unique elements accross all forms

/usr/lib/python3/dist-packages/ffc/analysis.py in _analyze_form(form, parameters)
    172                                       do_apply_geometry_lowering=True,
    173                                       preserve_geometry_types=(Jacobian,),
--> 174                                       do_apply_restrictions=True)
    175     elif r == "tsfc":
    176         try:

/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
    262     # reducing the number of operators later algorithms and form
    263     # compilers need to handle
--> 264     form = apply_algebra_lowering(form)
    265 
    266     # After lowering to index notation, remove any complex nodes that

/usr/lib/python3/dist-packages/ufl/algorithms/apply_algebra_lowering.py in apply_algebra_lowering(expr)
    184     """Expands high level compound operators (e.g. inner) to equivalent
    185     representations using basic operators (e.g. index notation)."""
--> 186     return map_integrand_dags(LowerCompoundAlgebra(), expr)

/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py in map_integrand_dags(function, form, only_integral_type, compress)
     56 def map_integrand_dags(function, form, only_integral_type=None, compress=True):
     57     return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
---> 58                           form, only_integral_type)

/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py in map_integrands(function, form, only_integral_type)
     37     if isinstance(form, Form):
     38         mapped_integrals = [map_integrands(function, itg, only_integral_type)
---> 39                             for itg in form.integrals()]
     40         nonzero_integrals = [itg for itg in mapped_integrals
     41                              if not isinstance(itg.integrand(), Zero)]

/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py in <listcomp>(.0)
     37     if isinstance(form, Form):
     38         mapped_integrals = [map_integrands(function, itg, only_integral_type)
---> 39                             for itg in form.integrals()]
     40         nonzero_integrals = [itg for itg in mapped_integrals
     41                              if not isinstance(itg.integrand(), Zero)]

/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py in map_integrands(function, form, only_integral_type)
     44         itg = form
     45         if (only_integral_type is None) or (itg.integral_type() in only_integral_type):
---> 46             return itg.reconstruct(function(itg.integrand()))
     47         else:
     48             return itg

/usr/lib/python3/dist-packages/ufl/algorithms/map_integrands.py in <lambda>(expr)
     55 
     56 def map_integrand_dags(function, form, only_integral_type=None, compress=True):
---> 57     return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
     58                           form, only_integral_type)

/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py in map_expr_dag(function, expression, compress)
     35     Return the result of the final function call.
     36     """
---> 37     result, = map_expr_dags(function, [expression], compress=compress)
     38     return result
     39 

/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py in map_expr_dags(function, expressions, compress)
     84                 r = handlers[v._ufl_typecode_](v)
     85             else:
---> 86                 r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
     87 
     88             # Optionally check if r is in rcache, a memory optimization

/usr/lib/python3/dist-packages/ufl/algorithms/apply_algebra_lowering.py in div(self, o, a)
    150     def div(self, o, a):
    151         i = Index()
--> 152         return a[..., i].dx(i)
    153 
    154     def nabla_div(self, o, a):

/usr/lib/python3/dist-packages/ufl/indexed.py in __getitem__(self, key)
    122 
    123     def __getitem__(self, key):
--> 124         error("Attempting to index with %s, but object is already indexed: %s" % (ufl_err_str(key), ufl_err_str(self)))

/usr/lib/python3/dist-packages/ufl/log.py in error(self, *message)
    170         "Write error message and raise an exception."
    171         self._log.error(*message)
--> 172         raise self._exception_type(self._format_raw(*message))
    173 
    174     def begin(self, *message):

UFLException: Attempting to index with (Ellipsis, Index(17)), but object is already indexed: <Indexed id=140295614367928>Blockquote

Thank you for your help
It seems that the error only concern the resolution, I noticed no error elsewhere, and the problem correctly works with a fixed gradient of pressure, and an unique component v of the velocity. I decided to solve the vectorial problem with velocity=(u, v) and the pressure p, in ordre to correctly see the flow conservation.

Hyanoa

2

There are many problems with your code, including several typos.
Here is a code that compiles, where I use a MixedElement to define the function space.
Also there are typos in your class Entree, and in your variational problem in solve you have issues as well.

from dolfin import *
from fenics import *

nu=1E-5
tol=1E-14
X=10
Y=2

from mshr import *
nx=100
ny=20
mesh2 = RectangleMesh(Point(0.0, -1.0), Point(X, Y-1.0), nx, ny, "left")
Ve = VectorElement("P",triangle, 2)
Qe = FiniteElement("P",triangle, 1)
W = FunctionSpace(mesh2, MixedElement([Ve, Qe]))


class Bord(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[1] - 1) < tol) or (abs(x[1] + 1) < tol))

class Entree(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (abs(x[0] - 0))

class Sortie(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (abs(x[0] - 10) < tol)

b_haut_bas = Bord()
bhbD = DirichletBC(W.sub(0), (0, 0), b_haut_bas)

b_in = Entree()
b_ex = Sortie()
inD = DirichletBC(W.sub(1), Constant(10), b_in)
exD = DirichletBC(W.sub(1), Constant(8), b_ex)

bcu = [bhbD,inD,exD]

class BordFlux(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[0] - 0) < tol) or (abs(x[0] - 10) < tol))

bfN = BordFlux()
boundaries = MeshFunction("size_t", mesh2, mesh2.topology().dim()-1, 0)
bfN.mark(boundaries, 0)
ds = Measure("ds", domain=mesh2, subdomain_data=boundaries)

v, q = TestFunctions(W)
u, p = TrialFunctions(W)

a1 = nu*inner(grad(u), grad(v))*dx - inner(p, div(v))*dx - inner(q, div(u))*dx
l1 = inner(Constant((0,0)), v)*dx

up = Function(W)
u, p = split(up)
solve(a1 == l1, up, bcu)
3

Thank you for your help dokken !
As I said, I am a beginner, and it had been very difficult to me to fin the exact syntax, especially for function’s space; now, I understand better !

Thank you again

4

You should check out all the demos and search the forum for other posts as there are plenty of examples here

5

Yes, I will take a look ! Thanks

7

There are multiple issues with your code, first of all, I believe you have inserted a really weird inlet condition:

First of all, it is tangential to the inlet, which is strange. Additionally, the inlet velocity is extremely large in comparison to the length scale of the rest of the problem. I also think expecting a stationary solution to this problem is flawed, as the geometry is the classical geometry used in the Turek benchmarks for unsteady flow.

Below I have changed your boundary condition to something a bit more reasonable (with respect to direction and magnitude, as well as writing it in such a way that it can take initial guesses.

from dolfin import *
from mshr import *

set_log_level(LogLevel.PROGRESS)

# Load mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
circle  = Circle(Point(0.2, 0.2), 0.05)
domain  = channel - circle
mesh    = generate_mesh(domain, 200)
tol        = 1E-8
mu         = 1.8E-5
rho        = 1.3
nu         = mu/rho

# Build function space
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
V = FunctionSpace(mesh, P2)
Q = FunctionSpace(mesh, P1)
W  = FunctionSpace(mesh, TH)

# Define boundaries
class walls(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[1] < tol or x[1] > 0.41-tol

class circle(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[1]>0.1 and x[1]<0.3 and x[0]>0.1 and x[0]<0.3

class inflow(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < tol and on_boundary

class outflow(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > 2.2 - tol and on_boundary

# Define boundary conditions

# Corrected boundary condition
bcu_inflow    = DirichletBC(W.sub(0), Constant((1, 0)), inflow())

bcu_walls     = DirichletBC(W.sub(0), Constant((0, 0)), walls())
bcu_circle    = DirichletBC(W.sub(0), Constant((0, 0)), circle())
bcp_outflow = DirichletBC(W.sub(1), Constant(0), outflow())
bcs = [bcu_inflow,bcu_walls, bcu_circle, bcp_outflow]

# Define variational problem
v, q   = TestFunctions(W)

def NS(nu, up_=None):
    if up_ is None:
        up_    = Function(W)

    u_, p_ = split(up_)
    f      = Constant((0, 0)) # Possible forcing term
    nu     = Constant(nu)
    F      = inner(dot(grad(u_), u_), v)*dx + nu*inner(grad(u_), grad(v))*dx \
             - inner(p_, div(v))*dx - inner(q, div(u_))*dx + inner(f, v)*dx

    solve(F == 0, up_, bcs, solver_parameters = {"newton_solver":{ "linear_solver" : "mumps"}})

    return up_

# Does not work
# nu = 5e-4
# up_ = NS(nu)

# Works
nu =1e-3
up_ = NS(nu)
nu = 5e-4
up_ = NS(nu, up_=up_)


u, p = up_.split(deepcopy=True)
XDMFFile("u.xdmf").write_checkpoint(u,"u", 0.0, append=False)

Note that it still doesn’t run for your particular mu, most likely as a combination of this not being a steady state problem, and an issue with mesh refinement and velocity magnitude.

9

Use function evaluation:

x = np.array([0.1,0.2,0]) 
print(u(x))
1 Like