### Problem in Krylov solver for splitting scheme

148
views
0
4 months ago by
Dear FEniCS community. Currently I am working on time dependent convective heat transfer in lid-driven cavity using consistent splitting scheme. However, it shows as follows;

 *** Warning: Krylov solver did not converge in 10000 iterations (PETSc reason DIVERGED_ITS, residual norm ||r|| = 1.914299e+07).
​

Here, I attach my coding for your reference.
from __future__ import print_function
from fenics import *

#define parameter
Tf = 1
Re = 100
Pr = 0.71
Gr = 10
dt = 0.1
num_steps = 10

#create mesh
mesh = UnitSquareMesh(129,129)

#define function space
V = VectorFunctionSpace(mesh,'CG', 2)
Q = FunctionSpace(mesh,'CG', 1)
L = FunctionSpace(mesh,'CG', 2)

#define boundary condition

#lid
bcu_lid = DirichletBC(V, Constant((0,1)), "x[0] < DOLFIN_EPS")
bcT_lid = DirichletBC(L, Constant(1), "x[0] < DOLFIN_EPS")

#right
bcu_right = DirichletBC(V, Constant((0,0)), "x[0] > 1.0 - DOLFIN_EPS")
bcT_right = DirichletBC(L, Constant(0), "x[0] > 1.0 - DOLFIN_EPS")

#top
bcu_top = DirichletBC(V, Constant((0,0)), "x[1] > 1.0 - DOLFIN_EPS")

#bottom
bcu_bottom = DirichletBC(V, Constant((0,0)), "x[1] < DOLFIN_EPS")

bcu = [bcu_lid, bcu_right, bcu_top, bcu_bottom]

bcT = [bcT_lid, bcT_right]

#define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
T = TrialFunction(L)
s = TestFunction(L)

# Create functions
u_ = Function(V)
u_n = Function(V)
p_ = Function(Q)
p_n = Function(Q)
T_ = Function(L)
T_n = Function(L)

#define expression used in variational forms
Re = Constant (Re)
Gr = Constant (Gr)
M = Constant (M)
Pr = Constant (Pr)
z = Gr*pow(Re,-2)
j = Constant((0, z))
k = Constant (dt)

# Step 1: Tentative velocity step
F1 = (1/k)*dot(u - u_n, v)*dx \
- T_n*dot(j,v)*dx

a1 = lhs(F1)
L1 = rhs(F1)

# Step 2:Pressure correction

# Step 3:Pressure update
a3 = p*q*dx
L3 = dot(p_,q)*dx + dot(p_n,q)*dx \
- pow(Re,-1)*div(u_)*q*dx

# Step 4: Temperature correction
F4 = (1/k)*inner(T - T_n, s)*dx \

a4 = lhs(F4)
L4 = rhs(F4)

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
A4 = assemble(a4)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
#[bc.apply(A2) for bc in bcp]
#[bc.apply(A3) for bc in bcp]
[bc.apply(A4) for bc in bcT]

# Create vtk files for visualization output

# Create progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

# Update current time
t += dt
print('current t is',t)

# Step 1
begin("Computing tentative velocity")
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1,'gmres', 'ilu')
end()

# Step 2
begin("Computing pressure correction")
b2 = assemble(L2)
#[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2,'gmres', 'hypre_amg')
end()

# Step 3
begin("Computing pressure update")
b3 = assemble(L3)
#[bc.apply(b3) for bc in bcp]
solve(A3, p_.vector(), b3,'gmres', 'ilu')
end()

# Step 4
begin("Computing temperature correction")
b4 = assemble(L4)
[bc.apply(b4) for bc in bcT]
solve(A4, T_.vector(), b4,'gmres', 'ilu')
end()

# Plot solution
plot(u_, title='Velocity', rescale=True)
plot(p_, title='Pressure', rescale=True)
plot(T_, title='Temperature', rescale=True)

# Save solution to file (vtk)
vtkfile_u << (u_, t)
vtkfile_p << (p_, t)
vtkfile_T << (T_, t)

# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
T_n.assign(T_)

# Update progress bar
progress.update(t / Tf)
print('u max:', u_.vector().array().max())

# Hold plot
interactive()​

I really need your help. Thank you.

Community: FEniCS Project