Issue with sudden 'nan' error

Hi Marcos,

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

> New question #190440 on CBC.PDESys:
> https://answers.launchpad.net/cbcpdesys/+question/190440
>
> 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], 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.0)\
> or near(x[0],1.0) or near(x[1],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

F = eps*inner(grad(q), grad(p))*dx - eps*inner(grad(q), dot(grad(u), u_))*dx

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] > 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 \
    + nu*inner(grad(v), grad(U) + grad(U).T)*dx \
    - inner(div(v), p)*dx - inner(q, div(u))*dx - inner(v, f)*dx \
    - eps*inner(grad(q), grad(p))*dx \
    - eps*inner(grad(q), dot(grad(u), u_))*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
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.
>

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?

Thanks for your help.

Marcos

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] > 1. - DOLFIN_EPS")
bc3 = DirichletBC(Q, 0., "on_boundary && x[1] > DOLFIN_EPS && x[0] > 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

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

Sorry!

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] > 1. - DOLFIN_EPS")
#bc3 = DirichletBC(Q, 0., "on_boundary && x[1] < DOLFIN_EPS && x[0] <
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 \
    + nu*inner(grad(v), grad(U) + grad(U).T)*dx \
    - inner(div(v), p)*dx - inner(q, div(u))*dx - inner(v, f)*dx \
    - eps*inner(grad(q), grad(p))*dx \
    - eps*inner(grad(q), dot(grad(u), u_))*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:
> https://answers.launchpad.net/cbcpdesys/+question/190440
>
> 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] > 1. - DOLFIN_EPS")
> bc3 = DirichletBC(Q, 0., "on_boundary && x[1] > DOLFIN_EPS && x[0] >
> 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
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.
>

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

On Tuesday, 13 March 2012, Marcos Samudio <
<email address hidden>> wrote:
> Question #190440 on CBC.PDESys changed:
> https://answers.launchpad.net/cbcpdesys/+question/190440
>
> 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
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.
>

Ok, I get it now! I was trying to do exactly what you say, fixing pressure to zero at the corner node. I will take a closer look at the normalization function.

Thanks again Mikael