# 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:** 1\
**Showing post:** 8

<div class="post-metadata">

**Author:** ![mekong](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mekong/32/7729_2.png) [@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")

</div>

> [@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).

---

_[View the full topic](https://fenicsproject.discourse.group/t/natural-boundary-condition-is-not-satisfied/16606)._