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 \
+ dot(dot(u_n,nabla_grad(u)),v)*dx \
+ pow(Re,-1)*inner(nabla_grad(u),nabla_grad(v))*dx \
+ dot(grad(p_n),v)*dx \
- T_n*dot(j,v)*dx

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

# Step 2:Pressure correction
a2 = inner(grad(p), grad(q))*dx
L2 = (1/k)*inner(u_- u_n, grad(q))*dx

# 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 \
- pow(Re,-1)*pow(Pr,-1)*dot(grad(T),grad(s))*dx \
+ inner(inner(grad(T),u_), 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
vtkfile_u = File('unsteady_free_lidD/u.pvd')
vtkfile_p = File('unsteady_free_lidD/p.pvd')
vtkfile_T = File('unsteady_free_lidD/T.pvd')

# 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
Please login to add an answer/comment or follow this question.

Similar posts:
Search »