Nonlinear variational solver isn't working

75
views
-1
13 days ago by
Dear FEniCS community. I need help. Could you please check whats the problem on my nonlinear variational solver? I don't know how to fix it. Thank you.

from __future__ import print_function
from fenics import *

#define dimensionless parameter
Re = 1
Gr = 1
Pr = 0.71
Tf = 10 #final time
num_steps = 100 #number of time steps
dt = 0.1

#create mesh
mesh = RectangleMesh(Point(0,0),Point(20,10),70,70)

#define function space
V = VectorElement('P', mesh.ufl_cell() , 2)
Q = FiniteElement('P', mesh.ufl_cell(), 1)
L = FiniteElement('P', mesh.ufl_cell(), 2)
element = MixedElement([V, Q, L])
W = FunctionSpace(mesh, element)

#define boundaries
walls = 'near(x[1],0)'
upper = 'near(x[1],1)'
inflow = 'near(x[0],0)'
outflow = 'near(x[0],1)'

#define boundary condition

#walls
bcu_walls = DirichletBC(W.sub(0), Constant((0,0)), walls)
bcT_walls = DirichletBC(W.sub(2), Constant(1), walls)

#inflow
bcu_inflow = DirichletBC(W.sub(0), Constant((0,0)), inflow)
bcT_inflow = DirichletBC(W.sub(2), Constant(0), inflow)

#upper
bcu_upper = DirichletBC(W.sub(0), Constant((0,0)), upper)
bcT_upper = DirichletBC(W.sub(2), Constant(0), upper)

bcu = [bcu_walls, bcu_inflow, bcu_upper]
bcT = [bcT_walls, bcT_inflow, bcu_upper]

bcs = bcu + bcT

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

# Create functions
u1 = Function(W)
u0 = Function(W)
p1 = Function(W)
p0 = Function(W)
T1 = Function(W)
T0 = Function(W)

#split functions
w = Function(W)
(u, p, T) = split(w)

w0 = Function(W)
(u0, p0, T0) = split(w0)

w1 = Function(W)
(u1, p1, T1) = split(w1)

y = TestFunction(W)
(v, q, s) = split(y)

#define expression used in variational forms
Re = Constant (Re)
Gr = Constant (Gr)
Pr = Constant (Pr)
k = Constant (dt)
j = Constant((0, 1))

#equation
F1 = div(u)*q*dx

F2 = (1/k)*dot(u - u0, v)*dx \
- Gr*pow(Re,-2)*T0*dot(j,v)*dx \

F3 = (1/k)*inner(T - T0, s)*dx \

F = F1 + F2 +F3

# 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)

#solve nonlinear variational problem
J = derivative(F, w)
problem = NonlinearVariationalProblem(F, w, bcs, J)
solver  = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm["newton_solver"]["absolute_tolerance"] = 1E-8
prm["newton_solver"]["relative_tolerance"] = 1E-7
prm["newton_solver"]["maximum_iterations"] = 1000
prm["newton_solver"]["relaxation_parameter"] = 1.0
prm["newton_solver"]["linear_solver"] = "gmres"
prm["newton_solver"]["preconditioner"] = "ilu"
prm["newton_solver"]["krylov_solver"]["absolute_tolerance"] = 1E-9
prm["newton_solver"]["krylov_solver"]["relative_tolerance"] = 1E-7
prm["newton_solver"]["krylov_solver"]["maximum_iterations"] = 1000
solver.solve()

# Plot solution
plot(u1, title='Velocity')
plot(p1, title='Pressure')
plot(T1, title='Temperature')

# Save solution to file (vtk)
(u, p, T) = w.split()
vtkfile_u << (u1, t)
vtkfile_p << (p1, t)
vtkfile_T << (T1, t)

# Update previous solution
w0.assign(w1)

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

# Hold plot
interactive()​
Community: FEniCS Project
This is clearly not a minimal example reproducing the error.
written 13 days ago by klunkean

0
13 days ago by
There are many issues in your code:

--Use of derivative for the Jacobian in the solve is wrong.
--For the split u1, p1, T1 are used, but in the weak form they are used as u, p, T.

Best, Emek