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 \
+ 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
13 days ago by
Emek  
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
Please login to add an answer/comment or follow this question.

Similar posts:
Search »