Dear all,

I am confused about the implementation of Adams Bashforth methods this Fenics tutorial on splitting schemes. See Section “Navier-Stokes equation”:

[…] and the convective term B_k^{n-\frac{1}{2}}=\mathbf{\bar{u}}\cdot \nabla \tilde u_k = (1.5 \mathbf{u}^{n-1}-0.5\mathbf{u}^{n-2})\cdot \nabla \tilde u_k is the implicit Adams-Bashforth discretization.

I implemented the splitting scheme in minimal_working_example.py (attached), which runs properly:

```
from __future__ import print_function
import math
from fenics import *
from mshr import *
import numpy as np
# from dolfin import *
import meshio
import ufl as ufl
#geometric parameters
L = 2.2
h = 0.41
r = 0.05
c_r = [0.2, 0.2]
# Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)
# 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 !
O = VectorFunctionSpace(mesh, 'P', 2, dim=2)
O3d = VectorFunctionSpace(mesh, 'P', 2, dim=3)
Q = FunctionSpace(mesh, 'P', 1)
#I will use Q4 for functions which involve high order derivatives
Q4 = FunctionSpace(mesh, 'P', 4)
# Define boundaries and obstacle
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2.2)'
walls = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'
# Define trial and test functions
#v[i] = v^i_{notes} (tangential velocity)
v = TrialFunction(O)
nu = TestFunction(O)
p = TrialFunction(Q)
q = TestFunction(Q)
#definition of scalar, vectorial and tensorial quantities
#latin indexes run on 2d curvilinear coordinates
i, j, k, l = ufl.indices(4)
#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))
T = 1.0
num_steps = 1024
dt = T / num_steps # time step size
# the Reynolds number, Re = \rho U l / \mu, Re_here = R_{notes fenics}
# Re = 1E3
mu = 0.001
print("c_r = ", c_r)
print("L = ", L)
print("h = ", h)
print("r = ", r)
print("T = ", T)
print("N = ", num_steps)
print("mu = ", mu)
# Create XDMF files for visualization output
xdmffile_v = XDMFFile("v.xdmf")
xdmf_file_p = XDMFFile("p.xdmf")
# Create time series (for use in reaction_system.py)
timeseries_v = TimeSeries("v_series")
timeseries_p = TimeSeries("p_series")
# Define functions for solutions at previous and current time steps
#this is v_{n-1}
v_n_1 = Function(O)
#this is v_{n-2}
v_n_2 = Function(O)
v_sav = Function(O)
v_ = Function(O)
#this is p_{n-1}
p_n_1 = Function(Q)
#this is p_{n-2}
p_n_2 = Function(Q)
p_ = Function(Q)
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')
bcv_inflow = DirichletBC(O, Expression(inflow_profile, degree=2), inflow)
bcv_walls = DirichletBC(O, Constant((0, 0)), walls)
bcv_cylinder = DirichletBC(O, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bc_v = [bcv_walls, bcv_cylinder, bcv_inflow]
bc_p = [bcp_outflow]
# Define expressions used in variational forms
V = 0.5 * (v_n_1 + v)
Deltat = Constant(dt)
mu = Constant(mu)
# Define variational problem for step 1
# step 1 for v
#CHOICE 1
# F1v = (dot((v - v_n_1) / Deltat, nu) \
# + (( 3.0/2.0 * v_n_1[j] * ((v_n_1[i]).dx(j)) - 1.0/2.0 * v_n_2[j] * ((v_n_2[i]).dx(j)) ) * nu[i]) \
# ) * dx \
# \
# + ( - ((nu[i]).dx(i)) * (p_n_1 + p_n_2) / 2.0 + mu * ((V[i]).dx(j) + (V[j]).dx(i) ) * ((nu[i]).dx(j)) ) * dx \
# + ( (p_n_1 + p_n_2) / 2.0 * nu[i] * ((n_inout())[i]) - mu * ( 0 * (V[i]).dx(j) + (V[j]).dx(i) ) * nu[i] * ((n_inout())[j]) ) * ds
#CHOICE 2
F1v = (dot((v - v_n_1) / Deltat, nu) \
+ (( 3.0/2.0 * v_n_1[j] - 1.0/2.0 * v_n_2[j] ) * ((V[i]).dx(j)) * nu[i]) \
) * dx \
\
+ ( - ((nu[i]).dx(i)) * (p_n_1 + p_n_2) / 2.0 + mu * ((V[i]).dx(j) + (V[j]).dx(i) ) * ((nu[i]).dx(j)) ) * dx \
+ ( (p_n_1 + p_n_2) / 2.0 * nu[i] * ((n_inout())[i]) - mu * ( 0 * (V[i]).dx(j) + (V[j]).dx(i) ) * nu[i] * ((n_inout())[j]) ) * ds
a1v = lhs(F1v)
L1v = rhs(F1v)
# step 2
F2 = ( \
- (( 1.0/2.0*(p - p_n_2) + mu * (v_[j]).dx(j) ).dx(i)) * (q.dx(i)) \
- (1.0 / Deltat) * (v_[i]).dx(i) * q \
) * dx
#THERE WOULD BE AN ADDITIONAL TERM ~ <q grad phi n>_{\partial \Omega} WHICH IS ZERO BECAUSE WE ASSUME THAT ON THE BOUNDARIES WHERE v IS SPECIFIED, WE HAVE d phi / d n = 0
a2 = lhs(F2)
L2 = rhs(F2)
# Define variational problem for step 3
# step 3 for v
F3v = (dot(v - v_ , nu) + Deltat * dot(grad( 1.0/2.0*(p_ - p_n_2) + mu*div(v_) ), nu)) * dx
a3v = lhs(F3v)
L3v = rhs(F3v)
# Assemble matrices
A1v = assemble(a1v)
A2 = assemble(a2)
A3v = assemble(a3v)
# Apply boundary conditions to matrices
[bc.apply(A1v) for bc in bc_v]
[bc.apply(A2) for bc in bc_p]
print("Starting time iteration ...", flush=True)
# Time-stepping
t = 0
for n in range(num_steps):
# Write the solution to file
xdmffile_v.write(v_n_1, t)
xdmf_file_p.write(p_n_1, t)
# Save nodal values to file
timeseries_v.store(v_n_1.vector(), t)
timeseries_p.store(p_n_1.vector(), t)
# Update current time
t += dt
# Step 1: Tentative velocity step
# step 1 for v
b1v = assemble(L1v)
[bc.apply(b1v) for bc in bc_v]
# this line solves for v^* and stores v^* in v_
solve(A1v, v_.vector(), b1v, 'bicgstab', 'hypre_amg')
# Step 2: surface_tension correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bc_p]
# this step solves for p^{n+1} and stores the solution in p_
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')
# Step 3: Velocity correction step
#compute v_n = v_{n-1}_{at the next step}
b3v = assemble(L3v)
solve(A3v, v_sav.vector(), b3v, 'cg', 'sor')
v_n_2.assign(v_n_1)
v_n_1.assign(v_sav)
p_n_2.assign(p_n_1)
p_n_1.assign(p_)
print("\t%.2f %%" % (100.0 * (t / T)), flush=True)
print("... done.", flush=True)
```

If I implement the exact same scheme as indicated in the webpage (Choice 2), the functional for the first step reads

```
#CHOICE 2
F1v = (dot((v - v_n_1) / Deltat, nu) \
+ (( 3.0/2.0 * v_n_1[j] - 1.0/2.0 * v_n_2[j] ) * ((V[i]).dx(j)) * nu[i]) ) * dx + [...]
```

This does *not* yield the same result as this Fenics sample code, for the exact same geometry and simulation parameters, nor it reproduces the turbulend, oscillatory wake behind the cylinder, but a laminar wake. The result of Choice 2 at t = 1 is shown in the image choice 2.png attached

On the other hand, I implemented the explicit Adams-Bashforth method as shown here (Choice 2). This follows the prescription of setting y_{n+1} = y_n + \Delta t [3/2 f(y_n)-1/2 f(y_{n-1})] for an ODE as d y/dt = f(y).

```
#CHOICE 1
F1v = (dot((v - v_n_1) / Deltat, nu) \
+ (( 3.0/2.0 * v_n_1[j] * ((v_n_1[i]).dx(j)) - 1.0/2.0 * v_n_2[j] * ((v_n_2[i]).dx(j)) ) * nu[i]) \
) * dx [...]
```

With choice 1, I do obtan the oscillatory, turbulent wake and the results of choice 1 agree quantitatively with those the Fenics sample code, and they are shown in choice 1.png

Why? I wonder whether the term B_k^{n-\frac{1}{2}}=\mathbf{\bar{u}}\cdot \nabla \tilde u_k = (1.5 \mathbf{u}^{n-1}-0.5\mathbf{u}^{n-2})\cdot \nabla \tilde u_k is correct.

Than you for your help !