Strange behavior of Fenics sample algorithm for Navier Stokes equations

[mekong]

Hello,
I have coded this algorithm to solve the Navier Stokes equations, and I observe a very strange behavior in the solution.

I created a minimal working example for you:

from __future__ import print_function
from fenics import *
import argparse
import ufl as ufl
from mshr import *
import numpy as np
import meshio

parser = argparse.ArgumentParser()
parser.add_argument("T")
parser.add_argument("N")
args = parser.parse_args()


#latin indexes run on 2d curvilinear coordinates
i, j = ufl.indices(2)

L = 2.2
h = 0.41
r = 0.05
c_r = [0.2, 0.2, 0]

# Create mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Ellipse(Point(0.2, 0.2), 0.05, 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 64)


# Define boundaries and obstacle
#CHANGE PARAMETERS HERE
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'
#CHANGE PARAMETERS HERE

T = (float)(args.T)
num_steps = (int)(args.N)
 
dt = T / num_steps  # time step size

mu = 0.001

# Create XDMF files for visualization output
xdmffile_v = XDMFFile("v.xdmf")
xdmffile_v_ = XDMFFile( "v_bar.xdmf" )
xdmf_file_p = XDMFFile( "p.xdmf" )
xdmf_file_phi = XDMFFile("phi.xdmf")

Q_v = VectorFunctionSpace( mesh, 'P', 2, dim=2 )
Q_p = FunctionSpace( mesh, 'P', 1 )

v_n_1 = Function( Q_v )
v_n_2 = Function( Q_v )
v_n = Function( Q_v )
v_ = Function( Q_v )
p_n_1 = Function( Q_p )
p_n_2 = Function( Q_p )
p_n = Function( Q_p )
phi = Function( Q_p )
nu = TestFunction( Q_v )
q = TestFunction( Q_p )


#trial analytical expression for a vector
class TangentVelocityExpression(UserExpression):
    def eval(self, values, x):
        values[0] = 0.0
        values[1] = 0.0
    def value_shape(self):
        return (2,)


# trial analytical expression for the  surface tension p(x,y)
class SurfaceTensionExpression(UserExpression):
        def eval(self, values, x):
            # values[0] = 4*x[0]*x[1]*sin(8*(norm(np.subtract(x, c_r)) - r))*sin(8*(norm(np.subtract(x, c_R)) - R))
            # values[0] = cos(norm(np.subtract(x, c_r)) - r) * sin(norm(np.subtract(x, c_R)) - R)
            values[0] = 0.0

        def value_shape(self):
            return (1,)


# set the initial conditions for all fields
v_n_1.interpolate( TangentVelocityExpression( element=Q_v.ufl_element() ) )
v_n_2.assign(v_n_1)
p_n_1.interpolate( SurfaceTensionExpression( element=Q_p.ufl_element() ) )
p_n_2.assign(p_n_1)

inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

bcv_inflow = DirichletBC( Q_v, Expression( inflow_profile, degree=2 ), inflow )
bcv_walls = DirichletBC( Q_v, Constant( (0, 0) ), walls )
bcv_cylinder = DirichletBC( Q_v, Constant( (0, 0) ), cylinder )
bcphi_outflow = DirichletBC( Q_p, Constant( 0 ), outflow )

# boundary conditions for the surface_tension p
bc_v = [bcv_walls, bcv_cylinder, bcv_inflow]
bc_phi = [bcphi_outflow]

# Define expressions used in variational forms
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_1 + p_n_2)/2.0 * ((nu[i]).dx(i)) + mu * ((V[i]).dx(j)) * ((nu[i]).dx(j))  \
    ) * dx

# step 2
F2 = ( (phi.dx(i)) * (q.dx(i)) +  (1.0 / Deltat) * (v_[i]).dx(i) * q ) * dx

# Define variational problem for step 3
F3 = ( (v_n[i] - v_[i]) + Deltat * (phi.dx(i)) ) * nu[i] * dx


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)
    xdmffile_v_.write(v_, t )
    xdmf_file_p.write(p_n_1, t )
    xdmf_file_phi.write(phi, t)

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    solve(F1 == 0, v_, bc_v)

    # Step 2: surface_tension correction step
    solve(F2 == 0, phi, bc_phi)
    # obtain p_n from phi by using the definition of phi
    p_n.assign(2*phi + p_n_2)

    # Step 3: Velocity correction step
    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_2.assign(p_n_1)
    p_n_1.assign(p_n)

    print("\t%.2f %%" % (100.0 * (t / T)), flush=True)

print("... done.", flush=True)

which you can run with python3 example.py 0.1 128, where 0.1 is the simulation time and 128 the number of time steps in which it is divided.

The solution for the pressure p is odd: it alternates between positive and negative signs at subsequent time steps, see the screenshots attached for two time steps as an example. The same alternating behavior occurs for other time steps.

May you please explain to me what I am doing wrong here ??

[dokken]

@mekong

p_n.assign(2*phi + p_n_2)

When was this the definition of p? One has that p_n= phi + p_n_2?

[mekong]

@dokken

When was this the definition of p? One has that p_n= phi + p_n_2?

Here you say that “\phi = p^{n-1/2}- p^\ast is a pressure correction” and that " A common choice is to use p^\ast = p^{n-3/2}".

So I used these definitions of \phi and p^\ast. Thus \phi = (p^n + p^{n-1})/2 - (p^{n-1}+p^{n-2})/2 = (p^n - p^{n-2})/2, thus p_n = 2\phi + p^{n-2}

[dokken]

Note that one doesn’t describe the half steps with an averaging in the code, but one instead solve for pressure at the half step in time (ie staggered)

[mekong]

I don’t understand. May you please clarify and point out where my code is wrong ?

[mekong]

@dokken

Note that one doesn’t describe the half steps with an averaging in the code, but one instead solve for pressure at the half step in time (ie staggered)

May you please show me exactly where this is done in the code of your page? I don’t see it. Also, given that the equations are solved to linear older in \Delta t, I don’t see the difference in solving at half time step, and solving at whole time steps and then taking the average.

Thank you!

[dokken]

@mekong

Also, given that the equations are solved to linear older in \Delta tΔt\Delta t, I don’t see the difference in solving at half time step, and solving at whole time steps and then taking the average.

As you can see from your equation, you are jumping from the n-2 step to the nth step, skipping a whole time step (assuming phi is constant over two consecutive steps).

You would have to read the source code of oasisx to see the implementation.

You can see it in:
https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html#variational-form
with:

 # Step 2: Pressure corrrection step
    with b2.localForm() as loc:
        loc.set(0)
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, phi.x.petsc_vec)
    phi.x.scatter_forward()

    p_.x.petsc_vec.axpy(1, phi.x.petsc_vec)
    p_.x.scatter_forward()

[mekong]

@dokken

As you can see from your equation, you are jumping from the n-2 step to the nth step, skipping a whole time step (assuming phi is constant over two consecutive steps).

I see what you are trying to say, but this is still unclear to me: what I have done is

1. Used your definition \phi = p^{n-1/2}- p^\ast
2. Used your definition of p^\ast = p^{n-3/2}
3. Used your definition of p^{n-1/2} = (p^n+p^{n-1})/2
4. Combined them to obtain \phi = (p^n + p^{n-1})/2 - (p^{n-1}+p^{n-2})/2 = (p^n - p^{n-2})/2.

Is one of these definitions wrong? May you please show where is the mistake in steps 1, 2, 3 and 4?

Note that the definition (1) itself jumps across two time steps.

[dokken]

@mekong

Used your definition of p^{n-1/2} = (p^n+p^{n-1})/2pn−1/2=(pn+pn−1)/2

I never made this definition for the pressure.
As you are assuming a know pressure in the first step, which you are later correcting, one doesn’t do the crank Nicholson averaging for this operator. Instead one leaves the pressure at n-1/2, always solving for a time step in between the velocity steps.

[mekong]

@dokken

I never made this definition for the pressure.

Not in the webpage, but you did in this post

@mekong

Also, is p^{n-3/2} defined as (p^{n-2}+p^{n-1})/2?

@dokken

Yes, that would be the idea. I’ve now updated the webpage to reflect this.

In an effort to contribute and improve this beautiful project, I suggest that the definition of p^{n-1/2} and p^{n-3/2} and the fact that p is solved at steps 1/2, 3/2, 5/2, … while the velocity is solved at steps 1, 2, 3, … should be made crystal clear in the webpage.

I will revise my code accordingly and see whether those odd jumps in sign go away!

[dokken]

@mekong

Not in the webpage, but you did in this post

Then I was mistaken back tehn.

@mekong

In an effort to contribute and improve this beautiful project, I suggest that the definition of p^{n-1/2}pn−1/2p^{n-1/2} and p^{n-3/2}pn−3/2p^{n-3/2} and the fact that ppp is solved at steps 1/2, 3/2, 5/2, … while the velocity is solved at steps 1, 2, 3, … should be made crystal clear in the webpage.

I will take this into account once i revise the text on Oasisx.