### Nonlinear variational solver isn't working

75

views

-1

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 \
+ dot(dot(u0,nabla_grad(u)),v)*dx \
+ pow(Re,-1)*inner(nabla_grad(u),nabla_grad(v))*dx \
- Gr*pow(Re,-2)*T0*dot(j,v)*dx \
F3 = (1/k)*inner(T - T0, s)*dx \
- pow(Re,-1)*pow(Pr,-1)*dot(grad(T),grad(s))*dx \
+ inner(inner(grad(T),u1), s)*dx
F = F1 + F2 +F3
# Create vtk files for visualization output
vtkfile_u = File('unsteady_coupled/u.pvd')
vtkfile_p = File('unsteady_coupled/p.pvd')
vtkfile_T = File('unsteady_coupled/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)
#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

### 1 Answer

0

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

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

Please login to add an answer/comment or follow this question.