Problem with solver

1

Hi! I am using FEniCS version 2017.2.0 and I try solver the lid-driven cavity problem with transient Navier-Stokes equation with stabilization SUPG/PSPG/LSIC.

My code:
from future import print_function
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import math

------------- Problem data (cavity) -----------------------

N = 16 # mesh division number
f = Constant((0.0,0.0)) # body force
Re = 1000.0 # Reynolds number
nu = 1.0/Re # kinematic viscosity
Vmax = 1.0 # maximum velocity
T = 2.5 # final time

def inicial_condition():
v0 = Constant((0.0, 0.0))
return v0

class Lid(SubDomain):
def inside(self, x, on_boundary):
return(on_boundary and near (x[1], 1.0))

class Walls(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[0], 0.0) or near(x[0], 1.0) or near(x[1], 0.0))

def zero(x):
return (near(x[0], 0.0) and near(x[1], 0.0))

def boundary_condition(W):
# Dirichlet value
g1 = Constant ((1.0, 0.0))
g2 = Constant ((0.0, 0.0))
g3 = Constant (0.0)
# conditions
bc1 = DirichletBC(W.sub(0), g1, Lid())
bc2 = DirichletBC(W.sub(0), g2, Walls())
bc3 = DirichletBC(W.sub(1), g3, zero, ‘pointwise’)
bcs = [bc1, bc2, bc3]
return bcs

--------------- code -------------------

mesh = UnitSquareMesh(N, N)
h = 1.0/N

time step

dt = 0.2*(h/Vmax)
n_aux = T/dt + 1.0
n = int(n_aux)
dt = T/n
t_range = np.linspace(0,T,n+1) [1:]

define function space

P1v = VectorElement(“Lagrange”, mesh.ufl_cell(), 1)
P1 = FiniteElement(“Lagrange”, mesh.ufl_cell(), 1)
TH = P1v*P1
W = FunctionSpace(mesh, TH)

trial and test function

vp = TrialFunction(W)
(v, p) = split(vp)
(w, q) = TestFunctions(W)
vp_ = Function(W)
#v_, p_ = split(vp_)

set conditions

v0 = inicial_condition()
vn = interpolate(v0, W.sub(0).collapse())
bcs = boundary_condition(W)

Galerkin formulation

F = (inner(v - vn, w)/dt + inner(dot(vn,nabla_grad(v)), w) + nu*inner(grad(w), grad(v)) - inner(p,div(w)) - inner(q,div(v)) - inner(f, w))*dx

Stabilization parameters

vnorm = sqrt(dot(v, v))

momentum residue

R = (1.0/dt)(v-vn) + dot(v,nabla_grad(v)) - nudiv(grad(v)) + grad§ - f

PSPG

tau_pspg = ((2.0/dt)2 + (2.0vnorm/h)*2 + 9.0(4.0nu/(h2))2)(-0.5)#1.0/sqrt( pow(2.0/dt,2) + pow(2.0vnorm)/h,2) + 9.0pow((4.0nu)/(h**2),2))
F_pspg = tau_pspginner(R, grad(q))*dx

SUPG

tau_supg = ((2.0/dt)2 + (2.0vnorm/h)*2 + 9.0(4.0nu/(h2))2)(-0.5)
F_supg = tau_supg*inner(R, dot(v, grad(w)))*dx

LSIC

tau_lsic = (vnormh)/2.0
F_lsic = tau_lsicinner(div(v), div(w))*dx

Add stabilization terms

F += -F_pspg + F_supg + F_lsic

Jacobian

F1 = action(F, vp_)
J = derivative(F1, vp_, vp)

time loop

for t in t_range:
print(’\nt = {}\n’.format(t))
problem = NonlinearVariationalProblem(F1, vp_, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm[“nonlinear_solver”]=“newton”
prm[“newton_solver”][“absolute_tolerance”] = 1E-4
prm[“newton_solver”][“relative_tolerance”] = 1E-4
prm[“newton_solver”][“maximum_iterations”] = 20
prm[“newton_solver”][“relaxation_parameter”] = 1.0
prm[“newton_solver”][“linear_solver”] = “gmres”
prm[“newton_solver”][“krylov_solver”][“absolute_tolerance”] = 1E-8
prm[“newton_solver”][“krylov_solver”][“relative_tolerance”] = 1E-8
prm[“newton_solver”][“krylov_solver”][“maximum_iterations”] = 1000
prm[“newton_solver”][“krylov_solver”][“monitor_convergence”] = True
prm[“newton_solver”][“krylov_solver”][“nonzero_initial_guess”] = False
#prm[“newton_solver”][“krylov_solver”][“gmres”][“restart”] = 40
prm[“newton_solver”][“preconditioner”] = “ilu”
#prm[“newton_solver”][“krylov_solver”][“preconditioner”][“structure”] = “same_nonzero_pattern”
#prm[“newton_solver”][“krylov_solver”][“preconditioner”][“ilu”][“fill_level”] =0

PROGRESS = 16
set_log_level(PROGRESS)
solver.solve()

v_ , p_ = vp_.split(True)
print(v_.vector().array())

t = t + dt
vn.assign(v_)

plt.figure()
plot(v_)
plt.figure()
plot(p_)
plt.show()

The parameters commented in solver don’t work.
The error:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = -nan (tol = 1.000e-04) r (rel) = -nan (tol = 1.000e-04)
PETSc Krylov solver starting to solve 867 x 867 system.
Traceback (most recent call last):
File “cavidade_SUPG.py”, line 124, in
solver.solve()
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** fenics-support@googlegroups.com

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0

*** DOLFIN version: 2017.2.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

I don’t know what is wrong, because appear -nan in first iteration of Newton.

1 Like
2

Check if it is working with a direct solver. It’s most probably a problem with your formulation and not the solver.

3

Since it failed to converge in 0 iterations it is probably a problem with your preconditioner. For solving SUPG problems in 3D I use the following solver parameters:

# Create the Newton solver
problem = NonlinearVariationalProblem(F_P, cp, bccp, J_P)
solver = NonlinearVariationalSolver(problem)
solverType = 'newton'
#solverType = 'snes' # PETSc SNES has some more options for line search etc.
prm = solver.parameters
solver.parameters['nonlinear_solver'] = solverType
nparam = solver.parameters[solverType+'_solver']
nparam['linear_solver'] = 'gmres'
nparam['preconditioner'] = 'jacobi' # or 'none'

Try chaningin the linear solver and preconditioner to ‘gmres’ and ‘jacobi’. Setting the preconditioner to ‘none’ also works for me, but sometimes gives me longer simulation times.

4

if you have

and

R = (1.0/dt) *(v-vn)

will it not yield a 1/0, thus nan at the first iteration?