# Natural boundary condition is not satisfied
**URL:** https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606
**Category:** General
**Created:** [December 22, 2024, 7:47pm UTC](https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606 "2024-12-22T19:47:05Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@mekong](https://fenicsproject.discourse.group/u/mekong)
#### Post date: [December 22, 2024, 7:47pm UTC](https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606/1 "2024-12-22T19:47:05Z")
Hello,
I am implementing Navier-Stokes equations in Fenics, following [this example](https://jsdokken.com/dolfinx-tutorial/chapter2/navierstokes.html). The manifold \Omega is a rectangle, and \partial \Omega\_{\rm IN} the left boundary of \Omega. Among the boundary conditions for the pressure increment, \phi, I impose
n^i \partial\_i \phi = 0
on \partial \Omega\_{\rm IN}. I implement these BCs by setting to zero the term
\langle n^i \partial\_i \phi \, \nu\_\phi \rangle\_{\partial \Omega\_{\rm IN}} = 0
in the variational functional, where n is the normal outwards vector and \nu\_{\phi} the test function related to \phi.
However, if I then compute
\int \_{\partial \Omega\_{\rm IN}} n^i \partial\_i \phi d S
I obtain a value which is significantly different from zero.
Here is a minimal working example (mesh files are [here](https://easyupload.io/qlihox)):
```auto
from __future__ import print_function
from fenics import *
from mshr import *
import meshio
import ufl as ufl
#CHANGE PARAMETERS HERE
tol = 1E-3
L = 2.2
h = 0.41
T = 0.001
num_steps = 4
#CHANGE PARAMETERS HERE
dt = T / num_steps # time step size
mu = 0.001
rho = 1.0
i, j, k = ufl.indices(3)
#create mesh
mesh=Mesh()
with XDMFFile("triangle_mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("line_mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
O = VectorFunctionSpace(mesh, 'P', 2, dim=2)
Q = FunctionSpace(mesh, 'P', 1)
v = Function(O)
v_ = Function(O)
sigma = Function(Q)
phi = Function(Q)
nu = TestFunction(O)
q = TestFunction(Q)
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2.2)'
walls = 'near(x[1], 0) || near(x[1], 0.41)'
mf = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
#test for surface elements
ds_l = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_r = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)
ds_t = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=4)
ds_b = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=5)
xdmffile_v = XDMFFile("v.xdmf")
xdmffile_sigma = XDMFFile( "sigma.xdmf")
v_n = Function(O)
v_ = Function(O)
sigma_n = Function(Q)
sigma_ = Function(Q)
inflow_profile_v = Expression(('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0'), degree=2, h=h)
# Define symmetric gradient
def epsilon(u):
# nabla_grad(u)_{i,j} = (u[j]).dx[i]
#sym(nabla_grad(u)) = nabla_grad(u)_{i,j} + nabla_grad(u)_{j,i}
# return sym(nabla_grad(u))
return as_tensor(0.5*(u[i].dx(j) + u[j].dx(i)), (i,j))
#the normal vector on the inflow and outflow
def n_inout():
x = ufl.SpatialCoordinate(mesh)
u = as_tensor([conditional(lt(x[0], L/2), -1.0, 1.0), 0.0] )
return as_tensor(u[k], (k))
# Define boundary conditions
bcv_inflow = DirichletBC(O, inflow_profile_v, inflow)
bcv_walls = DirichletBC(O, Constant((0, 0)), walls)
bc_phi_outflow = DirichletBC(Q, Constant(0), outflow)
bc_v = [bcv_inflow, bcv_walls]
bc_phi = [bc_phi_outflow]
V = 0.5 * (v_n + v_)
Deltat = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
# step 1
F1v = rho * ( (v_[i] - v_n[i]) / Deltat + v_n[j] *((v_n[i]).dx(j)) ) * nu[i] * dx \
+ ( 2.0 * mu * (epsilon(V))[i, j] * (nu[i]).dx(j) + sigma_n * (nu[i]).dx(i) ) * dx \
+ (- sigma_n * n_inout()[i] * nu[i] - mu * (V[j]).dx(i) * n_inout()[j] * nu[i]) * ds
F2 = ( (phi.dx(i)) * (q.dx(i)) + (rho / Deltat) * ((v_[i]).dx(i)) * q ) * dx
F3v = ((v_[i] - v[i]) - (Deltat / rho) * (phi.dx(i))) * nu[i] * dx
print("Starting time iteration ...", flush=True)
# Time-stepping
t = 0
for n in range(num_steps):
t += dt
solve(F1v == 0, v_, bc_v)
solve(F2 ==0, phi, bc_phi)
solve(F3v == 0, v)
# Update previous solution
v_n.assign(v)
sigma_n.assign(-phi + sigma_n)
# Write the solution to file
xdmffile_v.write(v_n, t)
xdmffile_sigma.write(sigma_n, t)
print("BC = ", assemble(((phi.dx(i)) * (n_inout())[i])**2 * ds_l))
print("\t%.2f %%" % (100.0 * (t / T)), flush=True)
print("... done.", flush=True)
```
The natural BC is imposed in `F2` and the printout of the natural BC is in `print("BC = ", assemble(((phi.dx(i)) * (n_inout())[i])**2 * ds_l))`.
The output is
```auto
Starting time iteration ...
[...] BC = 5746230.869323557
25.00 %
[...]
BC = 3392726.0511582345
50.00 %
[...]
```
Why is this residual so large ? Thank you !
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [December 22, 2024, 8:42pm UTC](https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606/2 "2024-12-22T20:42:18Z")
First of all, why do you need to use an external mesh to illustrate this? Your mesh is a square channel, that can easily be created with a `RectangleMesh`.
Secondly, your mesh is incredibly coarse, and you should check who the solution changes with mesh refinement.
This also boils down to solving a Poisson equation on a square domain, which you should start investigating first.
Try for instance the following:
```python
from fenics import *
try:
import ufl as ufl
except ImportError:
import ufl_legacy as ufl
# CHANGE PARAMETERS HERE
tol = 1e-3
L = 2.2
h = 0.41
T = 0.001
num_steps = 1
# CHANGE PARAMETERS HERE
dt = T / num_steps # time step size
mu = 0.001
rho = 1.0
i, j, k = ufl.indices(3)
# create mesh
mesh = RectangleMesh(Point(0, 0), Point(L, h), 200, 50)
mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
eps = 100 * DOLFIN_EPS
class left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0, eps)
left().mark(mf, 2)
class right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], L, eps)
right().mark(mf, 3)
class top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], h, eps)
top().mark(mf, 4)
class bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0, eps)
bottom().mark(mf, 5)
Q = FunctionSpace(mesh, "P", 1)
# test for surface elements
ds_l = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_r = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)
ds_t = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=4)
ds_b = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=5)
inflow_profile_v = Expression(
("4.0*1.5*x[1]*(h- x[1]) / pow(h, 2)", "0"),
degree=2,
h=h,
domain=mesh,
)
bc_phi_outflow = DirichletBC(Q, Constant(0), mf, 3)
bc_phi = [bc_phi_outflow]
Deltat = Constant(dt)
phi = TrialFunction(Q)
q = TestFunction(Q)
a = (phi.dx(i) * q.dx(i)) * dx
L = -(rho / Deltat) * (inflow_profile_v[i].dx(i) * q) * dx
phih = Function(Q, name="phi")
solve(a == L, phih, bc_phi)
n = FacetNormal(mesh)
expr = dot(grad(phih), n)
print("BC = ", assemble(inner(expr, expr) * ds_l))
ff = File("phi.pvd")
ff << phih
```
which will yield
```bash
BC = 7.481641223894098e-24
```
which enforces the solution correctly.
You should now work on making this closer and closer to your example to see where you have made a mistake.
For instance, using an iterative solver (with no options specified), is rarely a good idea.
---
### Author:  [@mekong](https://fenicsproject.discourse.group/u/mekong)
#### Post date: [December 23, 2024, 7:42am UTC](https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606/3 "2024-12-23T07:42:20Z")
> [@dokken](#):
>
> First of all, why do you need to use an external mesh to illustrate this?
Because, like I said in the OP, the natural BCs are enforced only on \partial \Omega\_\textrm{IN}, thus I need to extract `ds_l` and then only way I knew to do it is by generating the mesh externally, tagging the `ds` and reading it back.
> [@dokken](#):
>
> Secondly, your mesh is incredibly coarse, and you should check who the solution changes with mesh refinement.
Does the accuracy with which a BC is satisfied depend on the mesh size ?
> [@dokken](#):
>
> You should now work on making this closer and closer to your example to see where you have made a mistake.
I will, thank you.
> [@dokken](#):
>
> using an iterative solver (with no options specified), is rarely a good idea
Why ? What options should I specify? Please keep in mind that I need to implement this for a nonlinear variational problem, in my original PDE.
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [December 23, 2024, 8:09am UTC](https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606/4 "2024-12-23T08:09:38Z")
> [@mekong](#):
>
> Does the accuracy with which a BC is satisfied depend on the mesh size ?
It is a weak enforcement, ie in an integral sense, which in turn depends on your pde.
> [@mekong](#):
>
> Why ? What options should I specify? Please keep in mind that I need to implement this for a nonlinear variational problem, in my original PDE
Navier-Stokes is a non-linear PDE, but with a splitting scheme you split it into 3 linear PDEs.
As the second step, the pressure correction, is a Poisson problem, you should likely use a multigrid method (boomeramg).
You should also ensure that your system is assembled in a symmetric fashion (ie ensure that system assembler is called), so that you can use CG.
---
### Author:  [@mekong](https://fenicsproject.discourse.group/u/mekong)
#### Post date: [December 23, 2024, 4:32pm UTC](https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606/5 "2024-12-23T16:32:11Z")
> [@dokken](#):
>
> It is a weak enforcement, ie in an integral sense, which in turn depends on your pde
So, as I already asked, does the accuracy with which a BC is satisfied depend on the mesh size ?
> [@dokken](#):
>
> you should likely use a multigrid method (boomeramg).
> You should also ensure that your system is assembled in a symmetric fashion (ie ensure that system assembler is called), so that you can use CG.
May you please give me some references on how what a multigrid ‘boomerang’ method is and how to implement it in Fenics? (And yes, I looked up `boomerang` and `multigrid` combined with Fenics online, and the results are not intelligible).
How may I assemble in a symmetric way? I found [this](https://oldqa.fenicsproject.org/114/symmetric_assemble-in-dolfin-1-2-0/) post but, like I said before, mine is a toy example for this forum, and the problem in which I will actually use this is a nonlinear one.
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [December 23, 2024, 8:47pm UTC](https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606/6 "2024-12-23T20:47:18Z")
> [@mekong](#):
>
> So, as I already asked, does the accuracy with which a BC is satisfied depend on the mesh size ?
Yes, in some sense it does. If you have too few elements to resolve the physics you would like to adhere to, it might differ slightly.
However, as I showed with the built in mesh, for the Poisson problem at hand it should be possible to resolve the Neumann condition.
> [@mekong](#):
>
> May you please give me some references on how what a multigrid ‘boomerang’ method is and how to implement it in Fenics? (And yes, I looked up `boomerang` and `multigrid` combined with Fenics online, and the results are not intelligible).
Searching for boomeramg in the dolfin repo yields:
[https://bitbucket.org/search?q=repo%3Adolfin%20Boomeramg&account={795b0eed-436d-4242-83d7-a46a0b061354}](https://bitbucket.org/search?q=repo%3Adolfin%20Boomeramg&account=%7B795b0eed-436d-4242-83d7-a46a0b061354%7D)
ie set the solver to hypre\_amg to use boomeramg.
[https://bitbucket.org/fenics-project/dolfin/src/1c52e837eb54cc34627f90bde254be4aa8a2ae17/python/dolfin/fem/solving.py?at=master#lines-117](https://bitbucket.org/fenics-project/dolfin/src/1c52e837eb54cc34627f90bde254be4aa8a2ae17/python/dolfin/fem/solving.py?at=master#lines-117)
> [@mekong](#):
>
> How may I assemble in a symmetric way? I found [this](https://oldqa.fenicsproject.org/114/symmetric_assemble-in-dolfin-1-2-0/) post but, like I said before, mine is a toy example for this forum, and the problem in which I will actually use this is a nonlinear one.
You are making it quite hard to understand your problem. You say that the splitting scheme used for Navier Stokes is your toy problem. However, it is only for the second step of your splitting scheme you see an issue. Where is your non-linearity in this equation?
Regarding symmetry, it seems like solve(F==0) uses the symmetric system assembler, so at least that is good.
I’ve already shown you an implementation of the Poisson problem (where you use the divergence of a velocity field on the rhs of the eq) that satisfies the weak condition.
One question for your problem is
> [@mekong](#):
>
> ```auto
> #the normal vector on the inflow and outflow
> def n_inout():
> x = ufl.SpatialCoordinate(mesh)
> u = as_tensor([conditional(lt(x[0], L/2), -1.0, 1.0), 0.0] )
> return as_tensor(u[k], (k))
>
> ```
Why not just use `FacetNormal(mesh)`?
---
### Author:  [@mekong](https://fenicsproject.discourse.group/u/mekong)
#### Post date: [December 29, 2024, 8:04am UTC](https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606/7 "2024-12-29T08:04:52Z")
> [@dokken](#):
>
> ie set the solver to hypre\_amg to use boomeramg.
If I understand correctly, would this be a good setting for the solver ?
```auto
J = derivative( F, phi, J_phi )
problem = NonlinearVariationalProblem( F, phi, bc_phi, J )
solver = NonlinearVariationalSolver( problem )
params = {'nonlinear_solver': 'newton',
'newton_solver':
{
'linear_solver' : 'cg',
'preconditioner' : 'hypre_amg',
}
}
solver.parameters.update(params)
solver.solve()
```
I am currently digging into the other issues that you mentioned, I will get back to you soon about it. Thank you.
---
### Author:  [@mekong](https://fenicsproject.discourse.group/u/mekong)
#### Post date: [December 31, 2024, 11:02am UTC](https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606/8 "2024-12-31T11:02:06Z")
> [@dokken](#):
>
> You are making it quite hard to understand your problem. You say that the splitting scheme used for Navier Stokes is your toy problem. However, it is only for the second step of your splitting scheme you see an issue. Where is your non-linearity in this equation?
Before we get to the nonlinearity, I managed to reproduce the issue in the Fenics [example of Crank Nicholson method for Navier Stokes equations](https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html), where the natural BCs are not satisfied. I removed the cylinder for simplicity.
In the second step, the natural BCs are \nabla \phi \cdot n = 0 on \partial \Omega\_{\rm D} = \partial \Omega\_{\rm L} \cup \partial \Omega\_{\rm T} \cup \partial \Omega\_{\rm B}, where \partial \Omega\_{\rm L} is the boundary on the Left edge of the rectangle, and so on. Here is the minimal working example
```auto
from __future__ import print_function
import csv
from fenics import *
from mshr import *
import numpy as np
import ufl as ufl
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("T")
parser.add_argument("N")
args = parser.parse_args()
T = (float)( args.T )
num_steps = (int)( args.N )
L = 2.2
h = 0.41
dt = T / num_steps # time step size
mu = 0.001
def my_norm(x):
return (sqrt(np.dot(x, x)))
# create mesh
mesh = RectangleMesh(Point(0, 0), Point(L, h), 20, 20)
mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
eps = 100 * DOLFIN_EPS
class left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0, eps)
left().mark(mf, 2)
class right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], L, eps)
right().mark(mf, 3)
class top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], h, eps)
top().mark(mf, 4)
class bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0, eps)
bottom().mark(mf, 5)
facet_normal = FacetNormal( mesh )
# Define function spaces
#the '2' in ''P', 2)' is the order of the polynomials used to describe these spaces: if they are low, then derivatives high enough of the functions projected on thee spaces will be set to zero !
P = VectorFunctionSpace( mesh, 'P', 2, dim=2 )
Q = FunctionSpace(mesh, 'P', 1)
nu = TestFunction( P )
q = TestFunction(Q)
v_n_1 = Function( P )
v_n_2 = Function( P )
v_n = Function( P )
v_ = Function( P )
p_n_12 = Function( Q )
sigma_n_12 = Function( Q )
p_n_32 = Function( Q )
phi = Function( Q )
#test for surface elements
ds_l = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_r = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=3)
ds_t = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=4)
ds_b = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=5)
def my_norm(x):
return (sqrt(np.dot(x, x)))
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2.2)'
walls = 'near(x[1], 0) || near(x[1], 0.41)'
i, j = ufl.indices(2)
csvfile = open( 'solution/natural_bcs.csv', 'w', newline='' )
fieldnames = ['n', '<(n^i \partial_i phi)^2>_l', '<(n^i \partial_i phi)^2>_{tb}', '<(n^i \partial_i phi)^2>_c']
writer = csv.DictWriter( csvfile, fieldnames=fieldnames )
writer.writeheader()
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')
bc_v__inflow = DirichletBC( P, Expression( inflow_profile, degree=2 ), inflow )
bc_v__walls = DirichletBC( P, Constant( (0, 0) ), walls )
bc_phi_outflow = DirichletBC( Q, Constant( 0 ), outflow )
bc_v_ = [bc_v __inflow, bc_v__ walls]
bc_phi = [bc_phi_outflow]
V = (v_n_1 + v_) / 2.0
Deltat = Constant( dt )
mu = Constant( mu )
# Define variational problem for step 1
F1 = (( \
((v_[i] - v_n_1[i]) / Deltat \
+ (3.0 / 2.0 * v_n_1[j] - 1.0 / 2.0 * v_n_2[j]) * (V[i]).dx( j ))) * nu[i] \
- p_n_32 * ((nu[i]).dx( i )) + mu * ((V[i]).dx( j )) * ((nu[i]).dx( j )) \
) * dx
F2 = ((phi.dx( i )) * (q.dx( i )) + (1.0 / Deltat) * (v_[i]).dx( i ) * q) * dx
F3 = ((v_n[i] - v_[i]) + Deltat * (phi.dx( i ))) * nu[i] * dx
# Time-stepping
t = 0
for n in range( num_steps ):
# Update current time
t += dt
solve( F1 == 0, v_, bc_v_ )
solve( F2 == 0, phi, bc_phi )
p_n_12.assign( p_n_32 + phi )
solve( F3 == 0, v_n )
# update the fields for the next step
v_n_2.assign( v_n_1 )
v_n_1.assign( v_n )
p_n_32.assign( p_n_12 )
#write the residual of natural BCs on step 2 to file
writer.writerows( [
{'n': n, \
'<(n^i \partial_i phi)^2>_l': sqrt(assemble( (facet_normal[i] * (phi.dx( i ))) ** 2 * ds_l ) / assemble( Constant(1.0) * ds_l )), \
'<(n^i \partial_i phi)^2>_{tb}': sqrt(assemble( (facet_normal[i] * (phi.dx( i ))) ** 2 * (ds_t + ds_b) ) / assemble(Constant(1.0) * (ds_t + ds_b))) \
}] )
print( "\t%.2f %%" % (100.0 * (t / T)), flush=True )
print( "... done.", flush=True )
```
Run it with
```auto
rm -rf solution; mkdir solution; python3 mwe.py 1.0 16
```
where `1.0` the the maximal time and `16` the number of steps. The residuals of the natural BC at each step are printed in `solution/natural_bcs.csv`, and the natural BCs are not satisfied:
```auto
$ cat solution/natural_bcs.csv
n,<(n^i \partial_i phi)^2>_l,<(n^i \partial_i phi)^2>_{tb},<(n^i \partial_i phi)^2>_c
0,14.417332724119042,0.30002951772937814,
1,13.599912917977061,0.557155701398929,
[...]
15,0.011045988527805087,0.01579581210657168,
```
This issue comes from the following: when subtracting the equation for u^n and the one for u^\ast, one obtains \partial\_i \phi = \frac{u^n\_i - u^\ast\_i}{\Delta t} + O(\Delta t). By multiplying by n^i and using the Dirichlet BC for u^n and u^\ast on \partial \Omega\_{\rm D} we are left with n^i \partial\_i \phi = O(\Delta t) on \partial \Omega\_{\rm D}. The residual in `solution/natural_bcs.csv` is, indeed, not zero, but proportional to \Delta t. It is huge in the first steps though. What do you think ?
I think that this issue should be clarified in the [webpage](https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html).