Documentation 229 FV Framework

2.
29 Finite Volume MATLAB Framework Documentation
Matt Ueckermann, Jing Lin, Sydney Sroka, Pierre Lermusiaux
April 12, 2018
1
Contents
1 Problem Formulations: Incompressible N-S Equations 3
1.1 General 3-D N-S Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 N-S Equations with Boussinesq Approximation . . . . . . . . . . . . . . . . . . . . . 3
1.3 2-D Flow Idealization and Nondimensionalization . . . . . . . . . . . . . . . . . . . . 5
1.4 Variable Correspondence in the Finite Volume Framework . . . . . . . . . . . . . . . 6
1.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.5.1 Bottom Gravity Current . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Discretization: Schemes and Implementation 8

2.1 Advection Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.1.1 Central Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.2 Upwind Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.3 QUICK Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.4 TVD Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.5 Advection Schemes for Dynamically Orthogonal Equations . . . . . . . . . . 11
2.1.6 Mexed Routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.7 The CFL Condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Matrix Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.1 Gradient Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.2 Divergence Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.3 Laplace Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 Symmetry of the Discrete Laplace Operator . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Consistency of Discrete Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.5 Stokes Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.6 Solving Large Sparse Linear Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.7 Solving Poisson Equations with All Neumann Boundary Conditions . . . . . . . . . 13
3 Download and Installation 18

3.1 Downloading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Installing on Linux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.3 Installing on Mac . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.4 Installing on Windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4 Quick Start Guide 20

4.1 Brief Description of Provided Examples . . . . . . . . . . . . . . . . . . . . . . . . . 21
5 Setting Up New Examples 22

5.1 Naming Conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.2 Node Numbering Convention . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.3 Data Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.4 Domain, Masking and Boundary Conditions: General Rules and Guidelines . . . . . 26
5.4.1 Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
5.4.2 Masking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
5.4.3 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
5.5 Stokes Forcing Test Case: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.5.1 Setting Up Forcing Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2
5.5.2 Setting Up the Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.5.4 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
5.5.5 Misc. and Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
5.6 Sudden Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
5.6.1 Setting Up Forcing Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.6.2 Setting Up the Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.6.4 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.6.5 Misc. and Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
6 Troubleshooting 41
1 Problem Formulations: Incompressible N-S Equations
1.1 General 3-D N-S Equations
Before going into 2-dimensional flows which this finite volume framework is intended for, we first
state the governing equations in a general 3-dimensional setting so that users can have an idea of
in what ways a 3-D flow can be simplified to a 2-D one and be simulated by this code.
The general 3-D incompressible Navier-Stokes equations are:
∇ · u = 0, (1)
∂u 1
+ u · ∇u = − ∇p + ν∇2 u + g − 2Ω × u + F u (x, t), (2)
∂t ρ
∂φ
+ u · ∇φ = κφ ∇2 φ + Fφ (φ, x, t). (3)
∂t
Note that all boldface symbols represent 3-D vectors here.
u is the fluid velocity field while p is the fluid pressure field. ρ and ν, which are density and kinetic
viscosity, respectively, are given constant physical properties of the fluid. g is a given constant grav-
itational acceleration. Ω(x, t) is a given angular velocity vector field under the reference framework
(e.g. earth rotation). F u is a given source term representing other known exterior forces.
φ is an arbitrary tracer, namely a scalar field transported by advection and diffusion. κφ is a given
constant diffusivity of φ and Fφ is a given source term which can also depend on φ. For example, if
φ is the temperature, which is a measure of internal energy, then κφ is the thermal diffusivity and
Fφ represents heat sources and sinks due to, for instance, radiation and chemical reactions.
1.2 N-S Equations with Boussinesq Approximation
In general (1) and (2) already form a closed system and (3) is coupled to them only in a one-way
manner. However, there is one situation frequently-encountered in nature where a tracer also goes
into the momentum equation (2) and thus makes (1), (2) and (3) fully coupled. This is the case
called density driven flow such as natural convection in atmosphere and bottom gravity current in
ocean.
3
In these flows, the dominant driving effect is the buoyancy due to non-uniformly distributed fluid
density together with gravity. The density variation could be mainly due to either temperature
or salinity or even both, but hardly pressure. Moreover, the density fluctuation is usually small
compared to its absolute value (e.g. < %1 in ocean). In such cases, density variation is important
and can not be ignored completely, but on the other, compressibility effect, namely pressure wave,
is still immaterial and thus negligible. Hence, we want to approximate the flow in such a way that
it can still be treated as incompressible but at the same time, the effects of gravity exerted on fluid
with different densities are accounted.
The Boussinesq approximation achieves this goal. First, rewrite the momentum equation (2) for
incompressible flows as below:

∂u
ρ + u · ∇u = −∇p + ρν∇2 u + ρg − 2ρΩ × u + ρF u (x, t). (4)
∂t
Since the density ρ only has a small variation, we decompose it into the sum of a reference value
ρ0 and a small perturbation ρ0 :
ρ = ρ0 + ρ0 . (5)
Accordingly, we decompose the pressure p into the sum of an inactive part p0 1 due to ρ0 and an
active part p0 :
p = p0 + p0 (6)
where
p0 (x) = p0,ref + ρ0 g · (x − x0,ref ). (7)
Here x0,ref is an arbitrary reference point and p0,ref is an arbitrary reference pressure.
With (5), (6) and (7), the pressure and gravity terms in the momentum equation (4) can be
rearranged as
−∇p + ρg = −∇(p0,ref + ρ0 g · (x − x0,ref ) + p0 ) + (ρ0 + ρ0 )g = −∇p0 + ρ0 g (8)
Substituting (5) and (8) into (4) and dividing by ρ0 gives
ρ0 ρ0 ρ0

∂u 1 0
ν∇2 u − 2Ω × u + F u (x, t) .

1+ + u · ∇u = − ∇p + g + 1 + (9)
ρ0 ∂t ρ0 ρ0 ρ0
The Boussinesq approximation says if |ρ0 | ρ0 holds in a flow, the effects of ρ0 are negligible in all
terms but the gravity one. Hence (9) can be well approximated by
∂u 1 ρ0
+ u · ∇u = − ∇p0 + ν∇2 u + g − 2Ω × u + F u (x, t). (10)
∂t ρ0 ρ0
Since ρ0 is just a constant, for convenience we usually define a reduced pressure and a reduced
density:
p0 ρ0 g
P0 = , ρ̃0 = . (11)
ρ0 ρ0
1
Some authors (e.g. Cushman-Roisin and Beckers (2011) p. 85) call p0 hydrostatic pressure, which is actually
confusing. The term ”hydrostatic pressure” usually refers to another quantity that has a well-defined physical
meaning, which we will see later. The p0 here is more of a helper variable introduced just to simplify the equations
and p0 does not have realistic physical meaning.
4
Note that by Boussinesq approximation, the continuity equation (1) and tracer transport equation
(3) are unchanged.
For a more detailed discussion on the Boussinesq approximation, see Cushman-Roisin and Beckers
(2011) Sec. 3.7.
1.3 2-D Flow Idealization and Nondimensionalization
A 2-dimensional flow in a horizontal plane is governed by

∂u ∂v
+ = 0,
∂x ∂y
∂u ∂u ∂u ∂P ∂2u ∂2u
+u +v =− + ν 2 + ν 2 + f (x, y, t)v + Fu (x, y, t)
∂t ∂x ∂y ∂x ∂x ∂y
(12)
∂v ∂v ∂v ∂P ∂2v ∂2v
+u +v =− + ν 2 + ν 2 − f (x, y, t)u + Fv (x, y, t)
∂t ∂x ∂y ∂y ∂x ∂y
∂ρ ∂ρ ∂ρ 2
∂ ρ 2
∂ ρ
+u +v = κ 2 + κ 2 + Fρ (ρ, x, y, t).
∂t ∂x ∂y ∂x ∂y
A 2-dimensional flow in a vertical plane is governed by

∂u ∂w
+ = 0,
∂x ∂z
∂u ∂u ∂u ∂P 0 ∂2u ∂2u
+u +w =− + νxy 2 + νz 2 + Fu (x, z, t)
∂t ∂x ∂z ∂x ∂x ∂z
(13)
∂w ∂w ∂w ∂P 0 ∂2w ∂ 2 w ρ0 g
+u +w =− + νxy 2 + νz 2 − + Fw (x, z, t)
∂t ∂x ∂z ∂z ∂x ∂z ρ0
∂ρ0 ∂ρ0 ∂ρ0 ∂ 2 ρ0 ∂ 2 ρ0
+u +w = κxy 2 + κz 2 + Fρ (ρ0 , x, z, t).
∂t ∂x ∂z ∂x ∂z
If the hydrostatic pressure is removed, the equations become (with a fixed surface at z = b):
∂u ∂w
+ = 0,
∂x ∂z
∂P 0
Z b 0
∂u ∂u ∂u g ∂ρ ∂2u ∂2u
+u +w =− − dz + νxy 2 + νz 2 + Fu (x, z, t)
∂t ∂x ∂z ∂x ρ0 z ∂x ∂x ∂z
0 2 2
(14)
∂w ∂w ∂w ∂P ∂ w ∂ w
+u +w =− + νxy 2 + νz 2 + Fw (x, z, t)
∂t ∂x ∂z ∂z ∂x ∂z
∂ρ0 ∂ρ0 ∂ρ0 2
∂ ρ 0 2
∂ ρ 0
+u +w = κxy 2 + κz 2 + Fρ (ρ0 , x, z, t).
∂t ∂x ∂z ∂x ∂z
5
1.4 Variable Correspondence in the Finite Volume Framework
In this finite volume MATLAB framework, the following equations are solved:
∂u ∂v
+ = 0,
∂x ∂y
∂u ∂u ∂u ∂P ∂2u ∂2u
+u +v =− + ν1 2 + ν2 2 + f (x, y, t)v + Fu (x, y, t)
∂t ∂x ∂y ∂x ∂x ∂y
2 (15)
∂v ∂v ∂v ∂P ∂ v ∂2v
+u +v =− + ν1 2 + ν2 2 − ρg − f (x, y, t)u + Fv (x, y, t)
∂t ∂x ∂y ∂y ∂x ∂y
∂ρ ∂ρ ∂ρ 2
∂ ρ 2
∂ ρ
+u +v = κ1 2 + κ2 2 + Fρ (ρ, x, y, t).
∂t ∂x ∂y ∂x ∂y
The code accepts the parameters and variables through the following relations:
app.nu = ν1 , app.nu2 = ν2 , app.kappa = κ1 , app.kappa2 = κ2 ,

app.g = g, app.Fcor = f (x, y, t), (16)
u = u, v = v, P = P, rho = ρ.
1.5 Examples
1.5.1 Bottom Gravity Current
Based on Özgökmen and Chassignet (2002), if we assume the flow is isothermal and the variation
of density is entirely due to salinity with Boussinesq approximation, the governing RANS equations
of the flow read
∂u ∂w
+ = 0
∂x ∂z
∂u ∂u ∂u ∂P 0 ∂2u ∂2u
+u +w = − + νx 2 + νz 2
∂t ∂x ∂z ∂x ∂x ∂z
∂w ∂w ∂w ∂P 0 ∂2w ∂ 2 w ρ0 g
+u +w = − + νx 2 + νz 2 −
∂t ∂x ∂z ∂z ∂x ∂z ρ0
∂S 0 ∂S 0 ∂S 0 2
∂ S 0 2
νz ∂ S 0
+u +w = Kx 2 + Kx − k(S 0 − ∆S).
∂t ∂x ∂z ∂x νx ∂z 2
The density variation and salinity variation are defined as ρ0 = ρ − ρ0 , S 0 = S − S0 , respectively and
are related by ρ0 = ρ0 βS 0 (in accord with previous sections). Furthermore, p∗ is the non-hydrostatic
part of the mechanical pressure divided by the reference density ρ0 , i.e. P 0 = pm /ρ0 + gz (again,
as in previous sections). Last but not least, there is a relaxation term in the salinity equation
representing a salinity forcing where k is the relaxation frequency and ∆S is a given maximum
salinity variation. This term is only turned on near the inlet.
We assume the flow is going from left to right and down a slope. The boundary conditions are set
as below. At the bottom, the no-slip condition is used:
∂S 0
u = 0, w = 0, = 0.
∂n
6
At the top, the free-slip condition is used:
∂u ∂S 0
= 0, w = 0, = 0.
∂z ∂z
Both the inlet and the outlet use an open boundary condition
∂u ∂w ∂S 0
= 0, = 0, = 0.
∂x ∂x ∂x
Note that there is no driving effect in the boundary conditions. The gravity flow is driven by the
salinity forcing term turned on only in a region near inlet.
Now we first nondimensionalize the equation to facilitate both understanding the physics and
solving the equations numerically. We use the thickness of the gravity current h as a scale for length
and νx /h as a scale for velocity. Thus the time is scaled as h2 /νx and the pressure by (νx /h)2 .
Besides, we use ∆S to scale S 0 . Therefore, we have the following nondimensional variables:
hu hw h2 P 0 S0
u∗ = , w∗ = , P 0∗ = , S 0∗ = ,
νx νx νx2 ∆S
h2 t x z
t∗ = , x∗ = , z∗ = .
νx2 h h
After the nondimensionalization, the governing equations read

∂u∗ ∂w∗
+ = 0,
∂x∗ ∂z ∗
∂u ∗
∗ ∂u
∗
∗ ∂u
∗ ∂P 0 ∗ ∂ 2 u∗ ∂ 2 u∗
+ u + w = − + + r ,
∂t∗ ∂x∗ ∂z ∗ ∂x∗ ∂x∗ 2 ∂z ∗ 2
∂w ∗
∗ ∂w
∗
∗ ∂w
∗ ∂P 0 ∗ ∂ 2 w∗ ∂ 2 w∗ ∗
+ u + w = − ∗ + + r − GrS 0 ,
∂t∗ ∂x∗ ∂z ∗ ∂z
2 0 ∗∂x
∗ 2 ∂z ∗ 2
∗ ∗ 0∗
∂S 0
∗ ∂S
0
∗ ∂S 1 ∂ S ∂2S0∗ ∗
+ u + w = + r ∗ 2 − k ∗ (S 0 − 1).
∂t∗ ∂x∗ ∂z ∗ P r ∂x ∗ 2 ∂z
Four nondimensional numbers appear in the system:
νz gh3 β∆S νx kh2

r= , Gr = , Pr = , k∗ = .
νx νx2 Kx νx
r is the ratio of vertical to horizontal turbulent diffusivity. Gr is the Grashof number. P r is the
turbulent Prandtl number and k ∗ is the nondimensional relaxation frequency.
The boundary conditions in terms of non-dimensional variables are as below. At the bottom, the
no-slip condition is used:
∂S 0 ∗
u∗ = 0, w∗ = 0, = 0.
∂n
At the top, the free-slip condition is used:
∂u∗ ∂S 0 ∗
= 0, w∗ = 0, = 0.
∂z ∗ ∂z ∗
7
Both the inlet and the outlet use an open boundary condition
∂u∗ ∂w∗ ∂S 0 ∗
= 0, = 0, = 0.
∂x∗ ∂x∗ ∂x∗
Furthermore, the initial condition is set as

∗
u∗ = 0, w∗ = 0, S0 = 0
except at the salinity forcing region we have

∗ 1 1km − z
S0 = 1 − cos π .
2 1km − 0.8km
The code accepts the parameters and variables through the following relations:
1 r
app.nu = 1, app.nu2 = r, app.kappa = , app.kappa2 = , app.g = Gr,
Pr Pr
u = u∗ , v = w∗ , P = P 0∗ , ρ = S 0∗ .
2 Discretization: Schemes and Implementation
This section summarizes the numerical discretization utilized in the framework. Again, this is not
an exhaustive manual, but should be sufficient for students who have been introduced to numerical
fluid mechanics. In general the MATLAB functions are well commented, and users are encouraged
to read the codes for further study.
2.1 Advection Schemes
The different advection schemes implemented are described in this section. The functions discussed
can be found under <app root>/Src/Advect>, where <app root> is the root folder where these
MATLAB scripts were uncompressed.
The central, upwind and QUICK schemes for tracer and u, v-velocity advection are implemented.
The discrete finite volume problem is as follows:
∂φx,y n o
∆V = ux−∆x/2,y φ̂x−∆x/2,y − ux+∆x/2,y φ̂x+∆x/2,y + (17)
∂t n o
vx,y−∆y/2 φ̂x,y−∆y/2 − vx,y+∆y/2 φ̂x,y+∆y/2
Or using the East, West, North, South naming convention:

∂φx,y n o n o
∆V = uW φ̂w − uE φ̂e + vS φ̂s − vN φ̂n (18)
∂t
where the problem now is Rselecting appropriate values of the fluxes φ̂x±∆x/2,y±∆y , which are ap-
Ai φdAi
proximately equal to φ̂i ≈ Ai .
8
φ +φx,y
A logical choice would be φ̂x+∆x/2,y = x+∆x,y2 and this is known as the central scheme. How-
ever, in general, the central scheme is unstable.
For the following sections, we will use the east face flux as an example.
2.1.1 Central Scheme
The central difference scheme (CDS) uses a linear interpolation to determine the value of the
function at interfaces:
φP + φE
φ̂e = . (19)
2
The central scheme is implemented in SCAadvect CDS mex.m and UVadvect CDS mex.m.
2.1.2 Upwind Schemes
The upwind scheme works by choosing the upwind value of the function as the true value of the
function at the interface of a control volume. The discrete flux function then takes the value:

φP uE > 0
φ̂e = (20)
φE uE < 0
To implement this in the code, a trick involving absolute values is used to avoid using ”if” state-
ments:
1
φ̂e = {uE (φE + φP ) − |uE |(φE − φP )} (21)
2
The upwind tracer and u, v-velocity advection functions are implemented in SCAadvect UW.m and
UVadvect UW.m respectively.
For the tracer advection functions, the velocities at the center of the control surfaces are conveniently
available, hence only the upwind value of the function needs to be selected in two dimensions for
each of the four surfaces. To to this, the indices were carefully selected using integers ”w, e, s, n”
for the west, east, south, and north faces. Setting one of these integers equal to 1 and the rest to 0
then selects the correct indices for the corresponding face. Therefore, the only change in the code
from one flux to the next is setting these integers, as well as selecting either the u or v component
of velocity.
The final line selects only the active interior nodes, and returns the flux for those nodes. The
advection is performed for all nodes in the domain, including inactive boundary nodes in the
interior if present.
Essentially the same procedure is followed in the function performing the u, v-velocity advection,
however in this case we need to consider the staggering of the grid more carefully. That is, sometimes
it is necessary to average the vertical or horizontal velocity to obtain a value at the desired location.
The code for the u and v velocity advection is very similar, in fact, once the u-velocity code was
working it was copied, pasted, and with minor modifications was used for the v-velocity.
9
2.1.3 QUICK Scheme
The QUICK scheme interpolates the value of the function at the midpoint of the control surfaces
using a quadratic interpolation. This is similar to the unstable central flux scheme which uses a
linear interpolation, with the difference that an upwind bias is introduced. Since three points are
required for the interpolation, two points are chosen from the upwind direction, and only one is
chosen from the downwind direction. The QUICK scheme is then given by:

6/8φP − 1/8φW + 3/8φE uE > 0
φ̂e = (22)
6/8φE − 1/8φEE + 3/8φP uE < 0
Or using the naming convention upwind = U, up-upwind = UU, downwind = D:
φ̂e = 6/8φU − 1/8φU U + 3/8φD (23)
The code is essentially the same as the upwind case. The major difference is that first the interpo-
lated values are found. Then, instead of using the upwind and downwind values of the function given
by the CV centers, the interpolated upwind and downwind values of the functions are used.
The exact same procedure is used for the u, v-velocity advection, except (again) care must be taken
with the staggering of the control volumes. This scheme is implemented in SCAadvect mex.m and
UVadvect DO mex.m.
2.1.4 TVD Scheme
The TVD scheme uses a mixture between the CDS and UW schemes. TVD stands for Total
Variation Diminishing. Basically, what happens is the slope is ”‘limited”’ to ensure that the
maximum value of the function is never exceeded, and the minimum value of the function is never
undershot.
We use a standard Total Variation Diminishing (TVD) scheme, with an monotonized central (MC)
symmetric flux limiter (Van Leer, 1977). The scheme can be written for a variable η as:
ηi + ηi−1
F (ηi− 1 ) = ui− 1
2 2 2
η −η
∆t
(24)
i i−1
− ui− 1 1 − 1 − ui− 1 Ψ(r 1) ,

i− 2
2 2 2 ∆x
where the MC slope limiter Ψ(r) is defined as

1+r
Ψ(r) = max 0, min min , 2 , 2r ,
2
and the variable r as
h i
1
ui− 1 + ui− 1 (ηi−1 − ηi−2 ) + 12 ui− 1 − ui− 1 (ηi+1 − ηi )

2 2 2 2 2
ri− 1 = ,
2 ui− 1 (ηi − ηi−1 )
2
where ui− 1 is used without interpolation for the density advection (due to the grid staggering)
2
while a second-order linear interpolation is used for the non-linear u and v advection. For more
information about TVD schemes, we refer to LeVeque (2002).
10
The exact same procedure is used for the u, v-velocity advection, except (again) care must be taken
with the staggering of the control volumes. This scheme is implemented in SCAadvect TVD mex.m
and UVadvect TVD mex.m.
2.1.5 Advection Schemes for Dynamically Orthogonal Equations
For the implementation of the stochastic dynamically orthogonal (DO) schemes, we require more
flexibility in the advection schemes. For this reason, the DO specific routines UVadvect DO.m,
UVadvect DO UW.m and UVadvect DO CDS.m are provided. Basically, functions allow for either of
the following choices for the u advection.
∂uj ∂uj
ui + vi (25)
∂x ∂y
∂vj ∂vj
ui + vi (26)
∂x ∂y
2.1.6 Mexed Routines
To increase the efficiency of the solvers, it was found that the MATLAB implementations of the
advection schemes was slow compared to a C++ implementation. Therefore, most of the advection
schemes were coded in C, compiled, and linked to MATLAB (mexing).
Note, that you may have to compile the binaries yourself. This is a simple matter on a Linux or
Mac system with a
>> mex filename.cpp
whereas on a Windows system, first an appropriate C++ library needs to be installed and configured
to work with MATLAB. Refer to the Installation section for more details on this.
2.1.7 The CFL Condition
The Courant-Friedrichs-Lewy (CFL) condition is a stability criterion on the numerical solution of

the advection equations. This condition states that the numerical domain of dependence must
contain the physical domain of dependence. In other words, the maximum allowable information
propagation speed in a numerical scheme must exceed the physical advection speed. For an explicit
scheme for the advection terms, the CFL condition requires that
∆t
C ≤ 1, (27)
∆x
where C is the physical advection speed.
If this condition is violated, the numerical solution may become unstable. Therefore, if a simula-
tion blows up, we may need a smaller time step ∆t.
2.2 Matrix Operators
To solve the Navier-Stokes equations using the projection methods as used in this code, the Lapla-
cian and diffusion matrices need to be inverted, and gradients and divergences need to be taken.
11
The divergence and gradient operations do not require inversion, and are therefore implemented as
functions in a matrix-free format. The diffusion matrix is easily build from the Laplacian by adding
a diagonal component, hence the major matrix that needs to be built and stored is the Laplacian
matrix.
2.2.1 Gradient Operator
For solving the Navier-Stokes and Stokes equations, only the gradient of the pressure is required.
This is implemented in Grad2D.m. A simple second order central difference (finite difference) scheme
is used to find the gradient of the pressure:

∂φ φ(x + ∆x) − φ(x)
= (28)
∂x x+ ∆x ∆x
2
∂φ φ(y + ∆y) − φ(y)
= (29)
∂y y+ ∆y
∆y
2
The gradient is calculated for the entire domain, including any interior masked points, and then
only the active interior point are selected and outputted.
2.2.2 Divergence Operator
For solving the Navier-Stokes and Stokes equations, only the divergence of the velocity is required.
This is implemented in Div2D.m.The basic divergence operator implementation is essentially the
same as the Gradient operator, since it also use a second-order accurate finite difference approx-
imation of the derivatives. However, in the case where Neumann boundary conditions are used,
∂φ
additional steps are required to set the correct value of ∂n at the boundary. Therefore, in the code,
first the derivatives are found while ignoring the boundary conditions. Then the IDs of the Neu-
mann boundaries are found. Finally, on Neumann boundaries, the calculated derivative is replaced
by the specified values of the Neumann boundaries.
The implementation looks, perhaps, more complicated than it should be. However, this implemen-
tation is generally applicable to any interior masked domain.
2.2.3 Laplace Operator
The discrete version of the Laplacian operator ∇2 is built using the mk DiffusionOperator func-
tion. The base discrete matrix is built for either Neumann or Dirichlet boundary conditions.
However, due to staggering, the Dirichlet boundary conditions will be incorrect for some faces of
the u, v CVs, and all faces of the Pressure CVs. These are fixed using FixDiffuvDbcs. Also,
open boundary conditions are corrected for using FixDiffPObcs (open boundaries set the second
derivative equal to zero).
12
The second derivative is approximated using the following:
2
∂ φ ∂2φ

−φ(x + δx , y) + 2φ(x, y) − φ(x − δx , y)
− 2
+ 2 = (30)
∂x ∂y x,y ∆x2
−φ(x, y + δy ) + 2φ(x, y) − φ(x, y − δy )
+
∆y 2
note that the negative of this operator is implemented in the code.
2.3 Symmetry of the Discrete Laplace Operator
2.4 Consistency of Discrete Operators
2.5 Stokes Solvers
We use a projection method to solve the Stokes problem. That is, the time integration proceeds in
three steps, and is as follows:
uk

I
− ν∇ ũk+1 =
2
+ Fk+1 (31)
∆t ∆t
1
∇2 P k+1 = ∇ · ũk+1 (32)
∆t
uk+1 = ũk+1 − ∆t∇P k+1 (33)
where P in this case is only the “Pseudo-pressure,” and does not actually solve the correct pressure.
An improved algorithm is the Incremental Pressure correction scheme, and we use the rotational
form:
uk

I
− ν∇ ũk+1 =
2
− ∇P k + Fk+1 (34)
∆t ∆t
1
∇2 (q k+1 ) = ∇ · ũk+1 (35)
∆t
uk+1 = ũk+1 − ∆t∇q k+1 (36)
P k+1 = q k+1 + P k − ν∇ · ũk+1 (37)
For the Navier Stokes equations we treat the nonlinear terms explicitly. That is, to solve the Navier-
Stokes equations, we set Fk+1 ≈ −uk · ∇uk . For an excellent review of the projection methods, see
Guermond and Minev (2006).
2.6 Solving Large Sparse Linear Systems
2.7 Solving Poisson Equations with All Neumann Boundary Conditions
In the incompressible Navier-Stokes equations, the pressure is determined only up to a constant.

In discretizing the equations, this fact leads us to solving a Poisson equation with all Neumann
boundary conditions, whose solution is also determined up to a constant, when we either directly
13
solve the pressure Poisson equation or solve the derived Poisson equation for a pressure correction
variable in the projection methods.
The problem can be formulated in general as
∇ · ∇φ(x) = f (x), x ∈ Ω, (38)
∂φ(x)
= g(x), x ∈ ∂Ω, (39)
∂n
which has a solution when the compatibility condition
Z I
f dV = gdS (40)
Ω ∂Ω
is satisfied. This is because

Z Z I I I
∂φ
f dV = ∇ · ∇φdV = n · ∇φdS = dS = gdS. (41)
Ω Ω ∂Ω ∂Ω ∂n ∂Ω
The second equality comes from the Gauss’s Theorem and the third from the definition of a direc-
tional derivative. n is the unit outward normal vector at ∂Ω. When a solution exists, the solution
is not unique. The classical mathematical theory on Poisson equations guarantees us that the
solution, when it exists, is determined up to a uniform constant over Ω.
When we discretize the Poisson equation and the Neumann boundary condition and form a linear
system Ax = b, we would expect that the singular nature of this problem will yield a singular
coefficient matrix A and that the compatibility condition will yield a b that lies in the image of
A so that this linear system has an infinite number of solutions. Since φ is determined up to a
uniform constant over Ω, we would also expect that the null space (kernel) of A is one-dimensional
and is spanned by the uniform vector e = [1, . . . , 1]T .
However, our first question is: can these expectations really be met? That is, are the properties
of the continuum problem automatically inherited by its discrete counterpart, regardless of the
discretization schemes? If not, what are the constraints on the discretization that will guarantee
the transfer of the aforementioned properties to the discrete problem? One last question is: given
a discrete problem, which is a singular linear system, how can we stably, accurately and efficiently
obtain a specific solution? Note that even if the above expectations can be met analytically, round-
off errors can still ruin these properties and cause the algorithms to break down.
Although these questions sound really basic and there should have been quite standard answers to
them, to our surprise, we find really few expositions on these questions in textbooks and litera-
ture.
The reader might wonder: why not least squares? It is true that there are algorithms out there to
solve rank-deficient least-squares problems, which seems to provide an answer to our last question
and bypass all the others. However, this is not the case due to our context of solving PDEs. The
linear systems coming from discretizing PDEs are large and sparse ones, which makes unaffordable
any algorithms that do not take into account the sparsity to reduce the time and space complexity.
Unfortunately, QR and SVD, upon which almost all practical least-squares algorithms are based,
are extremely sparsity-unfriendly, because they require forming orthogonal matrices of the same
size as A. These orthogonal matrices are in general not sparse and have very few nonzero entries.
If we assume the size of A to be N -by-N , the storage requirement will be of order O(N 2 ), which
14
is prohibitively large. Moreover, the time complexity of QR and SVD is of order O(N 3 ), which is
again unacceptable.
With the QR- and SVD-based methods ruled out, we can only seek to manipulate A in a way that
the manipulated matrix A0 is invertible and the unique solution to A0 x = b is a specific solution to
Ax = b. Ideally, we also want to keep A0 as sparse as possible and let A0 share as many properties
with A as possible, so that we can apply our usual sparse linear solver to this system. Here we
assume the aforementioned properties of the continuum problem are all inherited by the discrete
one. Later we just to need to verify that the algorithms will not break down even when these
assumptions are not strictly true due to round-off errors.
The first idea is equation replacement or addition. Since the null space of A is one-dimensional,
there is one, though non-unique, redundant equation in Ax = b. The possible choices of this
redundant equation is certain. The singularity of A implies that there exists a subset of rows with
one nonzero coefficient associated with each of these rows such that this linear combination sums
to 0. In other words, each of these rows can be expressed as a linear combination of the others,
which makes it redundant provided the others. Note that this subset of rows is unique. If we will
only modify one row, this row has to be in this subset. For the case with all Neumann boundary
conditions, this subset is actually the full set of rows based on empirical results. However, if we
specify second order normal derivatives due to outflow boundary conditions at part of the boundary,
the set of rows that can be manipulated is restricted to those corresponding to cells right next to
these outflow boundaries. In some extreme cases with outflow boundaries, the rank deficiency of A
can even be more than 1, i.e. the dimension of the null space is higher than 1. We will discuss the
issue of proper outflow boundary conditions in a separate section.
Now, let’s come back to the case where A is singular with a one-dimensional null space spanned
by e = [1, . . . , 1]T and b lies in the image of A. We first find a new linear equation that can nail
down the arbitrary uniform constant in the solution and then just replace the chosen row in the
augmented matrix [A, b] with this new equation or add the new equation to the original row. In
principle, any new equation is valid as long as its coefficients do not sum to 0. We usually pick a new
equation that is simple and interpretable. We also want the equation to have few nonzero entries
and the nonzero entries to be close to the diagonal, so that we preserve as much as possible the
banded structure and properties such as symmetry and diagonal dominance of the original matrix.
For example, if we want to manipulate the i-th row, a good new equation would be αxi = 0, where
α is an arbitrary nonzero coefficients whose value can be picked according to our need. The idea
of setting the spatial average to 0 by x1 + · · · + xN = 0 will introduce a whole row of 1’s and
significantly change the matrix’s structure. Moreover, regarding the type of manipulation, adding
seems to be a better strategy than replacing because adding the equation αxi = 0 is equivalent to
simply adding a constant α to the i-th diagonal entry, which does not change at all the pattern of
the nonzero entries and the symmetry.
A second idea comes from Seibold (2008), whose general form only works for symmetric A but with
one choice of parameters, it reduces to adding a positive number to any diagonal element, which
somewhat coincides with the first idea. This second idea is based on the fact that since A is rank
deficient by 1, we can make A nonsingular by adding a certain rank-1 matrix wwT to it. Note that
this addition preserves the symmetry. Of course, the choice w is not arbitrary because it needs to
satisfy two constraints. First, A0 = A + wwT must be nonsingular. Second, the unique solution to
A0 x = b should also be a solution to Ax = b.
15
Here we will show that if wT e 6= 0 (e = [1, . . . , 1]T ), the above two constraints will be satisfied.
The proof will use the property of a symmetric matrix A that Ker(A) ⊥ Im(A). This can be
proved as below. If x ∈ Ker(A), we have Ax = 0. If y ∈ Im(A), there exists z such that y = Az.
Therefore,
xT y = xT (Az) = (xT A)z = (AT x)T z = (Ax)T z = 0,
and thus Ker(A) ⊥ Im(A). Note that the last but one equality uses the symmetry A = AT .
Now we want to prove A0 = A + wwT is nonsingular. We just need to show if A0 x = 0, x = 0 has
to be true. Since A0 x = Ax + wwT x = Ax + (wT x)w and eT Ax = 0 (e ∈ Ker(A), Ax ∈ Im(A)),
if A0 x = 0, we have 0 = eT A0 x = (wT x)(eT w). Since we pick w such that eT w 6= 0, we must
have wT x = 0. Therefore, A0 x = 0 implies Ax = 0, which means x ∈ Ker(A) = span{e}. Hence,
x = αe for some real number α. However, since eT w 6= 0, for wT x = αwT e = 0 to be true, the
only possibility is α = 0 and thus x = 0. This proves that A0 is nonsingular.
It remains to show that assuming b ∈ Im(A), A0 x = b implies Ax = b. If A0 x = b, we have
Ax = AA0−1 b and we want to show AA0−1 b = b. Since AA0−1 = (A0 −wwT )A0−1 = I −wwT A0−1 =
I − w(A0−1 w)T , it remains to show that w(A0−1 w)T b = 0, which is true if and only if (A0−1 w)T b =
0. Note that A0 e = Ae + (wT e)w = (wT e)w. Therefore, A0−1 w = (wT e)−1 e. On the other hand,
since b ∈ Im(A) and e ∈ Ker(A), we have eT b = 0. Hence, (A0−1 w)T b = (wT e)−1 eT b = 0, which
completes the proof.
This idea gives us plenty of flexibility because the only requirement on w is wT e 6= 0, i.e. the
orthogonal projection of w in Ker(A) is nonzero. However, if we want to preserve the sparsity of A,
w must only have few nonzero entries. A simplest choice is w = [0, . . . , 0, α, 0, . . . , 0]T , whose i-th
entry is the only nonzero entry and equals to α. In this case, A + wwT is equivalent to adding a
positive number α2 to the i-th diagonal element of A. This coincides with the first idea except for
a subtlety. If A is symmetric, using the argument of the second idea, we can add a positive number
to any of the diagonal element. But if A is asymmetric, the argument of the second idea does not
apply and we have seen that in some cases, such as the one with the second order normal derivative
specified at part of the boundary, the choice of the row to be manipulated can be restricted.
There is another perspective from which we can make sense of the second idea. 2 We know that
solving a symmetric positive definite linear system Ax = b is equivalent to solving the optimization
problem
1 T T
arg min x Ax − b x .
x∈RN 2
Now since A has a one-dimensional null space span{e}, the minimizer is not unique because any
two vectors whose difference lies in span{e} will share the same objective function value (Note
that bT e = 0). To nail down only one minimizer, we need an extra term in the objective function
to distinguish two vectors that are different by a proportion of e. A straightforward idea is to
find a vector w such that wT e 6= 0, so that the term wT x can produce the distinction. To make
the objective function bounded from below to guarantee the existence of a minimizer, we add a
quadratic penalization term 1/2 · (wT x)2 = 1/2 · xT wwT x to the objective function:
1 T 1 1
x Ax + xT wwT x − bT x = xT (A + wwT )x − bT x
2 2 2
and use its hopefully unique minimizer as an approximation to a solution to Ax = b. We say
“hopefully” and “approximation” because at first sight, we have no idea whether this penalization
2
Thanks to Florian Feppon for bringing this up.
16
term can nail down only one minimizer and whether it will shift the minimizer outside Ker(A).
However, note that minimizing this objective function is equivalent to solving (A + wwT ) = b. By
the previous argument for the second idea, we know that indeed the minimizer is unique and it is
a solution to Ax = b. This perspective from the optimization formulation of the problem gives us
further intuition and motivation for the second idea.
For both the new equation in the first idea and the vector w in the second, we have the freedom to
choose a nonzero proportion coefficient. Analytically this choice does not matter, but numerically
it will affect the condition number of the manipulated matrix A0 . When this coefficient is either
too big or too small, it will yield a ill-conditioned A0 . Typically all the nonzero entries of A are
of the same order of magnitude, so to minimize the condition number of A0 , we should choose the
proportion coefficient such that the numbers added to A are of the same order of magnitude as the
entries of A.
In the above two ideas, we assume that b ∈ Im(A). What R if this isH not the case? This brings us
back to the question whether the compatibility condition Ω f dV = ∂Ω gdS naturally gives rise to
b ∈ Im(A) for any specific discretization scheme. The answer is in general “no”, which can be easily
verified empirically. This is also intuitive because no matter which type of discretization is used,
whether finite difference, finite volume,
R or finite element, there is no mechanism
H that tries to match
the errors in numerically computing Ω f dV with those in computing ∂Ω gdS. This mismatch can
be reduced when the mesh is refined and the truncation errors diminish but its presence can never
be eliminated. Besides, round-off errors can also contribute to this mismatch. This mismatch of
numerical errors is what makes b ∈ / Im(A).
Then what should we do? If we assume A is symmetric, which is usually the case for most typical
discretization schemes, a straightforward fix is to manually project b onto Im(A) to enforce b ∈
Im(A), since the projection of b onto Im(A)⊥ = span{e} is totally due to numerical errors and
is uninformative and thus can be tossed away. Another motivation of this manipulation is that
this is exactly what we would do in solving a least squares problem. Since e = [1, . . . , 1]T , this
manipulation is nothing but removing the mean of the entries of b from each of b’s entry.
On the other hand, however, we know that after making A into A0 , which is invertible, the linear
system A0 x = b is indeed well-posed. Then why can’t we just solve for x̂ = A0−1 b? Is x̂ any
different from x = A0−1 bIm(A) , where bIm(A) is the projection of b onto Im(A)? If we denote the
projection of b onto span{e} by be , we can write b = bIm(A) + be . Therefore, x̂ = A0−1 b =
A0−1 bIm(A) + A0−1 be = x + A0−1 be . We know that x̂ is an equivalent solution to x if and only if
their difference A0−1 be is some proportion of e. Since be ∈ span{e} and be 6= 0, this is equivalent
to saying that e is an eigenvector of A0 . However, A0 e = Ae + wwT e = (wT e)w and hence A0 e k e
if and only if w k e, which is never the case in practice because this choice of w will ruin the
sparsity of matrix A. Therefore, in general, Ax̂ 6= bIm(A) . Indeed, we have Ax̂ = Ax + AA0−1 be =
bIm(A) + (A0 − wwT )A0−1 be = bIm(A) + be − w(wT A0−1 )be = b − (wT e)−1 (eT be )w. Since wT e 6= 0
and eT be 6= 0, we can see that Ax̂ is neither bIm(A) nor b, which makes solving for x̂ = A0−1 b not
a good idea.
There is yet another subtlety on how different is the x obtained here from the one, denoted by
x+ , obtained by solving a least squares problem, or equivalently, by taking a pseudo-inverse. If we
compute x+ = A+ b, where A+ is the classic Moore-Penrose pseudo-inverse of A, what we are doing
is essentially first projecting b onto Im(A) to ensure the existence of a solution, and then picking
the minimum-2-norm solution if the solution is not unique. The first step is the same, as has been
17
mentioned before. The difference here is that x is in general not the minimum-2-norm solution.
The minimum-2-norm solution should be orthogonal to the null space span{e}, i.e. orthogonal
to e. However, eT x = eT A0−1 b, assuming b ∈ Im(A), and A0−1 e ∦ e unless w k e. Therefore,
generically eT x 6= 0. Actually, we have wT x = wT A0−1 b = 0 because A0 e k w, so x is ensured
to be orthogonal to w rather than e. However, if we want, we can easily remove the entry mean
from x to obtain the minimum-2-norm solution x+ . Therefore, we can view the idea introduced
here as another way to apply a pseudo-inverse. This new way bypasses the formation of orthogonal
matrices and thus allows us to make use of the sparsity of the linear system to accelerate the
computation and reduce the storage requirement, but the cost to pay is that by not resorting to
orthogonality-based algorithms, the conditioning is worsened in computation and the numerical
solution is susceptible to larger errors.
3 Download and Installation
3.1 Downloading
The code can be downloaded from http://mseas.mit.edu/codes. You will have to register to get
access to the site. Once you are logged in, you may download the latest Framework from the “Down-
loads::Dynamic Models” link: http://mseas.mit.edu/codes/downloads/dynamic-models. The
code is also directly available on Stellar for 2.29 students.
3.2 Installing on Linux
Installation is simple on a Linux machine. Simply untar the file in the directory of your choice, for
example:
% cd $install-root\
% tar xvfz FV_MATLAB_Framework_vX-XXX.tar.gz
where $install-root is a directory of your choice, and X-XXX indicates the version number.
Some files are coded in C++ for efficiency and these need to be compiled. If you are lucky, the files
that come with the framework that we compiled on our machine will work on your machine. But
if this is not the case, you will need to ”mex” the routines yourself.
On Linux this is a simple task.
% cd $install-root\
% matlab -nodesktop -nosplash
>> MexSubroutines
and you are ready to go!
If this does not work, it is possible that you do not have a C++ compiler installed on your system.
In that case, you may have to install a compiler and set it up in MATLAB. For example, if you are
on an Ubuntu system, the following is a start:
% sudo apt-get install build-essential
% matlab -nosplash -nodesktop
18
>> mex -setup
at which point you are on your own, good luck!
3.3 Installing on Mac
Installation is a little more complicated a Mac. Start by uncompressing the file in the directory of
your choice.
Some files are coded in C++ for efficiency and they have to be compiled. If you are lucky, the files
that come with the framework that we compiled on our machine will work on your machine. You
can first try things out, but if our compiled routines don’t work on your machine, you will need to
”mex” the routines yourself.
On a mac, you first have to install Xcode, including the Developer tools. Once you have these
installed, you will have a ”‘/Developer/usr”’ directory on your Mac. Once properly installed, you
should be able to mex the required routines. Start by opening MATLAB, and from the command
window type:
1 >> MexSubroutines
and you are ready to go!

If this does not work, it is possible that the /Developer directory is not in you path. To check,
type
>> !echo $PATH
From MATLAB. If you do not see /Developer/usr/bin and /Developer/usr/lib/, then MATLAB
will not be able to find the compilers. To fix this, open a new X11 terminal. Start by checking if
the directories are in the path, if not add, them, then open MATLAB from the terminal:
bash-3.2$ !echo $PATH
bash-3.2$ export PATH=$PATH:/Developer/usr/bin
bash-3.2$ export PATH=$PATH:/Developer/usr/lib
bash-3.2$ cd /Applications/MATLAB_R2018a/bin
bash-3.2$ ./matlab
Then change folders until you get to the trunk/Src/Advect folder, and issue the commands above.
3.4 Installing on Windows
Installation on Windows is a little more complex than on a Linux machine. First the downloaded
file needs to be untarred by a program such as IZArc. Untar the file to a directory of your choice
$install-root (for example, $install-root=C:\Users\usrname\Documents\MIT229FV\).
Some files are coded in C++ for efficiency and they have to be compiled in a special way to be
used in MATLAB, which is called to be mexed. If you are lucky, the mexed files (.mexa or .mexw,
depending on operating systems) coming with the framework will work on your machine. But if this
19
is not the case, you will need to mex the routines yourself. You may be able to use the C-compiler
coming with MATLAB. To use it, start with:
1 >> mex −setup
Then you may compile the .cpp subroutines as follows:
1 >> MexSubroutines
If this does not work, you may need to install a compatible compiler. You can find relevant
information and download links for compatible compilers in different operating systems on the web
page https://www.mathworks.com/support/compilers.html. After doing this, you can mex the
necessary files as above. At this point, you should hopefully be ready to go!
4 Quick Start Guide
After starting MATLAB, browse to the root directory (directory containing Example.m).
To see if the Framework is correctly installed with all components working, you can run the 13
examples we have provided by issuing the following commands. From the root directory (directory
containing setup.m) run:
1 >> Example(?);
where ? should be replaced by any integer from 1-13.

The Example.m is basically a high-level interface to all the solvers and examples that we provide.
For detailed information about the inputs, outputs, different examples and solvers, type
1 >> help Example
at the MATLAB prompt. The examples in TestCase have been tested thoroughly, and should
work without error when using all default parameters.
Example may be invoked with 1, 2 inputs. By running Example.m with only 1 input, default values
are used for the solver, resolution, diffusivity, and other parameters. The second input is a structure
with which you may adjust many of the other solution parameters. If you adjust these values and
the solution no longer works, please see the Troubleshooting section of this document.
For example, you may want to run the Lock-Exchange example with the default solver, but a
different ν. This can be done as follows:
1 >> app.nu = 0.1

2 >> Example(6,app);
20
4.1 Brief Description of Provided Examples
1. Heated plate demo: This test-case demonstrates the different boundary conditions and geom-
etry that can be used on the different staggered grids (see node). Its main purpose, initially,
was to validate that the diffusion operator was working correctly.
2. Tracer advection demo with QUICK scheme: This example is used for a convergence study
to verify the correctness of the advection scheme. This demo uses the Quadratic Upwind
Interpolation for Convective Kinetics (QUICK) scheme (see QUICK). Its performance is
contrasted with the upwind scheme in the next example.
3. Tracer advection demo with Upwind scheme: This example is included to highlight the poor
performance of the low-order upwind scheme (see UW), which is still commonly used in many
CFD codes, compared to the previous example using the QUICK scheme.
4. Analytical Stokes Example: This example is used to verify the correctness of the Stokes solver
implementations. It is described in detail in Stokes.
5. Poiseuille flow test case: This example solves the classical Poiseuille flow in a pipe. It tests
the implementation of the open boundary condition.
6. Lock exchange experiment: This example tests Boussinesq solvers, that is, density driven
flows. The setup is as follows: initially a dense and light fluid are separated by a barrier or a
“lock.” At time 0, the barrier is instantaneously removed, and the flow evolves.
7. Sudden expansion: This example demonstrates the accuracy of a control-volume analysis com-
pared to a Bernoulli analysis for predicting the pressure drop across a sudden expansion in a
pipe. Students in MIT course 2.29 solve this problem analytically using the two approaches,
and verify the answer numerically at a later point in the course. The Bernoulli approach is
incorrect in this case, because the recirculation zones that develop violate the - no viscous
effects assumption - upon which Bernoulli is based. This example ensures that the numerical
solution remains symmetric up to a transition Reynolds number where numerical roundoff can
become sufficiently large to cause flow perturbations. It is also used to test the open boundary
condition, and it demonstrates that the domain masking works correctly. The setup script
for this example also shows how to create more complex geometries. This case can readily be
modified for the classical flow over a cylinder test-case. The setup file is described in detail
in expan.
8. Buoyant sudden expansion: This is nearly the same as the previous example, but a lighter
density fluid is allowed to enter the domain, causing fluid to rise upwards after the expansion.
This introduces an asymmetry into the problem, and vortex shedding is seen using the default
parameter values.
9. Warm rising bubble: This is another density driven flow example. Here, the solution is ini-
tialized with a small bubble of warm water surrounded with cold water in a tube with closed
ends. The bubble rises and diffuses/mixes with the cold water as the flow evolves. The plot
script for this example shows an elegant method for saving outputs (such as images of .mat
files). It involves adding additional fields to the app data-structure.
10. Deterministic Double-Gyre test case: This example makes use of the deterministic solver with
the Coriolis forces. The fluid is initially at rest, and at time 0, a wind instantaneously starts to
21
force the 2D horizontal flow. The 2D wind on top of the fluid is the shape of a cosine function,
with maximum amplitudes at the top, middle, and bottom of the 2D horizontal domain. The
direction of the wind at the top and bottom are to the left, whereas in the middle the direction
is to the right. This drives the flow to form a double-gyre, which becomes unstable after long
integration times for high enough Reynolds numbers. The default parameter values do no
reach the point of instability, since simulations that do take considerable time to complete.
11. Lid-driven cavity flow: This example is a standard test-case for CFD codes. Fluid initially at
rest is inside a closed box (or cavity). At the time 0, the lid (or top boundary) instantaneously
moves with a fixed velocity, and continues to do so until the end of the simulation. For low
Reynolds numbers, a steady state is reached, whereas for sufficiently high Reynolds numbers,
a steady state is not reached.
12. Lock-Exchange flow: In a lock exchange experiment, fluids of different densities initially at
rest are separated by a vertical barrier. When the fluids are allowed to mix, differences in the
hydrostatic pressure cause the denser fluid to flow in one direction along the bottom boundary
of the domain, while the lighter fluid flows in the opposite direction along the top boundary
of the domain.
13. Deterministic Rayleigh-Bernard Instability: This is a buoyancy-driven flow, where the top is
being cooled and the bottom heated – resulting in an unstable density stratification which
drives flow. This classic test-case can have multiple solutions depending on the input param-
eters (for example, clockwise versus counter-clockwise circulation).
5 Setting Up New Examples
This section is for users who wish to set up new problems using this framework or want to have a
better understanding of the inner working of this code. We remind the reader that this is not an
exhaustive manual, but should be sufficient for students who have been introduced to numerical
fluid mechanics. In general the MATLAB functions are well-commented, and for further study,
users are encouraged to read the code.
5.1 Naming Conventions
In the code two naming conventions are used. When dealing with a single control volume centered
at (x, y) we name the locations (x − ∆x/2, y),(x + ∆x/2, y),(x, y − ∆y/2), and (x, y + ∆y/2)
as [west/left], [east/right], [south/bottom], [north/top], respectively, using LOWERCASE letters.
Similarly, we name the locations (x − ∆x, y),(x + ∆x, y),(x, y − ∆y), and (x, y + ∆y) as [West/Left],
[East/Right], [South/Bottom], [North/Top], respectively, using UPPERCASE letters. Finally, it
is customary to label the value of the function centered at the present CV as ”Present” or just
”P”. We note that this may be confusing with the Pressure, however, Pressure is never used as
a subscript, so hopefully the meaning will be clear from the context. This naming convention is
illustrated in Fig. 1.
Usually a single letter of the alphabet is used for integer numbers used in for loops, and generally
the first for loop is started using the letter ”i”. Also, generally ”i” is used for the rows of matrices,
and ”j” is used for the columns of matrices.
22
Figure 1: Naming conventions
5.2 Node Numbering Convention
The tracer, u-velocity and v-velocity grids are shown superimposed in Fig. 2 and separately in
Fig. 3. The numbering convention used in the code is illustrated in Fig. 2. A node matrix will
be used and it will be in the form Node(i,j). Therefore, the ”i” index will be used for the y-
direction, and the ”j” index will be used for the x-direction. This convention was chosen such
that the way a matrix is printed on screen in MATLAB mirrors the numbering convention of the
physical domain. It is nearly ideal, since increasing the column number of the matrix j increases
the x spatial direction. However, to increase the y spatial direction, the row number of the matrix
i has to decrease. Note, nowhere are the node-numbering matrices required, but these are used for
coding convenience, and (hopefully) clarity.
5.3 Data Structures
Since structured grids are used, the data-structures are reasonably simple. Unknowns are stored
in 1-dimensional arrays. The first N bcs entries are reserved for boundary conditions, and are used
to specify the value of boundaries. Time varying boundary conditions are allowed in this way, and
coded through the updateBcs function.
The application data-structure app is used as the general input structure. By using a structure,
users are allowed to pass in whatever information they desire to be used in the setup scripts. For ex-
ample, one could include an app.bcsDu field that controls the u-velocity boundary condition.
The remaining data-structures are restricted to the matrix operators. All matrix operators are
stored as sparse two-dimensional arrays.
23
Figure 2: Grid setup
24
Figure 3: Grid setup
25
5.4 Domain, Masking and Boundary Conditions: General Rules and Guide-
lines
Here we describe the overall rules and guidelines required in the set-up of the domain, masking and
boundary conditions for all simulations and test cases. Some users may want to check sub-sections
on a forced Stokes flow (sub-section 5.5) and on a Sudden Expansion flows (sub-section 5.6) for
more detailed examples. In what follows, the description is a bit more generic.
5.4.1 Domain
The domain is set-up as a rectangle. The call to mesh the domain is:
1 [XP, YP, XU, YU, XV, YV, dx, dy] = meshdomain(app, [0, a], [0, b])
We note that domains are always rectangles/squares, so up to the number of finite volumes, the
set-up of the domain is pretty much the same for all examples. What varies from one example to
the next are the masking and boundary/initial conditions, and of course the equations used.
5.4.2 Masking
This mask is user-created. Zeros are unmasked, and nonzero integers indicate a masked region.
The four domain sides are not considered as masked, but they have a separate number assigned to
them in the Node numbering arrays.
To allow for different boundary conditions at the edges of each masked regions, the interior of each
masked region must have a different integer number. If two masked region have the same number,
they will be considered as a single masked region with its own definition of boundary condition at
the edges of the mask.
First, the whole mask is created that fills the whole rectangular domain. Then the user can define
additional masked regions using the position arrays XP and YP, and Boolean rules.
For example, let’s assume we would like to define a domain for sudden expansion, but with a
elliptical obstacle (Fig. 4). We then need to define three masked region, two rectangles and one
ellipse. Each one of these are defined using their positions and Boolean rules.
For the two rectangles, we have:
1 top rect = @(x, y) logical((y \ge 2 / 3) .* (x \le 4));

2 bot rect = @(x, y) logical((y \le 1 / 3) .* (x \le 4));
For the ellipse, we have:
1 ellipse = @(x, y) (((x − 10) .ˆ2 / 20 + (y − 0.5) .ˆ 2) < 0.1);
We then assign numbers to the three masked region, for example, associate 1 with the top domain,
2 for the bottom domain and 3 for the circle.
26
Figure 4: Sudden Expansion Boundary ID setup with an additional cirle region masked.
27
1 Mask(top rect(XP, YP)) = 1;
2 Mask(bot rect(XP, YP)) = 2;
3 Mask(ellipse(XP, YP)) = 3;
5.4.3 Boundary Conditions
The type of boundary conditions needs to be specified. The available types are:
1. ’D’: Dirichlet with a single value
2. ’Dmv’: Dirichlet with multiple values (i.e. spatially variable)
3. ’O’: Open boundary
4. ’N’: Neumann with a single value
5. ’Nmv’: Neumann with multiple values (i.e. spatially variable)
The boundary type is specified in a cell array bcid2type{i}, where i corresponds to the boundary
condition id number. For example, the ellipse has boundary id=3 (Fig. 4), so to set a Neu-
mann boundary on u-velocity for the ellipse, you would ensure that bcid2type u{e}=’N’. Each of
pressure, u-velocity, v-velocity, and optionally density should have its own array. The boundary
condition type is set using the set bcs function.
Initially, all BC values are set to zero everywhere. Some of these boundary conditions may perhaps
not be used in a specific example (for example, if density or temperature is not used), but it is still
fine to initialize them.
1 [bcsDP, bcsOP, bcsNP] = init bcs(nbcP);

2 [bcsDu, bcsOu, bcsNu] = init bcs(nbcu);
3 [bcsDv, bcsOv, bcsNv] = init bcs(nbcv);
4 [bcsDrho, bcsOrho, bcsNrho] = init bcs(nbcrho);
Hence, the above set-to-zero lines should in general not be changed in the set-up part of the
code.
Next, the user can modify these zero boundary conditions to the values that they wish to use. First
we obtain the coordinates for the boundary points. Each of the XY arrays contain two columns, the
first gives the x-coordinate, the second gives the y-coordinate associated with a particular boundary
(D, O, or N).
The four sides of the domains are set-up directly, without the need of using boolean selections
based on position arrays XP and YP, or XU and YU. However, for all masked regions, first, we
re-select the region of the masks that needs to be assigned to a value, this is done as in the mask
definition themselves. Then, a value is assigned to this. For example, to set the top half of the inlet
boundary condition to 1, and the bottom half to 0.5, you would first ensure that bcid2type rho{4}
= ’Dmv’, then:
1 %Get the boundary coordinates

2 [XYD, XYO, XYN] = get bc coord('u', XU, YU, Nodeu, bcsDu, bcsOu, bcsNu);
28
3 ids = (XYD(:, 1) == 0);
4 bcsDu(ids & XYD(:, 2) > 0.5) = 1;
5 bcsDu(ids & XYD(:, 2) ≤ 0.5) = 0.5;
For the masked regions in the interior, you can also add boundary values that vary in space. For
example, to add a Dirichlet parabolic density profile to the ellipse, you would first ensure that
bcid2type u{3} = ’Dmv’, then:
1 [XYD, XYO, XYN] = get bc coord('P', XP, YP, NodeRho, bcsDrho, bcsOrho, bcsNrho);
2 ids = ellipse(XYD(:, 1), XYD(:, 2));
3 bcsDrho(ids) = (XYD(ids, 1) − mean(XYD(ids, 1))).ˆ2;
We note that the script set bcs ensures that all conventions needed for the code are followed. In
most cases, the user, thus, does not need to worry about this: she/he can use his/her own mask
numbers and set bcs will reset them (it is only in some rare cases that it may need some updates).
The conventions that the arrays NodeP, Nodeu, Nodev, and NodeRho, should follow once they come
out of the script set bcs are:
1. No Neumann boundary condition can have a number smaller than a Dirichlet or Open bound-
ary condition. That is we require BCidNeumann > BCidDirichlet .
2. No Open boundary condition can have a number smaller than a Dirichlet condition, or larger
than a Neumann condition. That is we require BCidNeumann > BCidOpen > BCidDirichlet .
3. No Dirichlet boundary condition can have a number higher than an Open or Neumann bound-
ary condition. That is we require BCidNeumann > BCidOpen > BCidDirichlet .
4. The number of boundary conditions for velocity, pressure, and density have to be the same.
5. The bcsD arrays need to include open boundary conditions, and the bcsO arrays need to
indicate which of the Dirichlet boundaries are actually open boundaries.
6. Time-varying boundary conditions are implemented through the user-provided updateBcs
function. If this function is not provided, boundary conditions will not be time-dependent.
The function needs to have the following outputs, and take the following inputs:
1 [bcsDP, bcsNP, bcsDrho, bcsNrho, bcsDu, bcsNu,bcsDv, bcsNv]...

2 = updateBcs(time, XP, YP, XU, YU, XV, YV, app, P, rho, u, v,...
3 NodeP, NodeRho, Nodeu, Nodev);
The rest of this section walks through how to create the SuddenExpansion and StokesForcing input
files, including also the initial conditions.
29
5.5 Stokes Forcing Test Case:
This is an analytical test case used for testing the convergence of various Stokes solvers. The
solution is:
u = π sin(t) sin(2πy) sin2 (πx) (42)

2
v = −π sin(t) sin(2πx) sin (πy) (43)
P = sin(t) cos(πx) sin(πy) (44)
and is solved on a square domain of size [−1, 1] × [−1, 1].

The solver function requires 3 inputs: app,SetupScript,PlotScript. The first input is a structure,
where the fields app.Nx, app.Ny, app.T and app.dt define the discretization, giving the number
of INTERIOR cells in the x-direction, y-direction, the duration of the simulation, and the time step
size, respectively. The value of the kinematic viscosity (which is analogous to the inverse Reynolds
number for the full NS equations) is set through the app.nu field. The app.PlotIntrvl field is
an integer number that indicates the number of time integrations steps between calling the script
given by the string input PlotScript. The second input is the main focus of the quick-start guide,
and it is a string giving the name of a MATLAB script which sets up the necessary variables for
the problem to be solved. That is, the solver calls eval(SetupScript) to setup the problem at
hand. Similarly, the solver calls eval(PlotScript) every app.PlotIntrvl steps of the integration.
While PlotScript is normally used for plotting the solution, it can also be used for saving the
solution, and various other user-defined tasks. Also, any uncleared variables in SetupScript will
be available for use in PlotScript. Therefore, a function handle created in SetupScript can be
used in PlotScript without causing an error.
The setup file is in its entirety as follows:
1 % Stokes Flow forcing function test case

2
3 %% Forcing Function
4 app.nu = 1;
5 Fu = @(time,p) ...
6 pi.*cos(time+1).*sin(2.*pi.*p(:,2)).*sin(pi.*p(:,1)).ˆ2−...
7 2.*pi.ˆ3.*sin(time+1).*sin(2.*pi.*p(:,2)).*cos(pi.*p(:,1)).ˆ2+...
8 6.*pi.ˆ3.*sin(time+1).*sin(2.*pi.*p(:,2)).*sin(pi.*p(:,1)).ˆ2−...
9 sin(time+1).*sin(pi.*p(:,1)).*pi.*sin(pi.*p(:,2));
10 Fv = @(time,p) ...
11 −pi.*cos(time+1).*sin(2.*pi.*p(:,1)).*sin(pi.*p(:,2)).ˆ2−...
12 6.*pi.ˆ3.*sin(time+1).*sin(2.*pi.*p(:,1)).*sin(pi.*p(:,2)).ˆ2+...
14 sin(time+1).*cos(pi.*p(:,1)).*cos(pi.*p(:,2)).*pi;
15
16 %% Exact Solution (used in plotting function)
17 uexact = @(time,x) pi*sin(time+1).*(sin(2*pi*x(:,2)).*(sin(pi*x(:,1))).ˆ2);
18 vexact = @(time,x) −pi*sin(time+1).*(sin(2*pi*x(:,1)).*(sin(pi*x(:,2))).ˆ2);
19 pexact = @(time,x) sin(time+1).* cos( pi*x(:,1)).* sin(pi*x(:,2)) ;
20
21 %% Set domain size
22 a = 2; b = 2;
23 dx = a./(app.Nx); dy = b./(app.Ny);
24
25 %% Set boundary conditions
30
26 bcsDP= []; bcsNP= [0 0 0 0];
27 bcsDu= [0 0 0 0]; bcsNu= [];
28 bcsDv= [0 0 0 0]; bcsNv= [];
29
30 %% Set number of timesteps
31 Nt = app.T/app.dt;
32
33 %% create geometry mask
34 Nbcs = length(bcsDP) + length(bcsNP);
35 Mask = zeros(app.Ny, app.Nx);
36 %create Node matrices
37 [NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,[1 1 1 1],1);
38
39
40 %% Initialize fields
41 %create vectors of coordinates for p, u, and v control volumes
42 xp = linspace(dx/2/a−0.5, 0.5−dx/2/a, app.Nx)*a;
43 yp = fliplr(linspace(dy/2/b−0.5, 0.5−dy/2/b, app.Ny))*b;
44 xu = linspace(dx/a−0.5, 0.5−dx/a, app.Nx−1)*a;
45 yu = fliplr(linspace(dy/2/b−0.5, 0.5−dy/2/b, app.Ny))*b;
46 xv = linspace(dx/2/a−0.5, 0.5−dx/2/a, app.Nx)*a;
47 yv = fliplr(linspace(dy/b−0.5, 0.5−dy/b, app.Ny−1))*b;
48
49 %create matrix of x and y coordinates
50 [XP YP] = meshgrid(xp,yp);
51 [XU YU] = meshgrid(xu,yu);
52 [XV YV] = meshgrid(xv,yv);
53
54 %Initialize vectors
55 P = [bcsDP; bcsNP; pexact(0,[XP(idsP) YP(idsP)])];
56 u = [bcsDu; bcsNu; uexact(0,[XU(idsu) YU(idsu)])];
57 v = [bcsDv; bcsNv; vexact(0,[XV(idsv) YV(idsv)])];
58 Pex = [bcsDP; bcsNP; pexact(0,[XP(idsP) YP(idsP)])];
59 uex = [bcsDu; bcsNu; uexact(0,[XU(idsu) YU(idsu)])];
60 vex = [bcsDv; bcsNv; vexact(0,[XV(idsv) YV(idsv)])];
5.5.1 Setting Up Forcing Functions
The forcing functions can take time and space as the input variables. The spatial variable is
expected to be a matrix with two columns. The first column should contain the x locations of
control volume (CV) centres, and the second column should contain the y locations of CV centres.
To set up the forcing functions corresponding to the solutions (44) we proceed as follows:
1 Fu = @(time,p) ...
2 pi.*cos(time+1).*sin(2.*pi.*p(:,2)).*sin(pi.*p(:,1)).ˆ2−...
4 6.*pi.ˆ3.*sin(time+1).*sin(2.*pi.*p(:,2)).*sin(pi.*p(:,1)).ˆ2−...
5 sin(time+1).*sin(pi.*p(:,1)).*pi.*sin(pi.*p(:,2));
6 Fv = @(time,p) ...
7 −pi.*cos(time+1).*sin(2.*pi.*p(:,1)).*sin(pi.*p(:,2)).ˆ2−...
8 6.*pi.ˆ3.*sin(time+1).*sin(2.*pi.*p(:,1)).*sin(pi.*p(:,2)).ˆ2+...
10 sin(time+1).*cos(pi.*p(:,1)).*cos(pi.*p(:,2)).*pi;
31
You will note we used “time+1” in the forcing functions, and this is only to start the simulation
with a nonzero forcing.
5.5.2 Setting Up the Grid
The first step in creating the grid is setting the constant discretization sizes ∆x, ∆y, and number
of time-integration steps N t, which is accomplished simply as:
1 %% Set domain size

2 a = 2; b = 2;
4
5 %% Set number of timesteps
6 Nt = app.T/app.dt;
since here the domain is 2 × 2 units long, and we are integrating for T /∆t timesteps.
The Stokes solvers require 3 node numbering matrices, NodeP, Nodeu, and Nodev for the Pressure,
u-velocity, and v-velocity. A simple way to create these matrices is by providing a “masking”
matrix to the NodePad function. Since there is no interior geometry for this test-case (that is no
obstructions or holes in the interior of the domain), we can provide a mask with all-zero values the
size of the domain to the NodePad function as follows:
1 Mask = zeros(Ny, Nx);

3 [NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,[1 1 1 1],1);
In this case there are N x interior control volumes (or nodes) in the x-direction, and N y interior
control volumes in the y-direction. The mask is only created for the interior control volume.
Therefore NodePad “pads” the masking matrix to add the four boundaries conditions on each edge.
The second input to NodePad is a logical array, and tells the function which of the boundaries to
pad ([left, right, bottom, top]), and the last input is a logical scalar which checks for periodicity,
automatically adding periodic boundaries if required.
Note that the ID of the first interior degree of freedom is N bcs + 1, where N bcs is a scalar variable
containing the number of boundary conditions, in this case N bcs = 4. Therefore the first interior
node ID is 5, and this will be the first unknown that will be solved (that is, the first row in the A
matrix).
32
An example of NodeP, Nodeu, and Nodev for a N x = 5, N y = 4 domain is as follows:
 
1 4 4 4 4 4 2
 1 5 9 13 17 21 2 
 
 1 6 10 14 18 22 2 
N odeP =   1

 7 11 15 19 23 2  
 1 8 12 16 20 24 2 
1 3 3 3 3 3 2
 
1 4 4 4 4 2
 1 5 9 13 17 2 
 
 1 6 10 14 18 2 
N odeu =   1

 7 11 15 19 2 

 1 8 12 16 20 2 
1 3 3 3 3 2
 
1 4 4 4 4 4 2
 1 5 8 11 14 17 2 
 
N odev =   1 6 9 12 15 18 2 

 1 7 10 13 16 19 2 
1 3 3 3 3 3 2
Significant care must be taken with properly setting up the

boundary conditions
There are a number of conventions related to the boundary conditions (see section rules). First,
let us define the variables used for the boundary conditions. In each case, these variables are
one-dimensional double arrays:
1 bcsDP: Pressure Dirichlet boundary conditions

2 bcsNP: Pressure Neumann boundary conditions
3 bcsOP: ID of Pressure Dirichlet boundary that should be an open boundary
4 bcsDu: U−velocity Dirichlet boundary conditions
5 bcsNu: U−velocity Neumann boundary conditions
6 bcsOu: ID of U−velocity Dirichlet boundary that should be an open boundary
7 bcsDv: V−velocity Dirichlet boundary conditions
8 bcsNv: V−velocity Neumann boundary conditions
9 bcsOv: ID of V−velocity Dirichlet boundary that should be an open boundary
The numbering convention comes from how the node numbering matrices are used. The first N bcs
numbers are reserved for the boundary conditions. So, in this case, the first 4 numbers are reserved
for boundaries. bcsDP, then, is an array that contains the VALUE of the Dirichlet boundary
conditions. Since this problem contains only uniform Dirichlet boundaries for velocity and uniform
Neumann conditions for Pressure, bcsDu, bcsDv, and bcsNP are all 1 × 4 arrays with zero entries.
The remaining arrays are empty. This is implemented as follows:
1 bcsDP= []; bcsNP= [0 0 0 0];

2 bcsDu= [0 0 0 0]; bcsNu= [];
3 bcsDv= [0 0 0 0]; bcsNv= [];
33
5.5.4 Initial Conditions
To create the initial conditions, first “coordinate” matrices, which contain the coordinates of the
control volume centres, are created as follows:
2 %create vectors of coordinates for p, u, and v control volumes
3 xp = linspace(dx/2/a−0.5, 0.5−dx/2/a, app.Nx )*a;
4 yp = fliplr(linspace(dy/2/b−0.5, 0.5−dy/2/b, app.Ny ))*b;
5 xu = linspace( dx/a−0.5 ,0.5−dx/a , app.Nx−1)*a;
6 yu = fliplr(linspace(dy/2/b−0.5, 0.5−dy/2/b, app.Ny ))*b;
7 xv = linspace(dx/2/a−0.5, 0.5−dx/2/a, app.Nx )*a;
8 yv = fliplr(linspace(dy/b−0.5 , 0.5−dy/b, app.Ny−1))*b;
9
10 %create matrix of x and y coordinates
11 [XP YP] = meshgrid(xp,yp);
12 [XU YU] = meshgrid(xu,yu);
13 [XV YV] = meshgrid(xv,yv);
note here we are using the MATLAB meshgrid function to create the matrices, and that these
matrices are the size of the INTERIOR nodes only, that is, they do not go all the way up to the
boundaries. Also note, the vertical coordinate starts from the maximum value
and decreases, and this is due to the node numbering convention described in Section node.
When initializing the variables, the boundary condition values are also stored in the vector of un-
knowns. Therefore, when initializing the vector of unknowns, the values of the boundary conditions
have to be included as well.
For the general case there may be interior masked nodes, and then the created coordinate matrices
will have too many entries, and will not necessarily correspond to the IDs of the unknowns. Here
the idsP, idsu, and idsv outputs from the NodePad function is useful for selecting the correct
elements of the coordinate matrices.
Initializing the vector of unknowns using the analytical solution to the problem at t = 0 is accom-
plished as follows:
1 %% Solution
2 uexact = @(time,x) pi*sin(time+1).*(sin(2*pi*x(:,2)).*(sin(pi*x(:,1))).ˆ2);
3 vexact = @(time,x) −pi*sin(time+1).*(sin(2*pi*x(:,1)).*(sin(pi*x(:,2))).ˆ2);
4 pexact = @(time,x) sin(time+1).* cos( pi*x(:,1)).* sin(pi*x(:,2)) ;
5 %Initialize vectors
6 P = [bcsDP; bcsNP; pexact(0,[XP(idsP) YP(idsP)])];
7 u = [bcsDu; bcsNu; uexact(0,[XU(idsu) YU(idsu)])];
8 v = [bcsDv; bcsNv; vexact(0,[XV(idsv) YV(idsv)])];
5.5.5 Misc. and Plotting
To plot the solution we can simply use surf(P(NodeP)) in the plotting script. To plot the solution
at the correct (x, y) coordinates for the interior we can use
surf(XP,YP,P(NodeP(2:end-1,2:end-1))).
34
This problem requires ν = 1, and to ensure this happens, we explicitly set ν = 1 in the SetupScript
to guarantee this happens.
We also created vectors Pex, uex, vex in the SetupScript, and these are again used in PlotScript
for plotting the exact solution.
5.6 Sudden Expansion
Note: This section is now a bit out-dated and could be updated for easier use with the new boundary
condition interface. However, it still works as is, so we leave in the manual for now.
This example demonstrates the accuracy of a control-volume analysis compared to a Bernoulli
analysis for predicting the pressure drop across a sudden expansion in a pipe. Students in MIT
course 2.29 solve this problem analytically using the two approaches, and verify the answer numer-
ically at a later point in the course. Bernoulli’s principle does not apply in this case, because the
recirculation zones violate the inviscid flow assumption, upon which Bernoulli’s principle is based.
This example ensures that the numerical solution remains symmetric up to a transition Reynolds
number where numerical roundoff is sufficiently large enough to cause flow perturbations. It is also
used to test the open boundary condition, and it demonstrates that the domain masking works
correctly. The setup script for this example also shows how to create more complex geometries.
This case can readily be modified for the classical flow over a cylinder test-case.
The Navier Stokes solvers require 3 inputs: app,SetupScript and PlotScript. The first input is
a structure, where the fields app.Nx, app.Ny, app.T and app.dt define the discretization, giving
the number of INTERIOR points in the x-direction, y-direction, the duration of the simulation,
and the size of the time step, respectively. The value of the kinematic viscosity (which is analogous
to the inverse Reynolds number) is set through the app.nu field, whereas the kinematic diffusivity
for the density is set through app.kappa. The app.PlotIntrvl field is an integer number that
indicates the number of time integrations steps between calling the script given by the string input
PlotScript. The second input is the main focus of this section, and it is a string giving the
name of a MATLAB script which sets up the necessary parts of the problem to be solved. That
is, the solver calls eval(SetupScript) to setup the problem at hand. Similarly, the solver calls
eval(PlotScript) every app.PlotIntrvl steps of the integration. While PlotScript is normally
used for plotting the solution, it can also be used for saving the solution. Also, any uncleared
variables in SetupScript will be available for use in PlotScript. Therefore, a function handle
created in SetupScript can be used in PlotScript without causing an error.
The setup file is in its entirety as follows:
1 % Sudden Expansion test case

2
3 %# of interior CV's
4 %Check of app.Ny is a multiple of 3
5 if mod(app.Ny,3)
6 fprintf('NOTE: app.Ny was not a multiple of 3, change app.Ny from %g to ...
%g\n',...
7 app.Ny,app.Ny−mod(app.Ny,3));
8 app.Ny = app.Ny−mod(app.Ny,3);
9 end
10 if app.Ny/3≤2
35
11 display('app.Ny needs to be at least 9 for this test case, app.Ny changed to ...
9.')
12 app.Ny = 9;
13 end
14
15 %Define Forcing functions
16 Fu = @(time,p) 0;
17 Fv = @(time,p) 0;
18
19 %Set size of domain
20 a = 20; b = 1;
22
23 %Set boundary conditions
24 bcsDrho = [0 0 0 0]; bcsNrho = [0];
25 bcsDP = [0]; bcsNP = [0 0 0 0];
26 bcsOP = [1];
27 bcsDu = [0 1 0 0 0]; bcsNu= [];
28 bcsOu = [5];
29 bcsDv = [0 0 0 0]; bcsNv= [0];
30
31 %Set number of timesteps
32 Nt = ceil(app.T/app.dt);
33
36 bnd = [1 1 1 1];
38 Mask(1:round(app.Ny/3),1:round(app.Nx/5))=1;
39 Mask(end−round(app.Ny/3−1):end,1:round(app.Nx/5))=1;
40
42 [NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,bnd,1);
43 NodeRho=NodeP;
44 NodeP(find(NodeP==1))=−1;NodeP(find(NodeP==3))=1;NodeP(find(NodeP==−1))=3;
45 Nodeu(find(Nodeu==3))=−1;Nodeu(find(Nodeu==5))=3;Nodeu(find(Nodeu==−1))=5;
46 Nodev(find(Nodev==3))=−1;Nodev(find(Nodev==5))=3;Nodev(find(Nodev==−1))=5;
47 NodeRho(find(NodeRho==3))=−1;NodeRho(find(NodeRho==5))=3;NodeRho(find(NodeRho==−1))=5;
48
50 xp = linspace(dx/2, a−dx/2, app.Nx );yp = fliplr(linspace(dy/2, b−dy/2, app.Ny ...
));
51 xu = linspace( dx ,a−dx , app.Nx−1);yu = fliplr(linspace(dy/2, b−dy/2, app.Ny ...
));
52 xv = linspace(dx/2, a−dx/2, app.Nx );yv = fliplr(linspace(dy , b−dy, ...
app.Ny−1));
53 [XP YP] = meshgrid(xp,yp);[XU YU] = meshgrid(xu,yu);[XV YV] = meshgrid(xv,yv);
54 rhoinit=@(X,Y) 0*((X>0.5) − (X<0.5))*0.5;
55 pinit = @(X,Y) zeros(size(X));
56 uinit = @(X,Y) zeros(size(X));
57 vinit = @(X,Y) zeros(size(X));
58 rho=[bcsDrho' ;bcsNrho' ;rhoinit(XP(idsP), YP(idsP))];
59 P = [bcsDP; bcsNP; pinit(XP(idsP), YP(idsP))];
60 u = [bcsDu; bcsNu; uinit(XU(idsu), YU(idsu))];
61 v = [bcsDv; bcsNv; vinit(XV(idsv), YV(idsv))];
36
5.6.1 Setting Up Forcing Functions
The forcing functions can take time and space as the input variables. For this test case there are
no forcing function contributions, so we simply set:
1 Fu = @(time,p) 0;
2 Fv = @(time,p) 0;
5.6.2 Setting Up the Grid
The first step in creating the grid is setting the constant discretization sizes ∆x, ∆y, and number
of time integration steps N t, which is accomplished simply as:
1 %Set size of domain

2 a = 20; b = 1;
4
5 %Set number of timesteps
6 Nt = ceil(app.T/app.dt);
since here the domain is 20 × 1 units long, and we are integrating for T /∆t time steps.
The Navier-Stokes solvers require 4 node numbering matrices, Noderho, NodeP, Nodeu, and Nodev for the density,
Pressure, u-velocity, and v-velocity. A simple way to create these matrices is by providing a “masking” matrix to the
NodePad function. Since there we now have interior geometry for this test-case, we need to set parts of the masking
matrix equal to a nonzero number (for example 1, in this case, see NodePad.m for more information) to the NodePad
function as follows:

3 bnd = [1 1 1 1];
5 Mask(1:round(app.Ny/3),1:round(app.Nx/5))=1;
6 Mask(end−round(app.Ny/3−1):end,1:round(app.Nx/5))=1;
7
9 [NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,bnd,1);
10 NodeRho=NodeP;
In this case, there are N x interior nodes in the x-direction, and N y interior nodes in the y-direction, minus the nodes
that are masked. Note that only the first fifth (Nx /5) of the domain in the x-direction is masked, and the top and
bottom third (Ny /3) of the domain in the y-direction is masked. The mask is only created for the interior nodes.
Therefore NodePad ”pads” the masking matrix to add the four boundaries conditions on each edge. The second input
to NodePad is a logical array, and tells the function which of the boundaries to pad ([left, right, bottom, top]), and the
last input is a logical scalar which checks for periodicity, automatically adding periodic boundaries if required.
Note that the ID of the first interior degree of freedom is N bcs + 1, where N bcs is a scalar variable containing the
number of boundary conditions, in this case N bcs = 5. Therefore the first interior node ID is 6, and this will be the
first unknown that will be solved (that is, the first row in the A matrix).
37
An example of NodeP and Nodeu for a N x = 10, N y = 9 domain is as follows:
 
2 5 5 5 5 5 5 5 5 5 5 3
 2 1 1 12 21 30 39 48 57 66 75 3 
 
 2 1 1 13 22 31 40 49 58 67 76 3 
 
 2 1 1 14 23 32 41 50 59 68 77 3 
 
 2 6 9 15 24 33 42 51 60 69 78 3 
 
N odeP =   2 7 10 16 25 34 43 52 61 70 79 3 

 2 8 11 17 26 35 44 53 62 71 80 3 
 
 2 1 1 18 27 36 45 54 63 72 81 3 
 
 2 1 1 19 28 37 46 55 64 73 82 3 
 
 2 1 1 20 29 38 47 56 65 74 83 3 
2 4 4 4 4 4 4 4 4 4 4 3
 
2 5 5 5 5 5 5 5 5 5 3
 2 1 1 12 21 30 39 48 57 66 3 
 
 2 1 1 13 22 31 40 49 58 67 3 
 
 2 1 1 14 23 32 41 50 59 68 3 
 
 2 6 9 15 24 33 42 51 60 69 3 
 
N odeu =   2 7 10 16 25 34 43 52 61 70 3 

 2 8 11 17 26 35 44 53 62 71 3 
 
 2 1 1 18 27 36 45 54 63 72 3 
 
 2 1 1 19 28 37 46 55 64 73 3 
 
 2 1 1 20 29 38 47 56 65 74 3 
2 4 4 4 4 4 4 4 4 4 3
With the boundary conditions chosen for this problem, however, the arrangement given above is not possible (see
section rules), hence the boundary condition numbers are re-arranged using the following code:
1 %Re−arrange pressure boundary node numbers suchs that open boundary has lowest ...
number
2 NodeP(find(NodeP==1))=−1;
3 NodeP(find(NodeP==3))=1;
4 NodeP(find(NodeP==−1))=3;
5 %Re−arrange u−velocity boundary node numbers such that open boundary has ...
highest number
6 Nodeu(find(Nodeu==3))=−1
7 ;Nodeu(find(Nodeu==5))=3;
8 Nodeu(find(Nodeu==−1))=5;
9 %Re−arrange v−velocity boundary node numbers such that Neumann boundary has ...
highest number
10 Nodev(find(Nodev==3))=−1;
11 Nodev(find(Nodev==5))=3;
12 Nodev(find(Nodev==−1))=5;
13 %Re−arrange density boundary node numbers such that Neumann boundary has ...
highest number
14 NodeRho(find(NodeRho==3))=−1;
15 NodeRho(find(NodeRho==5))=3;
16 NodeRho(find(NodeRho==−1))=5;
38
and gives the following results for NodeP and Nodeu:
 
2 5 5 5 5 5 5 5 5 5 5 1
 2 3 3 12 21 30 39 48 57 66 75 1 
 
 2 3 3 13 22 31 40 49 58 67 76 1 
 
 2 3 3 14 23 32 41 50 59 68 77 1 
 
 2 6 9 15 24 33 42 51 60 69 78 1 
 
N odeP =   2 7 10 16 25 34 43 52 61 70 79 1 

 2 8 11 17 26 35 44 53 62 71 80 1 
 
 2 3 3 18 27 36 45 54 63 72 81 1 
 
 2 3 3 19 28 37 46 55 64 73 82 1 
 
 2 3 3 20 29 38 47 56 65 74 83 1 
2 4 4 4 4 4 4 4 4 4 4 1
 
2 3 3 3 3 3 3 3 3 3 5
 2 1 1 12 21 30 39 48 57 66 5 
 
 2 1 1 13 22 31 40 49 58 67 5 
 
 2 1 1 14 23 32 41 50 59 68 5 
 
 2 6 9 15 24 33 42 51 60 69 5 
 
N odeu =   2 7 10 16 25 34 43 52 61 70 5 

 2 8 11 17 26 35 44 53 62 71 5 
 
 2 1 1 18 27 36 45 54 63 72 5 
 
 2 1 1 19 28 37 46 55 64 73 5 
 
 2 1 1 20 29 38 47 56 65 74 5 
2 4 4 4 4 4 4 4 4 4 5
Significant care must be taken with properly setting up the

boundary conditions
There are a number of conventions related to the boundary conditions. First, let us define the variables used for the
boundary conditions, in each case these variables are one-dimensional double arrays:
1 bcsDP: Pressure Dirichlet boundary conditions

2 bcsNP: Pressure Neumann boundary conditions
3 bcsOP: ID of Pressure Dirichlet boundary that should be an open boundary
4 bcsDu: U−velocity Dirichlet boundary conditions
5 bcsNu: U−velocity Neumann boundary conditions
6 bcsOu: ID of U−velocity Dirichlet boundary that should be an open boundary
7 bcsDv: V−velocity Dirichlet boundary conditions
8 bcsNv: V−velocity Neumann boundary conditions
9 bcsOv: ID of V−velocity Dirichlet boundary that should be an open boundary
The numbering convention comes about from the way that the node numbering matrices are used. The first N bcs
numbers are reserved for the boundary conditions. So, in this case, the first 5 numbers are reserved for boundaries.
bcsDP, then is an array that contains the VALUE of the Pressure Dirichlet boundary conditions with one exception
(for open boundaries), and bcsNP contains the VALUE of the Pressure Neumann boundary conditions. bcsOP contains
the IDs of Dirichlet boundary conditions that are actually open boundary conditions, and this is where the exception
comes about for Dirichlet boundaries. This may seem strange, but due to the way open boundary conditions are
implemented, they are first treated as Dirichlet boundaries and the matrices are modified afterwards to open boundary
conditions. The numbering rules are as follow: The lowest numbered boundary conditions must be of Dirichlet or
Open type, and the highest numbered boundary conditions must be of Neumann type. That is, there can never be
a Dirichlet boundary that has a larger number than a Neumann boundary. This explains why the Node numbering
matrices had to be renumbered in the grid-creation section. Open boundary conditions must be the highest numbered
Dirichlet boundary condition, that is, no true Dirichlet boundary condition may have a higher number than an
open boundary condition. Therefore, the numbering order (from lowest-numbered to highest numbered) is then as
follows:
39
1. Dirichlet
2. Open
3. Neumann
2
∂ φ
The open boundary condition used sets ∂n 2 = 0, where φ is any quantity and n is in the normal direction. The
reason for this boundary condition is that it is a very weak condition on the solution. It essentially requires that the
second derivative vanishes as the exit, that is, no curvature.
To implement the following boundary conditions
ρ = 0 ∀x 6= 20 on ∂Ω
∂ρ
= 0 ∀x = 20 on ∂Ω
∂n
∂P
= 0 ∀x 6= 20 on ∂Ω
∂n
∂2P
= 0 ∀x = 20 on ∂Ω
∂n2
u = 0 ∀y = 0|y = 1 on ∂Ω
u = 1 ∀x = 0 on ∂Ω
∂u
= 0 ∀x = 20 on ∂Ω
∂n
v = 0 ∀x 6= 20 on ∂Ω
∂v
= 0 ∀x = 20 on ∂Ω
∂n
we use the code below:
1 %Density
2 bcsDrho = [0 0 0 0];
3 bcsNrho = [0];
4 %Pressure
5 bcsDP = [0];
6 bcsNP = [0 0 0 0];
7 bcsOP = [1];
8 %Velocities
9 bcsDu = [0 1 0 0 0];
10 bcsNu = [];
11 bcsOu = [5];
12 bcsDv = [0 0 0 0];
13 bcsNv = [0];
5.6.4 Initial Conditions
To create the initial conditions, first ”coordinate” matrices that contain the coordinates of the control volume centers
are created as follows:
2 xp = linspace(dx/2, a−dx/2, Nx );yp = fliplr(linspace(dy/2, b−dy/2, Ny ));
3 xu = linspace( dx ,a−dx , Nx−1);yu = fliplr(linspace(dy/2, b−dy/2, Ny ));
4 xv = linspace(dx/2, a−dx/2, Nx );yv = fliplr(linspace(dy , b−dy, Ny−1));
5 [XP YP] = meshgrid(xp,yp);[XU YU] = meshgrid(xu,yu);[XV YV] = meshgrid(xv,yv);
note here we are using the MATLAB meshgrid function to create the matrices, and that these matrices are of the
size of the INTERIOR nodes only, that is, they do not go all the way up to the boundaries. Also note, the vertical
40
coordinate starts from the maximum value and decreases, and this is due to the node numbering convention described
in .
When initializing the variables, the boundary condition values are also stored in the vector of unknowns. Therefore,
when initializing the vector of unknowns, the values of the boundary conditions have to be included as well.
In this case there are interior masked nodes, which means that the created coordinate matrices will have too many
entries, and will not necessarily correspond to the IDs of the unknowns. Here the idsP, idsu, and idsv outputs
from the NodePad function needs to be used to select the correct elements of the coordinate matrices in order to
properly initialize.
Initializing the vector of unknowns using a uniform u-velocity is accomplished as follows:
1 rhoinit=@(X,Y)zeros(size(X));
2 pinit = @(X,Y) zeros(size(X));
3 uinit = @(X,Y) zeros(size(X));
4 vinit = @(X,Y) zeros(size(X));
5
6 rho=[bcsDrho' ;bcsNrho' ;rhoinit(XP(idsP), YP(idsP))];
7 P = [bcsDP; bcsNP; pinit(XP(idsP), YP(idsP))];
8 u = [bcsDu; bcsNu; uinit(XU(idsu), YU(idsu))];
9 v = [bcsDv; bcsNv; vinit(XV(idsv), YV(idsv))];
5.6.5 Misc. and Plotting
To plot the solution we can simply use surf(P(NodeP)) in the plotting script. To plot the solution at the correct
(x, y) coordinates for the interior we can use surf(XP,YP,P(NodeP(2:end-1,2:end-1))).
6 Troubleshooting
1. I get an error that MATLAB cannot find certain files?
Did you remember to run setup at the start?
2. I changed app.Nx and app.Ny when running one of the examples in TestCase, but now my solution is
wrong/blows up, and I get NaNs. What’s going on?
You are probably violating the CFL condition (See section CFL). Try reducing app.dt until your solution is
no longer wrong/blows up.
3. My Boundary conditions aren’t working? What’s going on?
Did you carefully read the rules concerning boundary conditions? See section rules.
4. None of the above worked? What do I do now?
Make sure that you are not inputting any non-physical parameter (for example, a negative dt). Check over
your code, and try a few more things, but if you are stuck, please try the forum or submit a bug report on
http://mseas.mit.edu/codes.
41
References
Cushman-Roisin, B. and Beckers, J.-M. (2011). Introduction to Geophysical Fluid Dynamics: Physical and Numerical
Aspects. Academic Press, second edition.
Guermond, J. and Minev, P. (2006). An overview of projection methods for incompressible flow. Comput. Methods
Appl. Mech. Engrg., 195:6011–6045.
LeVeque, R. J. (2002). Finite Volume Methods for Hyperbolic Problems. Cambridge University Press.
Özgökmen, T. M. and Chassignet, E. P. (2002). Dynamics of two-dimensional turbulent bottom gravity currents.
Journal of Physical Oceanography, 32:1460–1478.
Seibold, B. (2008). A Compact and Fast Matlab Code Solving the Incompressible Navier-Stokes Equations on
Rectangular Domains. MIT Course 18.086 Computational Science and Engineering II.
Van Leer, B. (1977). Towards the ultimate conservative difference scheme. IV. A new approach to numerical convec-
tion. J. Comput. Phys., 23(3):276 – 299.
42

Documentation 229 FV Framework

Uploaded by

Copyright:

Available Formats

Documentation 229 FV Framework

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Documentation 229 FV Framework

Uploaded by

Copyright:

Available Formats

2.

29 Finite Volume MATLAB Framework Documentation

Matt Ueckermann, Jing Lin, Sydney Sroka, Pierre Lermusiaux

April 12, 2018

2 Discretization: Schemes and Implementation 8

3 Download and Installation 18

4 Quick Start Guide 20

5 Setting Up New Examples 22

1 Problem Formulations: Incompressible N-S Equations

1.1 General 3-D N-S Equations

1.2 N-S Equations with Boussinesq Approximation

−∇p + ρg = −∇(p0,ref + ρ0 g · (x − x0,ref ) + p0 ) + (ρ0 + ρ0 )g = −∇p0 + ρ0 g (8)

Substituting (5) and (8) into (4) and dividing by ρ0 gives

1.3 2-D Flow Idealization and Nondimensionalization

A 2-dimensional flow in a horizontal plane is governed by

A 2-dimensional flow in a vertical plane is governed by

app.nu = ν1 , app.nu2 = ν2 , app.kappa = κ1 , app.kappa2 = κ2 ,

1.5.1 Bottom Gravity Current

After the nondimensionalization, the governing equations read

Four nondimensional numbers appear in the system:

νz gh3 β∆S νx kh2

Furthermore, the initial condition is set as

except at the salinity forcing region we have

2 Discretization: Schemes and Implementation

2.1 Advection Schemes

Or using the East, West, North, South naming convention:

2.1.1 Central Scheme

2.1.2 Upwind Schemes

φ̂e = 6/8φU − 1/8φU U + 3/8φD (23)

2.1.4 TVD Scheme

where the MC slope limiter Ψ(r) is defined as

2.1.5 Advection Schemes for Dynamically Orthogonal Equations

2.1.6 Mexed Routines

2.1.7 The CFL Condition

The Courant-Friedrichs-Lewy (CFL) condition is a stability criterion on the numerical solution of

2.2 Matrix Operators

2.2.1 Gradient Operator

2.2.2 Divergence Operator

2.2.3 Laplace Operator

2.3 Symmetry of the Discrete Laplace Operator

2.4 Consistency of Discrete Operators

2.5 Stokes Solvers

2.6 Solving Large Sparse Linear Systems

2.7 Solving Poisson Equations with All Neumann Boundary Conditions

In the incompressible Navier-Stokes equations, the pressure is determined only up to a constant.

∇ · ∇φ(x) = f (x), x ∈ Ω, (38)

is satisfied. This is because

3 Download and Installation

3.2 Installing on Linux

3.3 Installing on Mac

and you are ready to go!

3.4 Installing on Windows

1 >> mex −setup

Then you may compile the .cpp subroutines as follows:

4 Quick Start Guide

where ? should be replaced by any integer from 1-13.

1 >> help Example

1 >> app.nu = 0.1