CONSTRAINED RUNS ALGORITHM AS A LIFTING OPERATOR
FOR THE BOLTZMANN EQUATION

arXiv:1410.4399v1 [cs.CE] 16 Oct 2014

YNTE VANDERHOYDONC∗ AND WIM VANROOSE†
Abstract. Lifting operators play an important role in starting a kinetic Boltzmann model from
given macroscopic information. The macroscopic variables need to be mapped to the distribution
functions, mesoscopic variables of the Boltzmann model. A well-known numerical method for the
initialization of Boltzmann models is the Constrained Runs algorithm. This algorithm is used in
literature for the initialization of lattice Boltzmann models, special discretizations of the Boltzmann
equation. It is based on the attraction of the dynamics toward the slow manifold and uses lattice
Boltzmann steps to converge to the desired dynamics on the slow manifold. We focus on applying the
Constrained Runs algorithm to map density, average flow velocity, and temperature, the macroscopic
variables, to distribution functions. Furthermore, we do not consider only lattice Boltzmann models.
We want to perform the algorithm for different discretizations of the Boltzmann equation and consider
a standard finite volume discretization.
Key words. Lifting operator, initialization, missing data, kinetic Boltzmann models, macroscopic partial differential equations, finite volume discretization, Constrained Runs
AMS subject classifications. 76P05, 82C40, 35K45, 35K57

1. Introduction. Fluid dynamics can typically be described by different levels of accuracy. One distinguishes between micro-, meso-, and macroscopic scales.
Macroscopic partial differential equations (PDEs) model only a few low order velocity
moments and are therefore not that accurate to describe interactions between particles. Well-known macroscopic equations for the modeling of fluid dynamics are the
Euler and Navier–Stokes equations. However, when a detailed description is necessary, a microscopic model is used. These individual-based models take the particle
collision physics into account. The Boltzmann equation can be used to describe kinetic models, that are ubiquitous to model at this scale, by modeling distribution
functions in phase space.
We focus on space, time, and velocity discretizations of the Boltzmann equation.
The initialization of such a discretized Boltzmann equation requires a lifting operator
when only macroscopic information is available. This operator defines a mapping between macroscopic and microscopic/mesoscopic variables such that the distribution
functions of the discretized Boltzmann equation can be built from the given macroscopic information. The concept of a lifting operator in a multiscale context was
introduced by Kevrekidis et al. in the equation-free framework to couple different
scales in a dynamical system [15]. Similar ideas were introduced in the heterogeneous
multiscale methods framework that starts from a predetermined incomplete form for
the macroscopic model and estimates the needed data for the incomplete macroscopic
model from the microscopic model [6].
Lots of previous work in the literature is based on the initialization of lattice
Boltzmann models (LBMs), special discretizations of the Boltzmann equation [19, 31].
Lifting operators for LBMs are considered in [22, 23, 24, 26, 27, 28, 30] where a given
initial density is mapped to distribution functions. The Chapman–Enskog expansion
is built for such model problems in [22]. This analytical expansion is first introduced
∗ Dept. Mathematics and Computer Science, Universiteit Antwerpen, Middelheimlaan 1, 2020
Antwerpen, Belgium (ynte.vanderhoydonc@uantwerp.be).
† Dept. Mathematics and Computer Science, Universiteit Antwerpen, Middelheimlaan 1, 2020
Antwerpen, Belgium (wim.vanroose@uantwerp.be).

1

in [4] to solve the Boltzmann equation. Furthermore, [22], [23], [24], and [26] apply
the Constrained Runs (CR) algorithm to LBMs. Originally, the CR algorithm is introduced by Gear et al. [7] to map macroscopic initial variables to missing microscopic
variables for stiff singularly perturbed ordinary differential equations (ODEs). This
algorithm is based on the attraction of the dynamics toward the slow manifold. The
dynamics on this slow manifold can be parameterized by only macroscopic variables
such as the density. The higher order velocity moments become slaved functionals of
the density in the LBM context. [27] compares the Chapman–Enskog expansion with
the CR algorithm for hybrid model problems with given initial density. [28] constructs
the numerical Chapman–Enskog expansion as a lifting operator for these model problems. This numerical Chapman–Enskog expansion is based on a combination of the
Chapman–Enskog expansion and the CR algorithm. It finds the coefficients of the
Chapman–Enskog expansion numerically based on the CR algorithm. This reduces
the number of unknowns in the lifting since it only finds the coefficients of the expansion rather than the full state of distribution functions. A review of these lifting
operators for LBMs with an initial density as macroscopic variable is given in [30].
A generalization of the numerical Chapman–Enskog expansion is made in [29] where
both density and momentum are considered to be given at initialization.
In this paper, we apply the Constrained Runs algorithm on a finite volume discretization of the one-dimensional Boltzmann equation which maps macroscopic information, namely density, average flow velocity, and temperature, to distribution
functions of the discretized Boltzmann equation. The intention of this generalization
is to deal with multiple given macroscopic variables. We do not conserve only the
density but also average flow velocity and temperature to work towards more realistic
physical systems. This lifting operator is constructed to deal, for example, with missing data in a hybrid kinetic and macroscopic PDE model for laser ablation. In laser
ablation, a complex interplay between various processes determines the outcome of a
laser beam hitting on a solid target [3]. First, the target is heated by the laser and the
heat is transported over the material. The material then starts to melt. Above the
surface of the melted material there is evaporation and the particles escape according
to a certain velocity distribution. Due to the strong evaporation, the particles right
above the melt/gas interface are not in an equilibrium state. This thin layer above the
surface is known as the Knudsen layer which cannot be described by a macroscopic
PDE. It requires a detailed kinetic model such as the Boltzmann equation. Once away
from the surface the plume expansion can be described by a macroscopic PDE model.
The challenge is to couple the models describing the different physical processes in a
mathematical correct way.
We first applied the generalization of the CR algorithm in [29] for a small number of velocities to deal with the conservation of density and momentum in lattice
Boltzmann problems.
The outline of this paper is as follows. Section 2 contains the different levels
of description that will be used throughout the paper. The mesoscopic level uses
the discretized version of the Boltzmann equation while the macroscopic equivalent
PDEs are the Euler and Navier–Stokes equations. Section 3 describes the Constrained
Runs algorithm. The origin of the Constrained Runs algorithm is included in section
3.1. Section 3.2 contains the application of Constrained Runs to lattice Boltzmann
models. The generalization to different discretizations of the Boltzmann equation and
multiple conserved moments is outlined in section 3.3. Numerical results can be found
in section 4 where the lifting operator is tested in a setting of restriction and lifting.
2

A conclusion and outlook are included in section 5.
2. Different levels of description. This section introduces models with different levels of description that will be used throughout this paper. The mesoscopic
scale is described by a finite volume discretization of the Boltzmann equation and will
be discussed in section 2.1. The macroscopic scale uses discretized partial differential
equations (PDEs) to describe the evolution of a few low order velocity moments in
section 2.2.
2.1. Mesoscopic description. Kinetic models make use of the Boltzmann equation [19] that describes the evolution of a distribution function f (x, v, t) (function
space CR2 (D)) that counts the number of particles or individuals in point x ∈ Dx ⊂ Rn ,
n ∈ N0 , with a velocity v ∈ Dv ⊂ Rn , at time t ≥ 0. The equation is given by (without
external forces, with Bhatnagar–Gross–Krook (BGK) collision term [14])
∂
∂
f (x, v, t) + v ·
f (x, v, t) = ω(f eq (x, v, t) − f (x, v, t)).
∂t
∂x

(2.1)

The Maxwell–Boltzmann equilibrium distribution is

f eq (x, v, t) = n



m
2πkB T

d/2



e

−

m(v − u)2

2kB T


,

where kB is the Boltzmann constant, d the number of spatial dimensions, m the
molecular mass, n the number density, T the temperature, and u the average flow
velocity. The macroscopic variables are defined by [18]
Z
n(x, t) = f (x, v, t)dv,
Z
ρ(x, t) = mf (x, v, t)dv,
Z
ρ(x, t)u(x, t) = mvf (x, v, t)dv,
Z
1
T (x, t) =
mkv − uk2 f (x, v, t)dv,
(2.2)
dρR
with k.k the two-norm and R = kB /m the specific gas constant.
The BGK collision term approximation represents a relaxation towards equilibrium with a relaxation parameter ω and an associated time scale τ = 1/ω. Atomic
collisions are taken into account by the relaxation frequency ω [1] given by
ω=ρ

kB T
,
m µ

(2.3)

with µ the viscosity of the gas which is determined by

ωµ
T
µ = µref
.
Tref
ωµ represents the viscosity index of the considered gas [2], and µref and Tref are the
reference viscosity and temperature determined experimentally.
3

Consider the one-dimensional Boltzmann equation in space
∂f
∂f
+v
= ω(f eq − f ).
∂t
∂x
The velocity discretization of the Boltzmann equation results into a set of Nv linear
advection equations, with Nv the number of velocity directions and vi = v0 + i∆v,
i ∈ {0, . . . , Nv − 1} and ∆v = (vmax − vmin )/Nv .
∂fi
∂fi
+ vi
= ω(fieq − fi ).
∂t
∂x

(2.4)

The discretization of the equilibrium distribution function is outlined below and
is based on the algorithm described in [9]. We look for the equilibrium distribution
function fieq in the general form

fieq = A exp −B 2 (vi − D)2 ,
(2.5)

which has to approach the Maxwell–Boltzmann distribution as the increment ∆v
tends to zero. A, B, and D are derived from the conservation equations of mass,
momentum, and energy, which are defined by
∆v

NX
v −1

fi = n = ∆v

∆v

vi fi = nu = ∆v

∆v

i=0

NX
v −1

vi fieq (A, B, D),

i=0

i=0

NX
v −1

fieq (A, B, D),

i=0

i=0

NX
v −1

NX
v −1

(vi − u)2 fi =

NX
v −1
nkB T
(vi − u)2 fieq (A, B, D).
= ∆v
m
i=0

From this follows
A=

n
,
R0

R1 = 0,

R2 −

R0 kB T
= 0,
m

where
Rj = ∆v

NX
v −1
i=0


(vi − u)j exp −B 2 (vi − D)2 ,

from which A, B, and D are obtained through p
a Newton algorithm on this set of
equations with the initial approximation of B = m/(2kB T ) and D = u.
For the space and time discretization, we consider the one discussed in [16] by
Lafitte and Samaey. It represents a finite volume discretization of the Boltzmann
equation. Consider a grid which is uniform in time with time step ∆t and a spatial
grid with space step ∆x:
Cj = [xj− 12 , xj+ 21 ),

1 ≤ j ≤ N,

centered in xj with xj = j∆x,
Tk = [tk , tk+1 ),

k ≥ 0,

4

tk = k∆t.

By integrating (2.4) on a cell Mj,k = Cj × Tk , this results in

∆t 
φ(f )ki,j+1/2 − φ(f )ki,j−1/2
∆x

+ ∆tω fieq (xj , tk ) − fi (xj , tk ) ,

fi (xj , tk+1 ) =fi (xj , tk ) −

with upwind fluxes

φu (f )ki,j+1/2

=

(

vi fi (xj , tk ), vi ≥ 0,
vi fi (xj+1 , tk ), vi < 0,

or centered fluxes
φc (f )ki,j+1/2 = vi

fi (xj+1 , tk ) + fi (xj , tk )
.
2

Different discretizations of the Boltzmann equation are given in [17]. It represents
a similar explicit finite volume scheme and a linearized implicit scheme to compute
steady states. It focuses on the positivity of solutions, the conservation of moments,
and the dissipation of entropy.
2.2. Macroscopic description. When the phenomena of fluid dynamics that
one wants to describe are macroscopic, the fluid is regarded as a continuous medium.
The fundamental equations to describe the motion of the fluid are the Euler and
Navier–Stokes equations [5]. It can be shown that these are macroscopic equivalents
of the Boltzmann equation. For this, the Chapman–Enskog expansion can be used
[4, 31].
The Euler equations represent conservation of mass, momentum, and energy while
the Navier–Stokes equations extends these equations to include the viscosity of the
fluid.
2.3. Initialization. Suppose there is some initial data based on macroscopic
variables, like density, average flow velocity, and temperature. An important question
is how to start simulating a Boltzmann model, which is based on distribution functions, given initial macroscopic information. Starting the Boltzmann scheme includes
some arbitrariness since there are many possible distribution functions which have the
same velocity moments. We need to create a mapping from the macroscopic velocity
moments to distribution functions. For this a lifting operator is necessary. Section
3 describes the Constrained Runs algorithm that is used in this paper as a lifting
operator for the initialization of Boltzmann models. It is a known algorithm for the
initialization of lattice Boltzmann models. This paper generalizes it to apply it for
different discretizations of the Boltzmann equation and multiple conserved velocity
moments.
3. Constrained Runs algorithm. The Constrained Runs (CR) algorithm is a
numerical method based on the attraction of the dynamics toward the slow manifold.
The origin of the CR algorithm in systems of ODEs is described in section 3.1. The
application to LBMs with given initial density is included in section 3.2. Here, we
immediately present the generalized CR algorithm for LBM problems. The generalization to different discretizations of the Boltzmann equation and multiple conserved
moments (density, average flow velocity, and temperature) is discussed in section 3.3.
5

3.1. Origin of Constrained Runs algorithm. The CR algorithm finds its
origin in systems of ODEs [7]. A review of this work is given in [30] and summarized
below. Given system
dr(t)
= p(r(t), s(t)),
dt
ds(t)
= q(r(t), s(t)),
(3.1)
dt
where only the initial condition for r, namely r(0) = r0 , is given. The aim is to find
s(0) = s0 such that the initial condition (r0 , s0 ) lies on (or close to) the slow manifold.
The latter can be formulated by the function s0 = s(r0 ).
Gear et al. [7] proposed to obtain the s-value from equation
dm+1 s(t = 0)
= 0,
(3.2)
dtm+1
the smoothness condition, with m = 0, 1, . . . that can be approximated by a forward
difference
dm+1 s(t)
.
(3.3)
dtm+1
It can be shown that this difference approximation used in the CR algorithm can
be interpreted as a backward extrapolation [25]. It corresponds with a backward
extrapolation in time based on a polynomial of degree m that passes through the
values sk with k = 1, . . . , m + 1 while the known variable r is reset to its original
initial value r0 . The used coefficients of the forward finite difference formulas at time
t are listed in Table 3.1 for different degrees of m.
∆m+1 s(t) ≈ ∆tm+1

m
0
1
2
3

t + ∆t
1
2
3
4

t + 2∆t
0
−1
−3
−6

t + 3∆t
0
0
1
4

t + 4∆t
0
0
0
−1

Table 3.1
Coefficients of the forward finite difference formulas at time t for different degrees of m in (3.3).
d2 s(t)

For example, for m = 1, this formula corresponds to ∆t2 dt2 ≈ s(t) − 2s(t + ∆t) + s(t + 2∆t).
Using (3.2), this leads to s(t) = 2s(t + ∆t) − s(t + 2∆t) from which the coefficients can be obtained.

The general CR algorithm for a constant extrapolation, m = 0, is given in Algorithm 1.
Similar algorithms can be constructed for higher degrees of m by advancing the
model during more time steps and using (3.2) for different values of m.
3.2. Constrained Runs algorithm applied to lattice Boltzmann models.
The CR scheme [7] is a fixed point iteration scheme that computes the full state of
a microscopic time simulator on (or close to) the slow manifold corresponding to the
given macroscopic variables.
Consider an LBM problem in a D1Q3 setting (one spatial dimension, three velocity directions vi = ci ∆x/∆t, ci = i, i ∈ {−1, 0, 1}). The lattice Boltzmann equation
(LBE) is
fi (x + ci ∆x, t + ∆t) = (1 − ω)fi (x, t) + ωfieq (x, t).
6

Algorithm 1 Constrained Runs for a constant extrapolation in time in the system
of ODEs (3.1)
Require: Initial condition r(0) = r0
Choose s0 , norm k.k and a tolerance θ
repeat
Advance the model with one time step ∆t:
r1 and s1 at time t = ∆t
Difference approximation ∆s0 = s1 − s0
s0 ← s0 + ∆s0
Reset r to r0 , the given initial condition
until k∆s0 k < θ
This is a special discretization of the Boltzmann equation (2.1) [19, 31]. The equilibrium distributions are given by fieq (x, t) = 31 ρ(x, t), i ∈ {−1, 0, 1} [21] in which
the particle densityPρ(x, t) is defined as the zeroth order moment of the distribution
functions ρ(x, t) = i∈{−1,0,1} fi (x, t). These equilibrium distributions correspond to
a local diffusive equilibrium.
We now describe the application of Constrained Runs to LBMs [22, 24] when we
assume that a slow manifold exists in this context.
The CR algorithm sets a few LBM steps after which the density is reset. The
number of LBM steps is related to the order m and determines the accuracy of the
resulting lifted distribution function. Doing only one LBM step corresponds to a
constant extrapolation in time.
The CR procedure for LBMs iterates upon the higher order moments φ and ξ,
momentum and energy, given a density ρ0 = ρ(x, 0). It is equivalent to determinef :=
{f1 ; f0 ; f−1 } or M := {ρ; φ; ξ} since the velocity moments are related via moment
matrix M ∈ R3×3 :

 




1 1 1
ρ
f1
f1
 φ  =  1 0 −1   f0  = M  f0  ,
(3.4)
1
1
ξ
f
f
0
−1
−1
2
2
or, in shorthand, M = M f. Denote M 0 as


1 1 1
M 0 =  0 0 0  = diag(1, 0, 0)M.
0 0 0

M 0 represents the part of matrix M that produces the density which should be
conserved in the CR algorithm.
To apply CR, the missing moments are written as


φ
s=
,
ξ

a long vector s ∈ R2N , the variable r0 = ρ0 ∈ RN denotes the known initial condition,
with N the number of spatial grid points.
The vector s(k) denotes the k-th iterate of the CR algorithm and the iterations
are related by


s(k+1) = Cm r0 , s(k) ,
(3.5)
7

where Cm denotes one step of the CR algorithm and m is related to the order of the
time derivative that is set to zero in equation (3.2).
Below, this paper presents the CR algorithm in a general setting. We write the
CR iteration in terms of the distribution functions f rather than missing moments
φ and ρ. The algorithm generates
a sequence
of moments M(k) with k = 0, 1, 2, . . .


(0)
1
1
1 1
1
and M
= {ρ0 ; 3 − 3 ρ0 ; 2 3 + 3 ρ0 } the initial guess corresponding to the
(k)
equilibrium distribution function. In this sequence, the zeroth moment M0 is always
(k)
(k)
the initial density ρ0 and M1 = φ(k) and M2 = ξ (k) converge to the slaved state.
The iteration consists of two steps. First, M(k) is transformed to f (k) using (3.4)
and this is used as an initial state of an LBM that is evolved over m + 1 time steps.
This gives states f (k) (∆t), f (k) (2∆t), . . . , f (k) ((m + 1)∆t) that can be converted to
M(k) (∆t), M(k) (2∆t), . . . , M(k) ((m + 1)∆t). With these states we can generate
(k)pre

Mi

=

m+1
X

(k)

wj Mi (j∆t)

for

i = 1, 2,

(3.6)

j=1

the improved higher order moments where we add a superscript “pre” to denote that
this is the value before resetting the density. The underlying idea is that M(k)pre is
chosen such that combined with M(k) (j∆t) it satisfies the smoothness condition (3.2)
by letting the weights wj correspond to the coefficients of the forward finite difference
formula in (3.3) given in Table 3.1.
For the ease of analysis, we also apply the backward interpolation formula to the
zeroth moment. Now, we can write for all moments that
M(k)pre =

m+1
X

wj M(k) (j∆t).

(3.7)

j=1

To get the next iterate we introduce a reset step that sets the zeroth moment to ρ0
such that the density is conserved.


ρ0
 (k)pre 
M(k+1) = M1
.
(k)pre
M2
This reset is the second step in the iteration. The two



0 0 0
1
M(k+1) = 0 1 0 M(k)pre + 0
0 0 1
0
 (k)pre
0
0 0
= M −M f
+M f ,

steps can be written as
 
0 0
ρ0
0 0  0  ,
0 0
0

where f 0 is an initial distribution function that has ρ0 as density. Reordering leads
to

M f (k+1) = M − M 0 f (k)pre + M 0 f 0 ,
i
h

f (k+1) = M −1 M − M 0 f (k)pre + M 0 f 0 ,

f (k+1) = I − M −1 M 0 f (k)pre + M −1 M 0 f 0 .
8

In conclusion we can write
f (k+1) = I − M −1 M


0




m+1
X
j=1



wj f (k) (j∆t) + M −1 M 0 f 0 .

(3.8)

The convergence rate and stability of this iteration depend on the order m and the
coefficients wj from Table 3.1.
Note that the operator P := I − M −1 M 0 that appears in the first term of (3.8) is
a projection operator. Indeed, we can write M −1 M 0 = M −1 diag(1, 0, 0)M and it is
now easy to see that (M −1 M 0 )2 = M −1 diag(1, 0, 0)M M −1diag(1, 0, 0)M = M −1 M 0 .
As a result P 2 = P .
The interpretation of (3.8) is that the first term projects the state on the space
orthogonal to the space spanned by the lowest order moments. The second term adds
the original components in the space of the low order moments, here the density.
In general, equation (3.5) is nonlinear and the fixed point can be found by a
Newton iteration [26]. This means solving
gm (r0 , s) := s − Cm (r0 , s) = 0,

(3.9)

for a given macroscopic value r0 . Newton’s method gives an update to the guesses as
follows
s(k+1) = s(k) + δs(k) ,
where the corrections δs(k) are found by solving the linear system



∂gm 
A r0 , s(k) δs(k) =
r0 , s(k) δs(k) ,
∂s

∂Cm 
δs(k) ,
= I−
r0 , s(k)
∂s


= −gm r0 , s(k) ,

∂Cm
m
with A = ∂g
∂s the linearization (Jacobian matrix) of gm and ∂s the linearization
of Cm . I represents the identity matrix. The linearization can be estimated with the
help of the approximation

Aei ≈

gm (r0 , s + εei ) − gm (r0 , s)
,
ε

(3.10)

with ei the unit vector, i = 1, . . . , 2n and ε small.
3.3. Constrained Runs algorithm applied to general discrete Boltzmann equations. The previous section outlined the application of the Constrained
Runs algorithm to a lattice Boltzmann model with three velocities and where only
the density, the lowest order velocity moment, is given.
One question is how the method should be applied to other discretizations of the
Boltzmann equation, for example the finite volume discretization discussed in section
2.1. These discrete models have many more discrete velocities than the LBM. Consider again a one-dimensional spatial domain. Instead of three, we have Nv velocity
directions with velocities represented as {v0 , . . . , vNv −1 }. Furthermore, the models
9

can have, besides the density, also momentum, and temperature given as macroscopic
variables.
In principle (3.8) is easily generalized to a problem with Nv velocoties. The matrix
M — the matrix that transforms the distribution functions into velocity moments —
is now a transposed Vandermonde matrix


1
1
...
1
1
 v0
v1
. . . vq−2 vq−1 
 2

2
2
2 
 v0
v
.
. . vq−2
vq−1
1
M =
,


..
..


.
.
q−1
q−1
q−1
q−1
. . . vq−2 vq−1
v1
v0
with q = Nv . Define M 0 as

1
 v0



 k−1
M0 = 
v0
 0





0

1
v1
..
.

...
...

v1k−1
0
..
.

1
vq−2
..
.

k−1
. . . vq−2
...
0
..
.

0

...

0

1



vq−1 




k−1 
vq−1
,
0 



0

with k the number of conserved macroscopic variables, 1 ≤ k < Nv . M 0 represents
the part of matrix M that produces the lower order velocity moments which should be
conserved in the CR algorithm. The part corresponding to the higher order velocity
moments is replaced with zeros. When k = 3 we conserve density, momentum, and
temperature. These macroscopic variables are linked to the three lowest order velocity
moments of the distribution functions.
Since M is a transposed Vandermonde matrix and the inverse of this matrix is
required to write down the projection operator in the first term of equation (3.8), we
should deal with the inverse of such a matrix.
A lot of literature elaborates on this problem in the context of polynomial interpolation. However, in realistic applications more than hundred velocity directions
need to be taken into account. Even with known techniques, like using the Lagrange
polynomials for inverting the Vandermonde matrix [20] or by defining the moments
in Chebyshev polynomials [12], it is still impossible to compute the inverse with lots
of velocity directions. Similarly, the method discussed by Golub and Van Loan in
[8] to solve a Vandermonde system by comparing it to polynomial interpolation and
considering the Newton representation of the interpolating polynomial has analogous
issues.
These issues are further outlined below. When M −1 is not computed correctly,
the conservation of the given macroscopic variables cannot be guaranteed. This is
demonstrated in Example 1 that shows that P is not a projection operator when
M −1 is computed incorrectly. We deal with these issues in the remainder of this
section.
Example 1 (Eigenvalues of P := I − M −1 M 0 ). This example shows one of
the major issues of inverting matrix M to obtain projection operator P in the CR
algorithm. Since P is a projection operator the eigenvalues should be either zero or
one. Figure 3.1 contains the computed eigenvalues of this operator with Nv = 56,
10

vmin = −9.9875 · 103 m/s, vmax = 9.9875 · 103 m/s, and ∆v = (vmax − vmin )/Nv . This
figure shows that we lose the property of the eigenvalues being either zero or one.
3
2.5

eig(P )

2
1.5
1
0.5
0
−0.5
0

10

20

30
1 : Nv

40

50

60

Fig. 3.1. Computed eigenvalues of operator P = I −M −1 M 0 with Nv = 56, vmin = −9.9875·103
m/s, vmax = 9.9875 · 103 m/s, and ∆v = (vmax − vmin )/Nv . P should be a projection operator thus
we expect the eigenvalues to be either zero or one. This figure shows that we lose this property when
using the inverse of matrix M in the numerical experiment.

Example 1 illustrates the difficulties that originate from inverting matrix M .
It is clear that the loss of accuracy leads to eigenvalues that deviate from zero or
one. This has important consequences for the Constrained Runs algorithm since the
conservation of the lower order velocity moments is required. And this is not longer
the case with incorrect projection operators.
Example 2 incorporates the Chebyshev polynomials to compute M −1 . Similarly,
we lose the property of eigenvalues being either zero or one.
Example 2 (Eigenvalues of P := I −M −1 M 0 by defining the moments in Chebyshev polynomials). This example checks the property that P should be a projection
operator when P is based on defining the velocity moments in Chebyshev polynomials. Figure 3.2 contains the computed eigenvalues of this operator with Nv = 56,
vmin = −9.9875 · 103 m/s, vmax = 9.9875 · 103 m/s, and ∆v = (vmax − vmin )/Nv . This
figure shows that we lose the property of the eigenvalues being either zero or one.
However, the use of M −1 can be circumvented. The projection on the space
of higher order moments can be written as a subtraction of the components in the
direction of the lower order moments. This then only requires orthogonal vectors that
span the space of lower order moments.
T
Let us look at the QR factorization of M 0 (:, 1 : k) where k is the number of
conserved macroscopic variables:
T

M 0 (:, 1 : k) = QR,
where Q contains k orthogonal columns and R is a k by k upper triangular matrix.
When density, momentum, and energy are conserved in a one-dimensional problem R
is a 3×3 matrix and Q contains three columns. Let us define M̂ 0 := [QT ; 0]. We now
have that M 0 (1 : k, :) = RT M̂ 0 (1 : k, :). The space spanned by the columns of M̂ 0 now
spans the same space as the original low order moments.
The CR iteration f (k+1) = P̂ f (k)pre + M̂ 0 f 0 , where P̂ projects on the space of
the higher order moments, can now be written as
T

f (k+1) = (1 − M̂ 0 M̂ 0 )f (k)pre + M̂ 0 f 0 ,
11

(3.11)

10
0
−10

eig(P )

−20
−30
−40
−50
−60
−70
0

10

20

30
1 : Nv

40

50

60

Fig. 3.2. Computed eigenvalues of operator P = I −M −1 M 0 with Nv = 56, vmin = −9.9875·103
m/s, vmax = 9.9875 · 103 m/s, and ∆v = (vmax − vmin )/Nv . P should be a projection operator thus
we expect the eigenvalues to be either zero or one. This figure shows that we lose this property when
using the inverse of matrix M in the numerical experiment when M is based on defining the velocity
moments in Chebyshev polynomials.

which avoids the inverse of the moment matrix.

(3.11) can be written because

T

M̂ T M̂ 0 = M̂ 0 (:, 1 : k)M̂ 0 (1 : k, :), where M̂ is the matrix that would have been
obtained by a QR factorization of the full moment matrix M .
We will demonstrate that an iteration using M̂ will preserve the same moments
as the original procedure based on M . However, the procedure based on M̂ is more
stable, which can be observed from the eigenvalues of the projection operator. Example 3 shows these eigenvalues for a similar numerical experiment as in Example 1
and Example 2.
Example 3 (Eigenvalues of P̂ := I − M̂ T M̂ 0 ). The eigenvalues of P̂ should be
either zero or one since P̂ is a projection operator. Figure 3.3 shows the computed
eigenvalues of this operator with Nv = 56, vmin = −9.9875 · 103 m/s, vmax = 9.9875 ·
103 m/s, and ∆v = (vmax − vmin )/Nv . The property of the eigenvalues being either
zero or one is fulfilled for this numerical experiment.
1.2
1

eig(P̂ )

0.8
0.6
0.4
0.2
0
−0.2
0

10

20

30
1 : Nv

40

50

60

Fig. 3.3. Computed eigenvalues of operator P̂ = I − M̂ T M̂ 0 with Nv = 56, vmin = −9.9875·103
m/s, vmax = 9.9875 · 103 m/s, and ∆v = (vmax − vmin )/Nv . P̂ should be a projection operator thus
we expect the eigenvalues to be either zero or one. This figure shows that operator P̂ fulfills this
property for the numerical experiment.
12

Example 3 shows that this alternative setting fulfills the property of eigenvalues
being zero or one for projection operator P̂ . We check the effect of this redefinition
— using P̂ instead of P — on the CR algorithm.
Now, we check whether this algorithm conserves the same lower order velocity
moments, the k conserved macroscopic variables



M 0 (1 : k, :)f (k) = M 0 (1 : k, :) P̂ f (k)pre + (I − P̂ )f 0 ,

= M 0 (1 : k, :)P̂ f (k)pre + M 0 (1 : k, :)f 0 − M 0 (1 : k, :)P̂ f 0 ,

= RT M̂ 0 (1 : k, :)P̂ f (k)pre + M 0 (1 : k, :)f 0 − RT M̂ 0 (1 : k, :)P̂ f 0 .

The first and last term of the right-hand side vanish because of the definition of the
projection space Ŵ of operator P̂ . It follows that
M 0 (1 : k, :)f (k) = M 0 (1 : k, :)f 0 .

This confirms the correspondence between the lower order moments in the CR algorithm and those of f 0 , the initial distribution function.

4. Numerical results. We illustrate the proposed method with the help of a
model problem coming from laser ablation. Here a Boltzmann model is studied with
a left boundary condition that models the evaporation of material from a heated
surface, while the right boundary models the ambient gas. We refer to [9] for a
detailed discussion of the physics and the Boltzmann model. This section contains
numerical results which test the generalized CR algorithm in a setting of restriction
and lifting.
Example 4 (One-dimensional Helium problem). As a model problem, we consider the one-dimensional laser ablation problem with Helium as a background gas.
The ambient gas parameters pa , Ta , na , ρa , and ua represent the pressure, temperature, number density, mass density, and average flow velocity while ps , Ts , ns , ρs , and
us are the surface parameters. These parameters, presented below in (4.1), are used
for initialization and to obtain boundary conditions. These are similar to the parameters presented in [9]. The boundary conditions at the surface are based on the
equilibrium distribution of the surface parameters while the outer boundary is placed
far enough from the interface such that the vapor is in equilibrium there.

pa = 0.10132500 · 106 P a,

Ta = 0.30000785 · 103 K,

ρa = mna ,

ua = 0,

pa
,
0.3

ρs = mns ,

us = 0.

ps =

Ts =

Ta
,
0.2

ns =

na =

pa
,
kB Ta

ps
,
kB Ts
(4.1)

13

We focus on the finite volume discretization, outlined in section 2.1, with discretization
parameters listed below in (4.2).
r
2kB Ts
N = 1600, Nv = 56, u0 =
, vmin = −4u0 , vmax = 4u0 ,
m


1
∆v
∆v
vmax − vmin
,
, λ= √
, v = vmin +
: ∆v : vmax −
∆v =
Nv
2
2
2πd2 ns


L
h
h
L = 30000λ, h = , x =
: h : + (N − 1)h ,
N
2
2
ρ(x, t)kB /mT (x, t)
,
µref = 0.19 · 10−4 P a · s, Tref = 273.15 K, ω(x, t) =
µref (T (x, t)/Tref )0.66
1
,
(4.2)
∆t = 0.9
max(v)/h + max(ω(x, 0))
with d the gas-kinetic molecular diameter and λ the mean free path. The time step
∆t is defined in such a way for stability reasons [17].
Since the length of the domain determines the discretization parameters, we can
also consider
L = 30λ,

(4.3)

which is a valid choice for the length since the thickness of the nonequilibrium layer
is about 10-20 mean free path lengths [9]. The remaining discretization parameters
are the same as those presented in (4.2).
These different domain lengths determine the possibility of using the CR algorithm as a lifting operator. This is one of the issues discussed in section 4.1, which is
related to the convergence rate of CR.
4.1. Comparison of P and P̂ in the CR algorithm. The following illustrates the issues of using the inverse of matrix M in the projection operator of the
Constrained Runs algorithm. Figure 4.1 shows both the spectra of the Jacobian matrices of the CR map with operators P , which uses the inverse of M , and P̂ , which
avoids this inverse, respectively in the left and right figure. The spectra correspond
to a Constrained Runs algorithm with m = 0. The eigenvalues should remain within
the unit circle for the algorithm to be stable. Since we are focusing on a generalization of Constrained Runs, we are not only taking the density as a given macroscopic
variable but also average flow velocity and temperature. This example includes 56
discrete velocities while preserving the three lowest order velocity moments for the
Helium problem presented in Example 4. The number of grid points is N = 200 as
opposed to the general parameters presented in (4.2). This reduction in the number
of variables is necessary to make the computation of the Jacobian matrix of the CR
map possible.
These figures show a clear stability with projection operator P̂ while an instability
occurs when using the CR algorithm with P due to the difficulties of computing the
inverse of M .
Another issue that is illustrated by the spectrum of the Jacobian matrix of the CR
algorithm is the convergence rate of Constrained Runs. Figure 4.2 shows the spectrum
of the Helium model problem presented in Example 4 with parameters listed in (4.3),
where the ratio of the grid distance ∆x to the mean free path is reduced. Again, the
CR algorithm with operator P shows a clear instability while P̂ stabilizes the method.
14

1

5

Im(µ)

Im(µ)

0.5

0

0

−0.5

−5
−12

−10

−8

−6
Re(µ)

−4

−2

−1

0

−1

−0.5

0
Re(µ)

0.5

1

Fig. 4.1. The spectrum of the Jacobian matrix of the Constrained Runs map with m = 0. The
left figure uses operator P while the right figure is based on projection operator P̂ . The parameters
are listed in (4.2) of Example 4 with a smaller number of grid points, namely N = 200, to make the
computation of the Jacobian matrix possible. The unit circle, which should contain the eigenvalues
to obtain a stable method, is shown in dashed lines. A clear instability occurs in the left figure when
the original projection operator P = 1 − M −1 M 0 is used in the CR algorithm.

However, the convergence is slower with these parameters since the spectral radius,
or equivalently the asymptotic convergence factor, is larger in Figure 4.2 compared to
Figure 4.1.
1

5

Im(µ)

Im(µ)

0.5

0

0

−0.5

−5
−12

−10

−8

−6
Re(µ)

−4

−2

−1

0

−1

−0.5

0
Re(µ)

0.5

1

Fig. 4.2. Similar as Figure 4.1 but now with parameters listed in (4.3), where there ratio of
∆x to the mean free path is much smaller.

4.2. Test Constrained Runs algorithm in Example 4. This section tests
the Constrained Runs algorithm on the model problem presented in Example 4. We
perform 10000 Boltzmann time steps on the initial state based on the ambient parameters presented in Example 4 and parameters listed in (4.2). The distribution
functions are rescaled with the mass m for numerical reasons due too small numbers
to perform the numerics on. This results in a reference distribution function fc that
is plotted in Figure 4.3. As can be seen, a traveling wave emerges in the domain. The
corresponding equilibrium distribution is f eq .
We test the CR algorithm as a lifting operator that maps density, average flow
velocity, and temperature to distribution functions. The CR algorithm is combined
15

Fig. 4.3. Plot of the reference distribution function fc obtained after 10000 time steps on the
initial state with parameters presented in Example 4. A traveling wave emerges in the domain. The
left boundary condition models the evaporating material in laser ablation. The right boundary is
determined by the parameters of the ambient gas. The vertical axis is rescaled with the mass for
numerical reasons as described in the text.

with Newton’s method to ensure stability. It uses a GMRES algorithm to estimate
the Jacobian matrix in Newton’s method. The parameters of the GMRES algorithm
are default parameters in Matlab. The convergence of GMRES might be improved if
a preconditioner is used similar to [13].
The lifting operator can be evaluated by restricting the reference distribution
function fc to its macroscopic variables and lift them back to a distribution function
f by using the lifting operator. The resulting f is compared with fc with the help of
the two-norm kf − fc k, which is shown in Table 4.1. These results can be compared
to kf eq − fc k which is equal to 6.4940e-007. They are based on Newton’s method
with a tolerance value of 1.0e-010. It might improve with a stricter tolerance.
Table 4.1
The error kf − fc k to test the CR algorithm as a lifting operator (combined with Newton’s
method) for various orders of accuracy. The reference distribution function fc is obtained for the
model problem in Example 4 by performing 10000 Boltzmann time steps on the initial state.
Order CR
0
1
2
3

kf − fc k
1.0428e-006
1.6413e-008
6.1629e-010
4.1965e-010

c
.
Figure 4.4 presents a log plot of the relative errors f f−fc (top left) and f −f
fc
c
These results are plotted in function of spatial grid points x and velocities v. f
corresponds to the distribution functions based on lifting with the Constrained Runs
algorithm of order m = 0 (top right), m = 1 (middle left), m = 2 (middle right),
m = 3 (bottom left), and m = 4 (bottom right). We show a narrow range of velocity
directions since the outer velocities bring the distribution functions to zero. Where
the figures show no results (white area in figures), |f − fc | is exactly equal to zero,
which makes it impossible P
to create the log plot.
P
Furthermore, we plot i |fieq − (fc )i | and i |fi − (fc )i | in Figure 4.5 to get an
idea on the construction of the hybrid domain.
This suggests the use of a Boltzmann model near the traveling wave while a PDE
eq

16

m=0

m=1

m=2

m=3

m=4

f eq −f

f −f

c
Fig. 4.4. Log plot of the relative error
(top left) and similarly of f c where fc repfc
c
resents the reference distribution function after 10000 time steps on the initial state with parameters
eq
presented in Example 4. f is the corresponding equilibrium distribution after 10000 steps and f
the distribution functions based on lifting with the Constrained Runs algorithm of order m = 0 (top
right), m = 1 (middle left), m = 2 (middle right), m = 3 (bottom left), and m = 4 (bottom right).
Where the figures show no results (white area in figures), |f − fc | is exactly equal to zero, which
makes it impossible to create the log plot.

17

domain can be used outside this part of the domain. Extra care is also needed at the
boundaries where the surface parameters determine the state.
Example 5 (Helium problem including radial velocities). This example considers the laser ablation problem from [9] with Helium as a background gas in a
three-dimensional domain described through radial velocities, which represent the axial
symmetry. pa , Ta , na , ρa , and ua again correspond to the ambient parameters pressure, temperature, number density, mass density, and average flow velocity while
ps , Ts , ns , ρs , and us are the surface parameters which are given below. These are
similar to the parameters presented in [9].
pa = 0.10132500 · 106 P a,

Ta = 0.30000785 · 103 K,

ρa = mna ,

ua = 0,

pa
,
0.3

ρs = mns ,

us = 0.

ps =

Ts =

Ta
,
0.2

ns =

na =

pa
,
kB Ta

ps
,
kB Ts

The discrete distribution function is defined on a two-dimensional grid with velocities
(vz )i = (vz )0 +i∆vz , i ∈ {0, . . . , Nz −1} and (vr )i′ = (vr )0 +i′ ∆vr , i′ ∈ {0, . . . , Nr −1},
specified in the cylindrical domain of the velocity space with (vz )min < vz < (vz )max
and 0 < vr < (vr )max . vr = [(vx )2 + (vy )2 ]1/2 is the radial velocity, which takes
into account the axial symmetry. The axial and radial velocity increments are ∆vz =
((vz )max − (vz )min )/Nz and ∆vr = (vr )max /Nr , respectively, (vz )0 = (vz )min + ∆vz /2
and (vr )0 = ∆vr /2. The discretization parameters now correspond to
r
2kB Ts
,
N = 1600, Nz = 56, Nr = 24, u0 =
m
(vz )max − (vz )min
,
(vz )min = −4u0 , (vz )max = 4u0 , ∆vz =
Nz


∆vz
∆vz
, (vr )min = 0,
: ∆vz : (vz )max −
v z = (vz )min +
2
2


(vr )max
∆vr
∆vr
(vr )max = 3u0 , ∆vr =
,
, vr =
: ∆vr : (vr )max −
Nr
2
2


L
h
h
1
, L = 30000λ, h = , z =
: h : + (N − 1)h ,
λ= √
2
N
2
2
2πd ns
ρ(z, t)kB /mT (z, t)
µref = 0.19 · 10−4 P a · s, Tref = 273.15 K, ω(z, t) =
0.66 ,
µref (T (z, t)/Tref )
1
.
(4.4)
∆t = 0.9
max(v z )/h + max(ω(z, 0))
Including radial velocities requires a redefinition of the velocity moment matrix M .
In general, tensors are necessary to describe higher order velocity moments [18]. The
lower order moments, density, average flow velocity, and temperature are given in
(2.2).
We start by checking the eigenvalues of P̂ with this extended moment matrix. The
eigenvalues should still be either zero or one since P̂ is a projection operator. Figure
4.6 shows the computed eigenvalues of this operator for Example 5 which shows only
eigenvalues equal to zero or one.
4.3. Test Constrained Runs algorithm in Example 5. This section tests
the Constrained Runs algorithm on the model problem presented in Example 5. It is
18

−7

8

−6

x 10

2

m=0

x 10

7

|fi − (fc )i |

1.5

5

1

i

4
3

P

P

i

|fieq − (fc )i |

6

2

0.5

1
0
0

1

2
x

−8

3.5

3

0
0

4
x 10

−9

m=1

x 10

1

−3

1.2

3

2
x

3

4
−3

x 10

m=2

x 10

1

|fi − (fc )i |

2

0.8
0.6

i

1.5

P

P

i

|fi − (fc )i |

2.5

0.4

1
0.2

0.5
0
0

1
−10

9

2
x

3

0
0

4
x 10

−9

m=3

x 10

1

−3

1.4

8

2
x

3

4
−3

x 10

m=4

x 10

1.2
1
|fi − (fc )i |

6
5
4

i

0.8
0.6

P

P

i

|fi − (fc )i |

7

3

0.4

2
0.2

1
0
0

1

2
x

3

4
−3

x 10

0
0

1

2
x

3

4
−3

x 10

P
P
eq
Fig. 4.5. Plot of the absolute differences
i fi − (fc )i (top left) and
i |fi − (fc )i | where
fc represents the reference distribution function after 10000 time steps on the initial state with
parameters presented in Example 4. f eq is the corresponding equilibrium distribution and f the
distribution functions based on lifting with the Constrained Runs algorithm of order m = 0 (top
right), m = 1 (middle left), m = 2 (middle right), m = 3 (bottom left), and m = 4 (bottom right).
Note the different scales in the figures. These figures are shown to get an idea on the construction
of the hybrid domain.

19

−15

1.5

1.2
1

x 10

1

0.5

0.6

eig(P̂ )

eig(P̂ )

0.8

0.4

0

0.2

−0.5
0
−0.2
0

200

400

600
800
1 : Nz ∗ Nr

1000

1200

1400

−1
1

1.5

2
1:3

2.5

3

Fig. 4.6. Computed eigenvalues of operator P̂ = I − M̂ T M̂ 0 (left) with Nz = 56 and Nr = 24
as presented in Example 5. P̂ should be a projection operator thus we expect the eigenvalues to be
either zero or one. We see, by investigating the region around zero, that operator P̂ fulfills this
property (right).

based on parameters suggested by Gusarov et al. [9] for the laser ablation problem.
The boundary conditions of the Boltzmann discretization are equilibrium distributions
of the surface parameters at the surface and equilibrium distributions of the ambient
parameters at the outer boundary. We perform 10000 Boltzmann time steps on the
initial state based on the ambient parameters presented in Example 5. The distribution functions are rescaled with the mass m and weights 2π∆vz ∆vr , corresponding to
the radial velocity directions. The traveling wave that emerges in the domain moves
from the surface, on the left, towards the ambient boundary, at the right. The resulting reference distribution function is denoted as fc . The corresponding equilibrium
distribution is denoted as f eq .
We test the performance of the CR algorithm as a lifting operator that maps
density, average flow velocity, and temperature to distribution functions. The CR
algorithm is combined with Newton’s method to ensure stability. It uses a GMRES
algorithm to invert the Jacobian matrix in Newton’s method and the matrix-vector
product of the Jacobian matrix is estimated using finite differences.
f −f
f eq −fc
(top left) and f c .
Figure 4.7 presents a log plot of the relative errors
fc
c
These results are plotted in function of spatial grid points z and velocities in the
z-direction v z . f corresponds to the distribution function based on lifting with the
Constrained Runs algorithm of order m = 0 (top right), m = 1 (bottom left), and
m = 2 (bottom right). It shows a narrow range of velocity directions v z since the
outer velocities bring the distribution functions to zero. Furthermore, the results
include radial velocities v r . The errors are based on sums of distribution functions
over the different radial velocity directions.
There are two remarks that illustrate the drawbacks of the Constrained Runs
algorithm.
Remark 1. A serious drawback is the dependence of the convergence rate of
the Constrained Runs algorithm on the parameters of the Boltzmann model. This is
illustrated in Figures 4.1 and 4.2. This makes it computationally not feasible to use
the CR algorithm as a lifting operator for certain choices of the parameters of the
Boltzmann model, especially when the grid resolution is smaller than the mean free
path. This is illustrated in Figure 4.8, where the number of outer GMRES iterations
20

m=0
0

0

−5

−5

−10

−10

−15

−15

−20
5000

−20
5000
8

vz

8

6

0

4
−5000 0

6

0

4

−3

x 10

2

vz

z

m=1

−5000 0

−3

x 10

2
z

m=2

0

0

−5

−5

−10

−10

−15

−15

−20
5000

−20
5000
8

vz

8

6

0

4
−5000 0

2

6

0

4

−3

x 10

z

f eq −f

vz

−5000 0

2
z

f −f

c
Fig. 4.7. Log plot of the relative error
(top left) and similarly of f c where fc repfc
c
resents the reference distribution function after 10000 time steps on the initial state with parameters
eq
presented in Example 5. f
is the corresponding equilibrium distribution and f the distribution
functions based on lifting with the Constrained Runs algorithm of order m = 0 (top right), m = 1
(bottom left), and m = 2 (bottom right). Furthermore, the results include radial velocities v r . The
errors are based on sums of distribution functions over the different radial velocity directions.

is shown in function of N , the number of spatial grid points for m = 0 (circle), m = 1
(asterisk), m = 2 (plus sign), m = 3 (point), m = 4 (cross), and m = 5 (triangle).
It is clear that the convergence rate of CR is determined by the parameters of the
Boltzmann model, especially the ratio of ∆x to the mean free path λ, is important.
To obtain a good convergence rate of Constrained Runs in the previous examples,
we chose the length of the domain L = 30000λ in (4.2) and (4.4). We would rather
like to have L = 30λ such that the nonequilibrium layer is well described. It typically
has a thickness of 10 to 20 mean free path lengths [9]. This means that the current
grid step equal to H = 30000λ/N has to be reduced to h = 30λ/N . This can be
interpreted as refinements from BoltzmannH to Boltzmannh .
However, to refine Boltzmann models additional operators are necessary that
transfer information between different scales. There is some literature on the coupling of lattice Boltzmann models with different grid resolutions. Guzik et al. present
a space-time interpolation method to couple different grid resolutions for lattice Boltzmann models [10, 11]. Here the first nonequilibrium term is matched between the two
Boltzmann models.
There is no literature for general discrete Boltzmann models with multiple velocities, to the best of our knowledge.
21

−3

x 10

1200

100

1000

80
GMRES

600

iter

iterGMRES

800

400

60
40
20
0

200

−20
0
0

500

1000
N

1500

2000

500

1000
N

1500

Fig. 4.8. Plot of the number of outer GMRES iterations in function of N , the number of spatial
grid points for m = 0 (circle), m = 1 (asterisk), m = 2 (plus sign), m = 3 (point), m = 4 (cross),
and m = 5 (triangle). with parameters presented in Example 4. The plot on the right represents a
zoom of the left figure. This shows that the convergence rate of CR depends on the parameters of
the Boltzmann model.

If it is possible to transfer between the Boltzmann subdomains with different grid
sizes, denoted BoltzmannH and Boltzmannh , then the following scenario would be
possible.
First, there is a lifting operator that transfers information between BoltzmannH
and PDEH where H is large. This has the advantage that this is better for the convergence rate of the Constrained Runs algorithm since the grid sizes determine the
spectral radius as highlighted in section 4.1. Then the Boltzmann model is systematically refined to reach the correct resolution to resolve the non-equilibrium layer. An
overall analysis is required to maintain the accuracy at each boundary.
Remark 2. Another drawback of the CR algorithm in general is the computational expense since it is based on performing Boltzmann steps to determine the
Jacobian matrix in Newton’s method. In a similar way, when a matrix-free method is
used such as GMRES the convergence is slow without preconditioning.
Remark 3. The Knudsen number Kn , which is the ratio of the mean free path
and the feature length, is often used to classify flow regimes [18]:
• Kn < 0.01: hydrodynamic regime,
• 0.01 < Kn < 0.1: slip flow regime,
• 0.1 < Kn < 10: transition regime,
• Kn > 10: free molecular flow.
This classification is used as a rule of a thumb to determine which mathematical
description is necessary.
In laser ablation, the surface gets hot and ejects particles. At the surface a
nonequilibrium layer emerges that requires a full Boltzmann description. While away
from the surface, a PDE model can be used to simulate the flow. Under what conditions can we lift information from the PDE model to information required by the
Boltzmann model? Can this be linked to the classification with the Knudsen number?
To know the link between the microscopic Boltzmann model and the macroscopic
PDE models, a Chapman–Enskog expansion is used. This assumes that the Knudsen
number is a small parameter (to achieve the Navier–Stokes equations). Indeed, in the
hydrodynamic regime, the Knudsen number is a small parameter.
This Knudsen number should help us to determine the conditions when to switch
22

from Boltzmann to PDE model in the hybrid domain. However, this is, so far, not
investigated in detail since the Knudsen number depends on the feature length that is
hard to estimate from region to region in the domain.
5. Conclusions. Many problems based on the Boltzmann equation require the
conversion of moments to the corresponding distribution function. One motivating
example is laser ablation where a material plume is ejected from a melting material
that is heated by a laser. To describe this process accurately a lifting operator is
required that maps hydrodynamic moments, namely density, average flow velocity,
and temperature, to distribution functions of the kinetic Boltzmann model at the
interface between the melted material and the gas.
This paper extends the applicability of the Constrained Runs (CR) algorithm
that lifts these moments to distribution functions to general discretizations of the
Boltzmann equation. Previously the CR algorithm was only used to lift lattice Boltzmann models, where it initializes or couples different models together. In this paper,
we have focused on lifting with CR in a finite volume discretization of the Boltzmann equation, but we believe that the results of the CR algorithm can be applied
to different discretizations.
The main difficulty encountered in this paper is that the straightforward formulation of CR for Boltzmann models requires the inverse of the moment matrix.
This is a very ill-conditioned matrix with the properties of a Vandermonde matrix
and traditional numerical methods to invert this matrix fail. In this paper we have
reformulated the algorithm such that this inversion can be avoided all together.
With this new formulation the method can conserve multiple macroscopic variables and include multiple velocity directions. The paper includes numerical results
that test the restriction and lifting. We perform 10000 Boltzmann time steps on the
initial state of the laser ablation problem to create a reference distribution function.
This reference solution is restricted to its macroscopic variables, density, average flow
velocity, and temperature, and lifted back to distribution functions to test the generalized CR algorithm. The error with the original reference distribution can now be
reduced by increasing the order of the CR algorithm.
A remaining drawback of the CR algorithm is the computational cost since it is
based on the time scale that a Boltzmann simulation needs to reach the slow manifold.
Another drawback is the dependence of the convergence rate of CR on the parameters
of the Boltzmann domain: length of the domain expressed as mean free path lengths,
step sizes, . . .
In [28] and [29] we have applied, for LBMs, the CR on distribution functions
represented as the first few terms of the Chapman–Enskog expansion. The iteration
then determines the coefficients of the expansion rather than the distribution function
itself. This significantly reduces the size of the problem and makes it easy to solve
the implicit, higher order problems. In the future, it might be possible to apply this
technique for general discretizations of Boltzmann models where multiple moments are
conserved since this paper generalized the underlying CR algorithm. However, when
a general Maxwell–Boltzmann equilibrium is used and multiple moments are given,
the Chapman–Enskog expansion becomes very complicated with the derivatives of the
various moments. At this moment it is unclear if this numerical Chapman–Enskog
technique is applicable in this context.
Acknowledgments. This work is supported by research project Hybrid macroscopic and microscopic modelling of laser evaporation and expansion, G.017008N,
funded by ‘Fonds Wetenschappelijk Onderzoek’ together with an ‘ID-beurs’ of the
23

University of Antwerp. Furthermore, we would like to thank Annemie Bogaerts and
David Autrique of the Plasmant group of the University of Antwerp for their help in
defining a model problem for the laser ablation problem.
REFERENCES
[1] K. Aoki, P. Degond, L. Mieussens, Numerical simulations of rarefied gases in curved channels:
thermal creep, circulating flow and pumping effect, Commun. Comput. Phys., 6 (2009),
pp. 919-954.
[2] G.A. Bird, Molecular gas dynamics and the direct simulation of gas flows, Oxford Sci. Publ.,
1994.
[3] A. Bogaerts, Z. Chen, R. Gijbels, A. Vertes Laser ablation for analytical sampling: what
can we learn from modeling?, Spectrochimica Acta Part B, 58 (2003), pp. 1867-1893.
[4] S. Chapman, T.G. Cowling, The mathematical theory of non-uniform gases, Cambridge University Press, Cambridge, 1953.
[5] A.J. Chorin, J.E. Marsden, A mathematical introduction to fluid mechanics, Springer-Verlag,
2000.
[6] W. E, B. Engquist, The heterogeneous multiscale methods, Comm. Math. Sci., 1 (2003),
pp. 87-132.
[7] C.W. Gear, T.J. Kaper, I.G. Kevrekidis, A. Zagaris, Projecting to a slow manifold: singularly perturbed systems and legacy codes, SIAM J. Appl. Dyn. Syst., 4 (2005), pp. 711-732.
[8] G.H. Golub, C.F. Van Loan, Matrix computations, The Johns Hopkins University Press,
1996.
[9] A.V. Gusarov, I. Smurov, Gas-dynamic boundary conditions of evaporation and condensation: numerical analysis of the Knudsen layer, Physics of Fluids, 14 (2002), pp. 4242-4255.
[10] S.M.J. Guzik, X. Gao, T. Weisgraber, B. Alder, P. Colella, An adaptive mesh refinement strategy with conservative space-time coupling for the lattice-Boltzmann method,
American Institute of Aeronautics and Astronautics, 2013.
[11] S.M. Guzik, T.H. Weisgraber, P. Colella, B.J. Alder, Interpolation methods and the
accuracy of lattice-Boltzmann mesh refinement, J. Comput. Phys., 259 (2014), pp. 461487.
[12] U.W. Hochstrasser, Orthogonal polynomials, In Handbook of mathematical functions: with
formulas, graphs, and mathematical tables, editors M. Abramowitz, I.A. Stegun, Dover
Publications, 1972.
[13] J. Huang, C. Yang, X.-C. Cai, Simulation of cavity flows by an implicit domain decomposition
algorithm for the discrete-velocity BGK model of the Boltzmann equation, Lect. Notes
Comput. Sci. Eng., 2014 (submitted).
[14] M. Junk, A. Klar, L. Luo, Asymptotic analysis of the lattice Boltzmann equation, J. Comput. Phys., 210 (2005), pp. 676-704.
[15] I.G. Kevrekidis, C.W. Gear, J.M. Hyman, P.G. Kevrekidis, O. Runborg, C. Theodoropoulos, Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level analysis, Commun. Math. Sci., 1 (2003), pp. 715-762.
[16] P. Lafitte, G. Samaey, Asymptotic-preserving projective integration schemes for kinetic equations in the diffusion limit, SIAM J. Sci. Comput., 34 (2012), pp. 579-602.
[17] L. Mieussens, Discrete velocity model and implicit scheme for the BGK equation of rarefied
gas dynamics, Math. Models Methods Appl. Sci., 10 (2000), pp. 1121-1149.
[18] H. Struchtrup, Macroscopic transport equations for rarefied gas flows: approximation methods in kinetic theory, Springer, 2005.
[19] S. Succi, The lattice Boltzmann equation for fluid dynamics and beyond, Oxford University
Press, Oxford, 2001.
[20] L.R. Turner, Inverse of the Vandermonde matrix with applications, Technical Report, NASA,
1966.
[21] R.G.M. Van der Sman, M.H. Ernst, Convection-diffusion lattice Boltzmann scheme for irregular lattices, J. Comput. Phys., 160 (2000), pp. 766782.
[22] P. Van Leemput, Multiscale and equation-free computing for lattice Boltzmann models, PhD
thesis, K.U. Leuven, 2007.
[23] P. Van Leemput, C. Vandekerckhove, W. Vanroose, D. Roose, Accuracy of hybrid
lattice Boltzmann/finite difference schemes for reaction-diffusion systems, Multiscale
Model. Simul., 6 (2007), pp. 838-857.
[24] P. Van Leemput, W. Vanroose, D. Roose, Mesoscale analysis of the equation-free Constrained Runs initialization scheme, Multiscale Model. Simul., 6 (2007), pp. 1234-1255.
24

[25] C. Vandekerckhove, Macroscopic simulation of multiscale systems within the equation-free
framework, PhD thesis, K.U. Leuven, 2008.
[26] C. Vandekerckhove, I. Kevrekidis, D. Roose, An efficient Newton–Krylov implementation
of the Constrained Runs scheme for initializing on a slow manifold, J. Sci. Comput., 39
(2009), pp. 167-188.
[27] Y. Vanderhoydonc, W. Vanroose, Lifting in hybrid lattice Boltzmann and PDE models,
Comput. Vis. Sci., 14 (2011), pp. 67-78.
[28] Y. Vanderhoydonc, W. Vanroose, Numerical extraction of a macroscopic PDE and a lifting
operator from a lattice Boltzmann model, Multiscale Model. Simul., 10 (2012), pp. 766-791.
[29] Y. Vanderhoydonc, W. Vanroose, Initialization of lattice Boltzmann models with the help
of the numerical Chapman–Enskog expansion, Procedia Computer Science, 18 (2013),
pp. 1036-1045.
[30] Y. Vanderhoydonc, W. Vanroose, C. Vandekerckhove, P. Van Leemput, D. Roose, Numerical lifting for lattice Boltzmann models. In Novel trends in lattice-Boltzmann methods,
volume 3 of Progress in Computational Physics, editor M. Ehrhardt, pp. 127-154, Bentham
Science Publishers, 2013.
[31] D.A. Wolf - Gladrow, Lattice-gas cellular automata and lattice Boltzmann models, Springer,
Berlin, 2000.

25

