Before we get to the nonlinearity, I managed to reproduce the issue in the Fenics example of Crank Nicholson method for Navier Stokes equations, 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
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
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:
$ 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.