Converting N-S example to use a mixed space

bosmacs · January 14, 2020, 3:20am

As an exercise in troubleshooting a more complex variational problem which I’m attempting to solve with fenics, I am trying to convert this incompressible N-S example to use a mixed function space to solve for velocity and pressure simultaneously, rather than using Chorin’s method.

I replaced the function space definitions with

# Define function spaces (P2-P1)
P2 = VectorElement("P", mesh.ufl_cell(), 2)
P1 = FiniteElement("P", mesh.ufl_cell(), 1)
TH = MixedElement([P2, P1])
W = FunctionSpace(mesh, TH)

# Define trial and test functions
w = TrialFunction(W)
u, p = split(w)

w = TestFunction(W)
v, q = split(w)

as well as the function definitions:

w0 = Function(W)
u0, p0 = split(w0)
w1 = Function(W)
u1, p1 = split(w1)

The weak formulation I’m attempting to implement is stated as
\begin{equation} \begin{cases} (D u_h^n, \vec{v}_h) + (\mathcal{I}(\vec{u}_h^n) \cdot \nabla \vec{u}_h^n, \vec{v}_h) + \frac{1}{Re}(\nabla \vec{u}_h^n, \nabla \vec{v}_h) - (p_{h}, \nabla \cdot \vec{v}_{h}) = (\vec{f}^n, \vec{v}_h) & \forall \vec{v}_h \in \vec{X}_{h} \\ %\label{eq:serial3}\\ (\nabla \cdot \vec{u}_h^n, r_h) = 0 & \forall r_h \in Q_{h} \end{cases}\end{equation}
where D is the discrete time derivative and \mathcal{I} is an interpolation operator; for the moment I’m using the previous timestep.

I couldn’t figure out how to solve these equations simultaneously, so I combined them in my implementation:

F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u, v)*dx + \
     nu*inner(grad(u), grad(v))*dx - inner(p, div(v))*dx - inner(f, v)*dx + div(u)*q*dx
a1 = lhs(F1)
L1 = rhs(F1)

I also changed the domain of the boundary conditions to use W.sub(0), etc. as well as the solve and assign in the time-stepping loop to use the combined vector, e.g.:

solve(A1, w1.vector(), b1, "gmres", "default")
...
w0.assign(w1)

When I try to execute the solver, my forms compile ok, but I get this error on the very first iteration:

*** -------------------------------------------------------------------------
*** 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_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  
*** -------------------------------------------------------------------------

Note that this is using the GMRES solver like the original example. If I remove “gmres” and use the default solver, I see no apparent update in the solution.

Can anyone point out anything I’m doing that’s obviously wrong? Other examples of solving coupled problems, involving N-S in particular, might also be helpful.

dajuno · January 14, 2020, 9:40am

Yes, it’s not possible to use GMRES to solve a saddle point problem without a suitable (block) preconditioner. Even then it’s a complicated problem. Better use direct solvers until you’re sure that your code runs fine.

You’d have to provide a minimal working example for us to see what’s wrong with your code using direct solvers.

For inspiration, check out:

https://fenicsproject.org/docs/dolfin/latest/python/demos/stokes-iterative/demo_stokes-iterative.py.html

Introduction to iterative solvers for saddle point problems: https://doi.org/10.1017/S0962492904000212

bosmacs · January 14, 2020, 10:27pm

Thank you for the references, the paper in particular looks like it will be quite helpful.

My full code follows, though as I noted above I have attempted to keep it as similar to the original example as possible:

from dolfin import *

# Print log messages only from the root process in parallel
parameters["std_out_all_processes"] = False;

# Load mesh from file
mesh = Mesh("lshape.xml")

# Define function spaces (P2-P1)
P2 = VectorElement("P", mesh.ufl_cell(), 2)
P1 = FiniteElement("P", mesh.ufl_cell(), 1)
TH = MixedElement([P2, P1])
W = FunctionSpace(mesh, TH)

# Define trial and test functions
w = TrialFunction(W)
u, p = split(w)

w = TestFunction(W)
v, q = split(w)

# Set parameter values
dt = 0.01
T = 3
nu = 0.01

# Define time-dependent pressure boundary condition
p_in = Expression("sin(3.0*t)", t=0.0, degree=2)

# Define boundary conditions
noslip  = DirichletBC(W.sub(0), (0, 0),
                      "on_boundary && \
                       (x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \
                       (x[0] > 0.5 - DOLFIN_EPS && x[1] > 0.5 - DOLFIN_EPS))")
inflow  = DirichletBC(W.sub(1), p_in, "x[1] > 1.0 - DOLFIN_EPS")
outflow = DirichletBC(W.sub(1), 0, "x[0] > 1.0 - DOLFIN_EPS")
bcu = [noslip]
bcp = [inflow, outflow]

# Create functions
w0 = Function(W)
u0, p0 = split(w0)  # the previous timestep solution
w1 = Function(W)
u1, p1 = split(w1)  # current pressure solution
#w1 = Function(W)
#u1, p1 = split(w1)  # current pressure solution

# Define coefficients
k = Constant(dt)
f = Constant((0, 0))

# velocity step
F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u, v)*dx + \
     nu*inner(grad(u), grad(v))*dx - inner(p, div(v))*dx - inner(f, v)*dx + div(u)*q*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Assemble matrices
A1 = assemble(a1)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

# Create files for storing solution
ufile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")

# Time-stepping
t = dt
while t < T + DOLFIN_EPS:

    # Update pressure boundary condition
    p_in.t = t

    # Compute tentative velocity step
    begin("Computing velocity and pressure")
    b1 = assemble(L1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, w1.vector(), b1) #, "gmres", "default")
    end()

    #u1, p1 = split(w1)

    # Plot solution
    plot(p1, title="Pressure", rescale=True)
    plot(u1, title="Velocity", rescale=True)

    # Save to file
    #ufile << u1
    #pfile << p1

    # Move to next time step
    w0.assign(w1)
    t += dt
    print("t =", t)

# Hold plot
import matplotlib.pyplot as plt
plt.show()

dajuno · January 15, 2020, 10:16am

You have to reassemble A1 at every time step, since a1 depends on u0.

bosmacs · January 15, 2020, 12:38pm

I added an assemble for A1, and it seems to have made no difference:

    begin("Computing velocity and pressure")
    b1 = assemble(L1)
    A1 = assemble(a1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, w1.vector(), b1) #, "gmres", "default")
    end()

Is there anything else obvious I’m missing?

dokken · January 15, 2020, 1:08pm

You are not applying the pressure conditions (bcp) to your problem which is supposed to drive your flow.

bosmacs · January 15, 2020, 10:29pm

I’m going to be kicking myself about that oversight – I overlooked those b.c. when attempting to adapt the example.

bosmacs · January 16, 2020, 1:34am

So after correcting the pressure boundary conditions, I’m seeing the solution change in time, which is progress. On the other hand, the solution seems to be pretty poor. I attempted to add a zero-mean pressure conditions by inserting

    u1, p1 = split(w1)

    # enforce zero-mean pressure condition
    p1 = p1 - assemble(p1 * dx)

into the time-stepping, but it had no noticeable effect. I will try playing with the interpolation operator, but are there any other suggestions on how to improve this direct scheme?

dokken · January 16, 2020, 5:49am

You are using a poor linearization, inner(grad(u0)*u. Have you considered solving the non-linear problem with a NewtonSolver?

bosmacs · January 20, 2020, 10:38pm

I was able to obtain a solution using a NewtonSolver as you suggested, with some significant instability in the pressure near the inlet:

I had great difficulty getting the newton iterations to converge with the problem as originally stated – specifically the inlet pressure condition that oscillates between positive and negative. I modified the problem to oscillate between a positive pressure and 0 instead to get the result pictured.

Are there other ways to improve the quality of the solution using a coupled solve as I’m trying to do here? I have tried enforcing a zero mean-pressure conditions like so:

pdofs = W.sub(1).dofmap()
w.vector()[pdofs.dofs()] -= assemble(p1 * dx)

But I’m not at all sure that’s correct. I intend to try a better linearization as well and see how that fares.

dokken · January 20, 2020, 11:27pm

I do not think that this is the correct condition to apply, as this is only used when the pressure can only be determined up to a constant. Since you use pressure Dirichlet conditions, this is not the case. However, due to the use of pressure driven flow, you need to do integration by parts carefully, as the common boundary condition related to pressure is du/dn=pn.
I will attach what I have found to be a very good source on Navier-Stokes.

bosmacs · January 28, 2020, 2:36am

Whether I apply the zero-mean pressure condition (which I agree doesn’t make sense in this setting) or not, it seems to have little effect on my solution using the Newton solver.

Thank you for the reference; it was helpful in general but did not appear to discuss the fully coupled scheme attempted here aside from a few mentions in the introduction (unless I missed it).

Are there any more fundamental problems with my usage of FEniCS? How would I apply du/dn = pn as you mentioned?

nate · January 28, 2020, 10:53pm

(\nabla {\bf u} - p \mathrm{I}) \cdot {\bf n} = {\bf g}_N, where {\bf g}_N = {\bf 0} is the natural boundary condition you should recover from integrating the elliptic term by parts. I.e: \nabla {\bf u} \cdot {\bf n} = p {\bf n}.

bosmacs · January 29, 2020, 1:47am

Thanks Nate; what I mean to be asking is, how does one apply this condition in FEniCS?

bosmacs · January 29, 2020, 3:09am

I think I may have figured it out. I added

n = FacetNormal(mesh)
I = Identity(v.geometric_dimension())

and the term

- inner(nu*dot(grad(u), n) - dot(p*I, n), v)*ds

which seems to have helped tremendously. Does this look right?

bilenemek · February 3, 2020, 8:06pm

Dear bosmacs

The so-called do nothing boundary condition looks correct.

There are many different ideas how to implement viscous fluid flow in the FEM. You may find this work useful (here is the FEniCS code)

Topic		Replies	Views
H(div) function spaces in Navier-Stokes problem General	2	114	April 13, 2024
Newton solver with mixed function space variational formulation	1	472	June 7, 2022
Some questions about the mixed function space General	1	110	November 13, 2023
Darcy law mixed formulation with give pressure gradient variational formulation	1	691	April 13, 2021
Weak formulation for square of gradient variational formulation	2	1314	August 15, 2019

Converting N-S example to use a mixed space

Related Topics