Convection Flow in a Channel Problem

Errors
1

Hi. I need a favor in pointing the problem in my coding. I don’t know what’s missing in my source code, but for some reason, the solver kept on diverging.

Now I am working on convective flow in a channel with an open trapezoidal enclosure. I’m working for weeks but could not run it successfully.

I have tried to solve it using Newton Solver and SNES Solver, but it seems both produces the same error. The error appear:

*** Error: Unable to solve nonlinear system with PETScSNESSolver.
*** Reason: Solver did not converge.
*** Where: This error was encountered inside PETScSNESSolver.cpp.
*** Process: 0

Please help me. Here I attach the code:

# Import libraries
from __future__ import print_function
from fenics import *
from mshr import *
import matplotlib.pyplot as plt

# Create the main channel
channel = Rectangle(
                    Point(0., 1), # Bottom-Left
                    Point(5, 1.5) # Top-Right
                   )

# Create the trapezoidal polygon
trapezoid = Polygon([
                     Point(1.5, 0.0),  # Bottom-left
                     Point(2.5, 0.0),  # Bottom-right
                     Point(3, 1),  # Top-right
                     Point(1, 1)   # Top-left
                    ])

domain = channel + trapezoid

mesh = generate_mesh (domain, 128)
#mesh = generate_mesh (domain, 84.16)

#visualize the mesh
plot(mesh)

# Define vector element and finite element according to their order
P1 = VectorElement("CG", mesh.ufl_cell(), 2)
P2 = FiniteElement("CG", mesh.ufl_cell(), 1)
P3 = FiniteElement("CG", mesh.ufl_cell(), 2)

# Define function space
W = FunctionSpace(mesh, MixedElement([P1,P2,P3]))

# Define test function and trial function
(v,q,s) = TestFunctions(W)

w = Function(W)
(u,p,T) = (as_vector ((w[0], w[1])), w[2], w[3])

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1.5)
        
class ChannelExit(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 5.)

class ChannelEntrance(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.)

class BottomChannel(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1.)

class SideWallCavity(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[1] >= 0. and x[1] <= 1. and x[0] <= 3. and x[0] >= 1.
        
class BottomCavity(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.)

class Heater(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.) and x[0] >= 1.75 and x[0] <= 2.25

boundaries = MeshFunction ("size_t", mesh, 1)

boundaries.set_all(0)

chexit = ChannelExit()
chexit.mark(boundaries, 1)

chtop = Top()
chtop.mark(boundaries, 2)

chent = ChannelEntrance()
chent.mark(boundaries, 3)

chbot = BottomChannel()
chbot.mark(boundaries, 4)

cavsw = SideWallCavity()
cavsw.mark(boundaries, 5)

cavbot = BottomCavity()
cavbot.mark(boundaries, 6)

heater = Heater()
heater.mark(boundaries, 7)

#Velocity
v_chent = DirichletBC (W.sub(0), (1,0), boundaries, 3)

v_chtop = DirichletBC (W.sub(0), (0,0), boundaries, 2)
v_chbot = DirichletBC (W.sub(0), (0,0), boundaries, 4)
v_cavsw = DirichletBC (W.sub(0), (0,0), boundaries, 5)
v_cavbot = DirichletBC (W.sub(0), (0,0), boundaries, 6)
v_heater = DirichletBC (W.sub(0), (0,0), boundaries, 7)

#Temperature
t_chent = DirichletBC(W.sub(2), 0, boundaries, 3)
t_heater = DirichletBC(W.sub(2), 1, boundaries, 7)

#Momentum
p_chexit = DirichletBC(W.sub(1), 0, boundaries, 1)

bc = [v_chtop, v_chent, v_chbot, v_cavsw, v_cavbot, v_heater, t_chent, t_heater, p_chexit]

Re = 100.0
Pr = 1.3
Ri = 5

invRe = 1.0/Re
invPr = 1.0/Pr

Ri = Constant(Ri)
e = Constant((0,1))

# Coupled Equations
F1 = div(u)*q*dx

F2 = inner(grad(u)*u, v)*dx\
    + invRe*inner(grad(u),grad(v))*dx\
    - p*div(v)*dx\
    - T*Ri*dot(e,v)*dx

F3 = dot(u, grad(T))*s*dx\
    + invRe*invPr*dot(grad(T), grad(s))*dx

F = F1 + F2 + F3

# Run solver
J = derivative(F, w)

#problem = NonlinearVariationalProblem(F, w, bc, J)
#solver = NonlinearVariationalSolver(problem)

# Non linear solver parameters
#prm = solver.parameters

#prm['newton_solver']['linear_solver'] = 'mumps'
#prm['newton_solver']['maximum_iterations'] = 500

# Solve the problem
#solver.solve()

# Create the nonlinear problem and solver
problem = NonlinearVariationalProblem(F, w, bc, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters

# Set SNES solver parameters

prm['nonlinear_solver'] = 'snes'
prm['snes_solver']['method'] = 'newtonls'  # Choose SNES method
prm['snes_solver']['linear_solver'] = 'mumps'  # Linear solver for the Jacobian

# Solve the problem
solver.solve()