Where this residual on SteadyNS equations comes from?

variational formulation
1

Hi everyone,
I’m new in this forum and I’m facing a problem that I’m struggling to understand.
My goal il quite simple: I’m solving the Incompressible Steady Navier Stokes Equations with no forcing trying to obtain the base flow.

from dolfin import *

Re = 40
nu = Constant(1 / Re)

mesh = Mesh()
with HDF5File(MPI.comm_world, "./Mesh.h5", "r") as h5file:
    h5file.read(mesh, "mesh", False)
    facet = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    h5file.read(facet, "facet")

VelocityElement = VectorElement("CG", mesh.ufl_cell(), 2)
PressureElement = FiniteElement("CG", mesh.ufl_cell(), 1)
Space = FunctionSpace(mesh, VelocityElement * PressureElement)

w = Function(Space)
u, p = split(w)

F_base = Function(Space)
F_base.interpolate(Constant((0.0, 0.0, 0.0)))
f_base, _ = split(F_base)

v, q = TestFunctions(Space)

bcs= []
bcs.append(DirichletBC(Space.sub(0), Constant((1.0, 0.0)), facet, 1))
bcs.append(DirichletBC(Space.sub(0), Constant((0.0, 0.0)), facet, 2))
bcs.append(DirichletBC(Space.sub(0).sub(1), Constant(0.0), facet, 3))
bcs.append(DirichletBC(Space.sub(1), Constant(0.0), facet, 4))

G = (+nu*inner(grad(u), grad(v))
     + inner(grad(u)*u, v)
     - inner(p, div(v))
     + inner(q, div(u))
     + inner(f_base, v)
     ) * dx
J = derivative(G, w)
problem = NonlinearVariationalProblem(G, w, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

By doing so, I obtain the base flow of this case (Note that I enforced a source (or forcing) term f_base to be 0 everywhere: norm(F_base) = 0.)

Now, If I simply try to compute f_base back from the Incompressible Steady Navier Stokes Equations
f_base = project( -nu * div(grad(u)) + dot(grad(u), u) + grad(p))

I get something which is no more 0 (norm(f_base != 0)).
What It’s wrong here?

Thank you all for you help.
Have a great day!

2

You do not expect the residual of the strong form of the PDE to be zero for finite element solutions. (In fact, the mismatch between the finite element solution and the strong problem is essential to many stabilized methods, which use the strong problem’s residual to modulate artificial dissipation in a way that maintains high-order accuracy.) However, it also looks like there’s a sign error in the definition of G, since the f_base term should be subtracted to have the usual interpretation as a body force.

1 Like
3

Great, thanks for the reply, got your point!.
Now, using the weak formulation of the SteadyNS (as before), I obtain the finite element solution for the base flow.
If I try to solve the same problem (using the just computed baseflow as input) to obtain the finite element solution for f_base I get the following error:

bcs2= []
#enforcing the forcing term to be 0 on every boundary
bcs2.append(DirichletBC(Space.sub(0), Constant((0.0, 0.0)), facet, 1))
bcs2.append(DirichletBC(Space.sub(0), Constant((0.0, 0.0)), facet, 2))
bcs2.append(DirichletBC(Space.sub(0), Constant((0.0, 0.0)), facet, 3))
bcs2.append(DirichletBC(Space.sub(0), Constant((0.0, 0.0)), facet, 4))

T = (+nu*inner(grad(u), grad(v))
     + inner(grad(u)*u, v)
     - inner(p, div(v))
     + inner(q, div(u))
     - inner(f_base, v)
     ) * dx
J = derivative(T, F_base)
problem = NonlinearVariationalProblem(T, F_base, bcs2, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

*** -------------------------------------------------------------------------
*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  149cf7bbfffe157a560d262289a63cfe2d1e7244
*** -------------------------------------------------------------------------

Have you any idea? Am I doing something wrong?
Thank you again

4

If the residual T is already nearly zero for the initial guess, it may be difficult to attempt to converge a second solve for F_base, since the default convergence tolerance is defined relative to the residual computed from the initial guess. Also, the Jacobian J of T with respect to F_base is going to be singular, since the derivative with respect to F_base’s “pressure” component is zero.

5

Great! Thank you… completely got the problem and your explaination has been really complete.
Have a nice day mate.