Sinusoidal cavity problem

1

Hi. I need a favor in pointing the problem in my coding. I don’t know whats wrong with my boundary conditions. Now I am working on convective flow in a sinusoidal cavity. I’m working for weeks but could not run it successfully. The error appear:

*** Warning: Found no facets matching domain for boundary condition.

Please help me. Here I attach the code:

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

# Define dimensionless parameter
Ra = 10000
Pr = 6.26
n = 1
M = 0
Nr = 0.1
Nb = 0.1
Nt = 0.1
Le = 10
o = 0
k = 0

# Create mesh
domain_n_points = 128
domain_points = list()
for n in range(domain_n_points + 1):
    y = n*1./domain_n_points
    domain_points.append( Point(0.2-0.2*cos(6.*pi*y), y) )
domain_points.append( Point(1., 1.) )
domain_points.append( Point(1., 0.) )
domain_points = domain_points[::-1] # Polygon vertices must be given in counter clockwise order.
domain = Polygon(domain_points)
mesh = generate_mesh(domain, 16)
plot(mesh)

# Create boundaries
class Noslip(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1]) < DOLFIN_EPS

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 1) < DOLFIN_EPS

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary #and abs(x[0]) < DOLFIN_EPS

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1] - 1) < DOLFIN_EPS

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
noslip = Noslip() 
noslip.mark(boundaries,1)
top = Top() 
top.mark(boundaries,2)
bottom = Bottom()
bottom.mark(boundaries, 3)
right = Right()
right.mark(boundaries, 4)
left = Left()
left.mark(boundaries, 5)
plot(boundaries)

# Define function space
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
L = FiniteElement("CG", mesh.ufl_cell(), 2)
X = FiniteElement("CG", mesh.ufl_cell(), 2)

element = MixedElement([V, Q, L, X])
W = FunctionSpace(mesh, element)

# Define boundary condition
no_slip  = DirichletBC(W.sub(0), (0,0), boundaries,1)
#lid  = DirichletBC(W.sub(0), (1,0), "x[1] > 1.0 - DOLFIN_EPS")
cold = DirichletBC(W.sub(2), 0,boundaries, 4)
hot = DirichletBC(W.sub(2), 1,boundaries, 5)
notconc = DirichletBC(W.sub(3), 0,boundaries, 4)
conc = DirichletBC(W.sub(3), 1,boundaries, 5)

bc = [no_slip, hot, cold, notconc, conc]


# Define test functions
(v,q,s,r) = TestFunctions(W)

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

# Define expression used in variational forms
miu = 1.0/Pr
theta = 1.0/Le
e = Constant ((0,1))
M = Constant(M) 
Nr = Constant(Nr)
Nb = Constant(Nb)
Nt = Constant(Nt)
phi = Nt/Nb
o = Constant(o)
k = Constant(k)

def epsilon(w):
	return (grad(w)+(grad(w).T))

def secondinvariant(w):
	return 0.5*inner(epsilon(w),epsilon(w))+ DOLFIN_EPS

def powerlaw(w):
	return pow(secondinvariant(w),0.5*(n-1))

# Define variational problem
F1 = div(u)*q*dx 

F2 = inner(dot(u,grad(u)), v)*dx - div(v)*p*dx + Pr*powerlaw(u)*inner(epsilon(u),grad(v))*dx \
+ M*dot(u,v)*dx - T*Pr*Ra*dot(e,v)*dx + c*Nr*Pr*Ra*dot(e,v)*dx

F3 = dot(u,grad(T))*s*dx + dot(grad(T),grad(s))*dx \
- Nb*inner(grad(c),grad(T))*s*dx - Nt*inner(grad(T),grad(T))*s*dx- o*T*s*dx

F4 = dot(u,grad(c))*r*dx + theta*dot(grad(c),grad(r))*dx \
+ theta*phi*dot(grad(T),grad(r))*dx - k*c*r*dx

F = F1 + F2 + F3 + F4

# Solve 
J = derivative(F, w)
problem = NonlinearVariationalProblem(F, w, bc, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve

# Create VTK files for visualization output
vtkfile_u = File('test/u.pvd')
vtkfile_p = File('test//p.pvd')
vtkfile_T = File('test//T.pvd')
vtkfile_c = File('test//c.pvd')

# Save solution to file (VTK)
(u, p, T, c) = w.split()
vtkfile_u << u
vtkfile_p << p
vtkfile_T << T
vtkfile_c << c
	
# Plot solution
plot(u)
plot(T)
plot(c)

# Hold plot
plt.show()
2

Hi,
your Left subdomain tags all the boundary, consider changing it to

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] < 0.5 and x[1] > DOLFIN_EPS and x[1] < 1-DOLFIN_EPS

to get only the sinusoidal left part. You can then check that boundaries.array() contains all index from 0 to 5.

3

Thank you bleyerj for your reply. Can I ask why x[0] < 0.5?

4

Hi, since your left boundary is curved it is not possible to have an exact equation such as x[0]==0 for instance. Hence, the condition gets all the facet which are less than x=0.5 (your sinusoid maximum being less than 0.5), which are on the boundary and which are not on the bottom and top boundary. Then you get only those lying on the boundary of the sinusoid.

5

I understood. Thank you so much for your help. The problem of b.c settled, however, unable to solve linear system. Before this, I have run this problem successfully using square cavity. However, it failed using sinusoid. Hope u can give me some ideas on how to solve it. Thank you.

6

Could you supply a minimal code that is failing?

7

This is the error:

Error:   Unable to solve nonlinear system with NewtonSolver.
Reason:  Newton solver did not converge because maximum number of iterations reached.
Where:   This error was encountered inside NewtonSolver.cpp.
Process: 0
8

When Im using linear solver = gmres, and preconditioner = ilu, there come the error:

Error:   Unable to solve linear system using PETSc Krylov solver.
Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PCSETUP_FAILED, residual norm ||r|| = 0.000000e+00).
Where:   This error was encountered inside PETScKrylovSolver.cpp.
Process: 0
9

Full code is attached above. The iteration start with nan.

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
10

I would recommend trying snes instead of Newton, but it looks like that there could be a problem with the nonlinear residual itself.

Are you sure that the expression for F is computed correctly? Also, for marking the mesh properly for assigning dirichlet boundary conditions, I would use Gmsh instead of mshr

11

Yes, the expression of F is already tested and solved successfully using square cavity. I never try to use Gmsh, will try to study it.