# Issue with sudden 'nan' error

Asked by Marcos Samudio on 2012-03-12

Hello everyone

As I stated in my previous question, I am having issues when running transient simulations. This is a basic code that shows the problem I go through when trying to simulate unsteady NS equations, with PSPG stabilization and no RANS model. It must be placed in the upper directory for import statements to work properly.

import cbc.cfd.icns as icns
from cbc.cfd.icns import solver_parameters
from NSProblem import *

class LidDrivenCavity(NSProblem):

def __init__(self, parameters, mesh_size):
parameters['problem_name'] = self.__class__.__name__
NSProblem.__init__(self, parameters=parameters)

self.mesh = mesh = UnitSquare(mesh_size, mesh_size)
self.mf = FacetFunction("uint", mesh)
self.mf.set_all(0)

#Set some constants
self.nu = self.prm['viscosity'] = 1./self.prm['Re']
self.U = 1.0 #Max. velocity

#Determine timestep size
if self.prm['time_integration'] == 'Transient':
self.prm['dt'] = self.timestep()

#Initialize u and p
self.q0 = Initdict({'p': '0', 'u': ('0', '0')})

#Create boundary conditions
self.boundaries = self.create_boundaries()

#Define lid boundary function
def lid_func(self):
if self.prm['time_integration'] == 'Transient':
#We set a large t. This will be changed if we choose to use smoothly
#growing BC. Additional code has been written for this purpose and
# and is not included here
func = {'u': Expression(("1-exp(-0.5*t)", "0."), t=1e10)}
else:
func = {'u': Constant((1., 0.))}
return func

def create_boundaries(self):

inlet = FlowSubDomain(lambda x, on_boundary: near(x, 1.) and on_boundary,
bc_type='Wall',
func = self.lid_func(),
mf = self.mf)

walls = FlowSubDomain(lambda x, on_boundary: on_boundary and (near(x,0.0)\
or near(x,1.0) or near(x,0.0)) ,
bc_type = 'Wall',
mf = self.mf)

return [inlet, walls]

def __info__(self):
return "Lid Driven Cavity"

#############################################
####problem_parameters <== NSProblem.py <== Problem.py####
#############################################

problem_parameters['cfl'] = 0.5
problem_parameters['T'] = 50.
problem_parameters['dt'] = 0.2
problem_parameters['Re'] = 500.
problem_parameters['viscosity'] = None
problem_parameters['Re_tau'] = Re_tau = 395.0
problem_parameters['utau'] = utau = 0.05
problem_parameters['turbulence_model'] = 'Laminar'
problem_parameters['save_solution'] = False
problem_parameters['plot_velocity'] = True
problem_parameters['plot_pressure'] = False
problem_parameters['max_iter'] = 300
problem_parameters['iter_first_timestep'] = 25 #To improve accuracy on first step
problem_parameters['time_integration'] = 'Transient'
problem_parameters['smooth_transient_bcs'] = False
problem_parameters['max_err'] = 1e-07
problem_parameters['save_steps'] = 1 # Timesteps between saves

############################################
####solver_parameters <== NSSolver.py <== PDESystem.py####
############################################

#DICTIONARY UPDATE
solver_parameters = recursive_update(solver_parameters,
dict(degree=dict(u=1, u0=1, u1=1, p=1),
pdesubsystem=dict(u=101, p=101, velocity_update=101, up=1),
linear_solver=dict(u='bicgstab', p='gmres', velocity_update='bicgstab', up='lu'),
precond=dict(u='jacobi', p='amg', velocity_update='ilu')))
#SOLVER PARAMETERS SETUP
solver_parameters['degree']['u']=1
solver_parameters['degree']['u0']=1
solver_parameters['degree']['u1']=1
solver_parameters['degree']['p']=1
solver_parameters['pdesubsystem']['u']=101
solver_parameters['pdesubsystem']['p']=101
solver_parameters['pdesubsystem']['velocity_update']=101
solver_parameters['pdesubsystem']['up']=1
solver_parameters['linear_solver']['u']='bicgstab'
# solver_parameters['linear_solver']['u']='lu'
solver_parameters['linear_solver']['p']='gmres'
# solver_parameters['linear_solver']['p']='lu'
solver_parameters['linear_solver']['velocity_update']='bicgstab'
# solver_parameters['linear_solver']['velocity_update']='lu'
solver_parameters['linear_solver']['up']='lu'
solver_parameters['precond']['u']='jacobi'
# solver_parameters['precond']['u']='none'
solver_parameters['precond']['p']='amg'
# solver_parameters['precond']['p']='none'
# solver_parameters['precond']['velocity_update']='ilu'
solver_parameters['precond']['velocity_update']='jacobi'
solver_parameters['omega'].default_factory = lambda : 0.8
solver_parameters['convection_form'] = 'Standard'
solver_parameters['max_iter'] = 1
solver_parameters['familyname'] = 'Navier-Stokes'
solver_parameters['iteration_type'] = 'Picard'
solver_parameters['max_err'] = 1e-07

###############################
# PROBLEM FORMULATION #
###############################
problem= LidDrivenCavity(problem_parameters, 20)
NS_solver = icns.NSCoupled(problem, solver_parameters)
problem.solve(logging = True)

Running this code yields a 'nan' error at time 0.2749, without showing any signs of divergence. It just keeps happening with different values of Re and different meshes.

I really think this is a bug, since changing slightly some parameters (like taking Re to be 500.00001 instead of 500., and other similar changes) makes the simulation run longer.

I hope someone can find what is going wrong here. If you have any clues or suggestions, please let me know.

Marcos

## Question information

Language:
English Edit question
Status:
Solved
For:
CBC.PDESys Edit question
Assignee:
No assignee Edit question
Solved by:
Marcos Samudio
Solved:
Last query:
 Mikael Mortensen (mikael-mortensen) said on 2012-03-12: #1

Hi Marcos,

On 12 March 2012 16:25, Marcos Samudio <<email address hidden>
> wrote:

> New question #190440 on CBC.PDESys:
>
> Hello everyone
>
> As I stated in my previous question, I am having issues when running
> transient simulations. This is a basic code that shows the problem I go
> through when trying to simulate unsteady NS equations, with PSPG
> stabilization and no RANS model. It must be placed in the upper directory
> for import statements to work properly.
>
> import cbc.cfd.icns as icns
> from cbc.cfd.icns import solver_parameters
> from NSProblem import *
>
> class LidDrivenCavity(NSProblem):
>
> def __init__(self, parameters, mesh_size):
> parameters['problem_name'] = self.__class__.__name__
> NSProblem.__init__(self, parameters=parameters)
>
> self.mesh = mesh = UnitSquare(mesh_size, mesh_size)
> self.mf = FacetFunction("uint", mesh)
> self.mf.set_all(0)
>
> #Set some constants
> self.nu = self.prm['viscosity'] = 1./self.prm['Re']
> self.U = 1.0 #Max. velocity
>
> #Determine timestep size
> if self.prm['time_integration'] == 'Transient':
> self.prm['dt'] = self.timestep()
>
> #Initialize u and p
> self.q0 = Initdict({'p': '0', 'u': ('0', '0')})
>
> #Create boundary conditions
> self.boundaries = self.create_boundaries()
>
> #Define lid boundary function
> def lid_func(self):
> if self.prm['time_integration'] == 'Transient':
> #We set a large t. This will be changed if we choose to use
> smoothly
> #growing BC. Additional code has been written for this purpose
> and
> # and is not included here
> func = {'u': Expression(("1-exp(-0.5*t)", "0."), t=1e10)}
> else:
> func = {'u': Constant((1., 0.))}
> return func
>
> def create_boundaries(self):
>
> inlet = FlowSubDomain(lambda x, on_boundary: near(x, 1.) and
> on_boundary,
> bc_type='Wall',
> func = self.lid_func(),
> mf = self.mf)
>
> walls = FlowSubDomain(lambda x, on_boundary: on_boundary and
> (near(x,0.0)\
> or near(x,1.0) or near(x,0.0)) ,
> bc_type = 'Wall',
> mf = self.mf)
>
> return [inlet, walls]
>
> def __info__(self):
> return "Lid Driven Cavity"
>

It's much easier to follow if we strip this down to neccessary parameters:

problem_parameters['cfl'] = 0.5
problem_parameters['T'] = 50.
problem_parameters['Re'] = 500.
problem_parameters['plot_velocity'] = True
solver_parameters = recursive_update(solver_parameters,
dict(degree=dict(u=1, p=1),
pdesubsystem=dict(up=1),
stabilization_prm=0.01,
max_iter=10))

# Solve with
problem= LidDrivenCavity(problem_parameters, 20)
NS_solver = icns.NSCoupled(problem, solver_parameters)
problem.solve()

No need for the rest.

You are using a coupled Navier-Stokes solver with the same degree on
pressure and velocity functionspaces. For this reason a term

is added to the regular Navier-Stokes equation and this seems to cause
problems. I stripped the problem down to a much easier to follow script
that still creates the error:

from cbc.pdesys import *

mesh = UnitSquare(20, 20)
V = VectorFunctionSpace(mesh, 'CG', 1)
Q = FunctionSpace(mesh, 'CG', 1)
W = V * Q

u, p = TrialFunctions(W)
v, q = TestFunctions(W)
up_ = Function(W)
up_1 = Function(W)
up_2 = Function(W)
u_, p_ = up_.split()
u_1, p_1 = up_1.split()
u_2, p_2 = up_2.split()

nu = Constant(1./500.)
f = Constant((0., 0.))
dt = Constant(0.035)

bc1 = DirichletBC(V, (0., 0.), "on_boundary")
bc2 = DirichletBC(V, (1., 0.), "on_boundary && x > 1. - DOLFIN_EPS")
bcs = [bc1, bc2]
normalization = extended_normalize(W, 2)
eps = Constant(0.01)

U = 0.5*(u + u_1)
U1 = 1.5*u_1 - 0.5*u_2
F = (1./dt)*inner(u - u_1, v)*dx + inner(v, dot(U1, nabla_grad(U)))*dx \
- inner(div(v), p)*dx - inner(q, div(u))*dx - inner(v, f)*dx \

tmp = Function(W)
a, L = lhs(F), rhs(F)
t = 0.
while t < 2.5:
t += dt(0)
i = 0; error=1.
while i < 10 and error > 1e-7:
i += 1
A = assemble(a)
b = assemble(L)
[bc.apply(A, b) for bc in bcs]
tmp.vector()[:] = up_.vector()[:]
solve(A, up_.vector(), b)
normalization(up_.vector())
tmp.vector().axpy(-1, up_.vector())
error = tmp.vector().norm('l2')
print 'error ', i , error
up_2.vector()[:] = up_1.vector()[:]
up_1.vector()[:] = up_.vector()[:]
print 'Time = ', t

plot(u_, title='Velocity')
plot(p_, title='Pressure')

I get the same 'nan' error about the same time. Not sure why yet, but I
think this is a better place to start debugging. I'll have a closer look
tomorrow.

Best regards

Mikael

>
> #############################################
> ####problem_parameters <== NSProblem.py <== Problem.py####
> #############################################
>
> problem_parameters['cfl'] = 0.5
> problem_parameters['T'] = 50.
> problem_parameters['dt'] = 0.2
> problem_parameters['Re'] = 500.
> problem_parameters['viscosity'] = None
> problem_parameters['Re_tau'] = Re_tau = 395.0
> problem_parameters['utau'] = utau = 0.05
> problem_parameters['turbulence_model'] = 'Laminar'
> problem_parameters['save_solution'] = False
> problem_parameters['plot_velocity'] = True
> problem_parameters['plot_pressure'] = False
> problem_parameters['max_iter'] = 300
> problem_parameters['iter_first_timestep'] = 25 #To improve accuracy on
> first step
> problem_parameters['time_integration'] = 'Transient'
> problem_parameters['smooth_transient_bcs'] = False
> problem_parameters['max_err'] = 1e-07
> problem_parameters['save_steps'] = 1 # Timesteps between saves
>
> ############################################
> ####solver_parameters <== NSSolver.py <== PDESystem.py####
> ############################################
>
> #DICTIONARY UPDATE
> solver_parameters = recursive_update(solver_parameters,
> dict(degree=dict(u=1, u0=1, u1=1, p=1),
> pdesubsystem=dict(u=101, p=101, velocity_update=101, up=1),
> linear_solver=dict(u='bicgstab', p='gmres',
> velocity_update='bicgstab', up='lu'),
> precond=dict(u='jacobi', p='amg', velocity_update='ilu')))
> #SOLVER PARAMETERS SETUP
> solver_parameters['degree']['u']=1
> solver_parameters['degree']['u0']=1
> solver_parameters['degree']['u1']=1
> solver_parameters['degree']['p']=1
> solver_parameters['pdesubsystem']['u']=101
> solver_parameters['pdesubsystem']['p']=101
> solver_parameters['pdesubsystem']['velocity_update']=101
> solver_parameters['pdesubsystem']['up']=1
> solver_parameters['linear_solver']['u']='bicgstab'
> # solver_parameters['linear_solver']['u']='lu'
> solver_parameters['linear_solver']['p']='gmres'
> # solver_parameters['linear_solver']['p']='lu'
> solver_parameters['linear_solver']['velocity_update']='bicgstab'
> # solver_parameters['linear_solver']['velocity_update']='lu'
> solver_parameters['linear_solver']['up']='lu'
> solver_parameters['precond']['u']='jacobi'
> # solver_parameters['precond']['u']='none'
> solver_parameters['precond']['p']='amg'
> # solver_parameters['precond']['p']='none'
> # solver_parameters['precond']['velocity_update']='ilu'
> solver_parameters['precond']['velocity_update']='jacobi'
> solver_parameters['omega'].default_factory = lambda : 0.8
> solver_parameters['convection_form'] = 'Standard'
> solver_parameters['max_iter'] = 1
> solver_parameters['familyname'] = 'Navier-Stokes'
> solver_parameters['iteration_type'] = 'Picard'
> solver_parameters['max_err'] = 1e-07
>
> ###############################
> # PROBLEM FORMULATION #
> ###############################
> problem= LidDrivenCavity(problem_parameters, 20)
> NS_solver = icns.NSCoupled(problem, solver_parameters)
> problem.solve(logging = True)
>
> Running this code yields a 'nan' error at time 0.2749, without showing any
> signs of divergence. It just keeps happening with different values of Re
> and different meshes.
>
> I really think this is a bug, since changing slightly some parameters
> (like taking Re to be 500.00001 instead of 500., and other similar changes)
> makes the simulation run longer.
>
> I hope someone can find what is going wrong here. If you have any clues or
> suggestions, please let me know.
>
> Marcos
>
> --
> contact for CBC.PDESys.
>

 Marcos Samudio (marcossamudio) said on 2012-03-12: #2

A lot easier to debug indeed! I will take a look at it too. I noticed you did not add the boundary terms that are present in the pspg stabilization, is there a reason for that?

Marcos

 Marcos Samudio (marcossamudio) said on 2012-03-13: #3

I have tried the following and got acceptable results:

bc1 = DirichletBC(V, (0., 0.), "on_boundary")
bc2 = DirichletBC(V, (1., 0.), "on_boundary && x > 1. - DOLFIN_EPS")
bc3 = DirichletBC(Q, 0., "on_boundary && x > DOLFIN_EPS && x > DOLFIN_EPS")
bcs = [bc1, bc2, bc3]

This way I avoid the use of extended_normalize. The error is quite large, but the simulation yields good looking velocity and pressure fields.

My guess is that extended_normalize could be causing the this problem, so I will take a further look at it tomorrow.

Marcos

 Marcos Samudio (marcossamudio) said on 2012-03-13: #4

It should be
bc3 = DirichletBC(Q, 0., "on_boundary && x < DOLFIN_EPS && x < DOLFIN_EPS")

Sorry!

 Mikael Mortensen (mikael-mortensen) said on 2012-03-13: #5

The problem is not extended_normalize, but the fact that the pressure is
underresolved as it can contain an arbitrary constant due to only Neumann
boundary conditions. This seems to bother the PETSc direct solver, whereas
the problem runs just fine for me with iterative solvers, try bicgstab with
ilu. I can't make your example with three boundary conditions work, but I
see what you're trying to do. There's another really nice solution to this
problem, by constraining pressure through the introduction of an additional
space of real numbers. Try this code and look at neumann-poisson demo:

from cbc.pdesys import *

mesh = UnitSquare(20, 20)
V = VectorFunctionSpace(mesh, 'CG', 1)
Q = FunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh, 'R', 0)
W = MixedFunctionSpace([V, Q, R])

u, p, c = TrialFunctions(W)
v, q, d = TestFunctions(W)

up_ = Function(W)
up_1 = Function(W)
up_2 = Function(W)
u_, p_, c_ = up_.split()
u_1, p_1, c_1 = up_1.split()
u_2, p_2, c_2 = up_2.split()

nu = Constant(1./500.)
f = Constant((0., 0.))
dt = Constant(0.035)

bc1 = DirichletBC(V, (0., 0.), "on_boundary")
bc2 = DirichletBC(V, (1., 0.), "on_boundary && x > 1. - DOLFIN_EPS")
#bc3 = DirichletBC(Q, 0., "on_boundary && x < DOLFIN_EPS && x <
DOLFIN_EPS", 'pointwise')
bcs = [bc1, bc2]
#bcs = [bc1, bc2, bc3]
normalization = extended_normalize(W, 2)
eps = Constant(0.01)

U = 0.5*(u + u_1)
U1 = 1.5*u_1 - 0.5*u_2
F = (1./dt)*inner(u - u_1, v)*dx + inner(v, dot(U1, nabla_grad(U)))*dx \
- inner(div(v), p)*dx - inner(q, div(u))*dx - inner(v, f)*dx \
+ inner(d, p)*ds + inner(c, q)*ds

# Or use + inner(d, p)*dx + inner(c, q)*dx, or inner(d, p)*ds(1) +
inner(c, q)*ds(1)
# Point is that the pressure now is completely determined

tmp = Function(W)
a, L = lhs(F), rhs(F)
t = 0.
while t < 2.5:
t += dt(0)
i = 0; error=1.
while i < 10 and error > 1e-7:
i += 1
A = assemble(a)
b = assemble(L)
[bc.apply(A, b) for bc in bcs]
tmp.vector()[:] = up_.vector()[:]
#solve(A, up_.vector(), b, 'bicgstab', 'ilu')
solve(A, up_.vector(), b)
#normalization(up_.vector())
tmp.vector().axpy(-1, up_.vector())
error = tmp.vector().norm('l2')
print 'error ', i , error
up_2.vector()[:] = up_1.vector()[:]
up_1.vector()[:] = up_.vector()[:]
print 'Time = ', t

plot(u_, title='Velocity')
plot(p_, title='Pressure')

Nice solution, but one row in A will have many non-zeros which is bad for
the condition number and thus not ideal for iterative solvers. For direct
solvers it's great though.

Mikael

On 13 March 2012 04:05, Marcos Samudio <<email address hidden>
> wrote:

> Question #190440 on CBC.PDESys changed:
>
> Marcos Samudio posted a new comment:
> I have tried the following and got acceptable results:
>
> bc1 = DirichletBC(V, (0., 0.), "on_boundary")
> bc2 = DirichletBC(V, (1., 0.), "on_boundary && x > 1. - DOLFIN_EPS")
> bc3 = DirichletBC(Q, 0., "on_boundary && x > DOLFIN_EPS && x >
> DOLFIN_EPS")
> bcs = [bc1, bc2, bc3]
>
> This way I avoid the use of extended_normalize. The error is quite
> large, but the simulation yields good looking velocity and pressure
> fields.
>
> My guess is that extended_normalize could be causing the this problem,
> so I will take a further look at it tomorrow.
>
> Marcos
>
> --
> contact for CBC.PDESys.
>

 Marcos Samudio (marcossamudio) said on 2012-03-13: #6

Both gmres and bicgstab work ok, so in fact this is something related to the direct solver.

What I do not quite get is how to set the pressure values in the original code (solvable with iterative solvers). As you stated, pressure is arbitrary up to a certain value, and I do not see when the solver assigns this arbitrary value at some given point...

Thanks,

Marcos

 Mikael Mortensen (mikael-mortensen) said on 2012-03-13: #7

On Tuesday, 13 March 2012, Marcos Samudio <
> Question #190440 on CBC.PDESys changed:
>
> Marcos Samudio posted a new comment:
> Both gmres and bicgstab work ok, so in fact this is something related to
> the direct solver.
>
> What I do not quite get is how to set the pressure values in the
> original code (solvable with iterative solvers). As you stated,
> pressure is arbitrary up to a certain value, and I do not see when the
> solver assigns this arbitrary value at some given point...

That's what normalization does, it gives pressure the average value of
zero. If you use the trick with the real functionspace then the integral of
p over the entire domain (or surface) will be zero. Most codes I know fix
the pressure to zero in one single node, which is somewhat tricky to do in
Fenics.

Mikael

>
> Thanks,
>
> Marcos
>
> --