Overview

Chorin’s method is one of the standard ways of solving incompressible Navier Stokes. The basic idea is to explicitly compute a trial step for the velocity, then compute the pressure that corrects the divergence of the velocity, and finally apply that pressure gradient to copmute the final velocity. The FEniCS tutorial explains it well with an FEM implementation. However, the standard phrasing of the method assumes using a forward Euler time discretization, $\partial_t u \approx (u-u _0)/\Delta t$. It’s good enough for solving Navier Stokes efficiently, but what if we are trying to couple the fluid with an ODE that has different stability requirements?

How do we rephrase the method as a general Runge Kutta?

Formulation

The partial differential algebraic equation for Navier Stokes is \begin{equation} \partial_t u + \nabla u \cdot u - \mathrm{Re}^{-1} \nabla^2 u + \nabla p = f \end{equation} subject to the constraint \begin{equation} \nabla\cdot u = 0. \end{equation} The standard derivation of Chorin’s method shoves in the forward Euler approximation for the time rate immediately. I like to do that as the final step to make it easier to cycle through time steppers.

Let us define a forcing term $r$: \begin{equation} r = f - \nabla u\cdot u + \mathrm{Re}^{-1} \nabla \cdot \nabla u. \end{equation} This lets us rewrite the ODE as \begin{equation} \partial _t u = r - \nabla p \end{equation} which is still subject to the same constraint as before. We’re allowed to take the time derivative of the constraint, and exchanging with the spatial divergence gives us a new equation, \begin{equation} \partial _t \nabla \cdot u = 0 \quad \rightarrow \quad \nabla \cdot \partial _t u = 0 \quad \rightarrow \quad \nabla \cdot (r-\nabla p) = 0 \end{equation} We can use this equation to solve for $p$ initially for $r$, and then apply our time integrator to $\partial _t u$.

Now we want to phrase this in a variational formulation to weaken the derivations to apply the finite element method. The key to this reconsideration is that we recast the trial velocity step as a projection of the equation for $r$ onto the velocity basis functions. Let $v,\delta v \in \mathcal{V}$ be functions and test functions in a suitable space for velocities, and $p,\delta p \in \mathcal{P}$ be the same for the pressures.

How do we phrase this as an explicit integrator?

  1. Project $r(t)$ onto $\mathcal{V}$ by solving \begin{equation} (\delta u, r) = (\delta u, f) - (\delta u, \nabla u \cdot u) - \mathrm{Re}^{-1} (\nabla \delta u, \nabla u) \end{equation}
  2. Solve for $p(t+h)$ using $r$: \begin{equation} (\nabla \delta p, \nabla p) = (\delta p, -\nabla\cdot r) \end{equation}
  3. Apply the time stepper to the ODE: \begin{equation} (\delta u, \partial_t u) = (\delta u, r) - (\delta u, \nabla p) \end{equation} which comes out to the discrete system that can integrated normally: \begin{equation} \mathbf{M} \dot{\mathbf{u}} = \mathbf{R}(t) \end{equation}

The process is phrased as a projection step as an optimization; it would be possible to skip 1 and have steps 2 and 3 reperform the calculations inside of two different element assembly routines. However, the force coupling term $(\delta u,f)$ may be a very complicated interaction between two different numerical methods that we do not want to repeat. For example, as I’m writing this I’m implementing a particle-fluid code where that term is a discrete summation of point-integrals using my FEniCS_Particles library.

Implementation in FEniCS

We use the standard Taylor-Hood elements (quadratic velocities and linear pressures). The forms needed to perform the above steps look as follows in the FEniCS UFL:

f_v_M = inner(tu,Du)*dx
f_r_proj = - inner(tu,dot(grad(u),u))*dx \
           - mu*inner(grad(tu),grad(u))*dx
f_p_K = inner(grad(tp),grad(Dp))*dx
f_p_r = inner(tp,-div(r))*dx
f_v_dot = inner(tu, r - grad(p))*dx

This projection step requires us to implement a new class. We merge the step into the implicit pressure

class RK_field_chorin_pressure(RK_field_fenics):
    def __init__(self, r, p, f_v_M, f_r_proj, f_p_K, f_p_R, bcs=None, **kwargs):
        self.r, self.p = r, p
        self.f_r_proj = f_r_proj
        self.f_p_R = f_p_R
        self.bcs = bcs
        self.K_p = assemble(f_p_K)
        self.M_v = assemble(f_v_M)
        pyrk.RKbase.RK_field_dolfin.__init__(self, 0, [p.vector()], None, **kwargs)
    def sys(self,time,tang=False):
        solve(self.M_v, self.r.vector(), assemble(self.f_r_proj))
        R = assemble(self.f_p_R) + self.K_p*self.p.vector()
        return [R, self.K_p] if tang else R

The big addition to this class from the standard RK_field_fenics class is that it performs the $r$ projection step at the first line of sys() before return $\mathbf{R}$ and $\mathbf{K}$. Yes, there’s some weird inheritance going on where I extend a parent but use the grandparent’s constructor.

We can then initialize this special RK_field class and the standard RK_field_fenics class to have two modules representing the pressure prediction and the velocity ODE, and insert them into an explicit Runge Kutta time stepper object with any tableau we desire:

rkf_p = RK_field_chorin_pressure(r,p,f_v_M, f_r_proj, f_p_K,f_p_r,
bcs_p)
rkf_p.maxnewt = 1
rkf_v = RK_field_fenics(1, [ u ], f_v_M, f_v_dot, [], bcs_v )
Tfinal = 0.01
DeltaT = Tfinal/100.0
step = pyrk.exRK.exRK(DeltaT, pyrk.exRK.exRK_table['RK4'], [rkf_p, rkf_v] )

Setting the maxnewt field to one tells the stepper class that the pressure step is linear; it assumes everything is nonlinear.

The full implementation is located at https://github.com/afqueiruga/chorin_rk. It requires afqsrungekutta and afqsfenicsutil in the $pythonpath.

Order of accuracy

We use the lid driven cavity as a standard test problem. There are no analytical solutions, but it’s been studied to death as a benchmark because predictable vortices arise. Reference solutions can be found in “Numerical Simulation in Fluid Dynamics” by Griebel, Dornsheifer, and Neunhoeffer. Let $w$ and $h$ denote the domain size, which are both set to 1. A steady state solution we computed using RK4 with $Re=1/100$ is below:

To determine converge accuracy, the pressure and velocity fields are probed at $(0,h/4)$ where the origin is the center of the cavity. We already know that Taylor-Hood elements work, so we only use one $40\times40$ mesh for the study and only vary the time step.