Simulating Viscous Fingering with a Timespace Method and Anisotropic Mesh
Adaptation
Kristian Ejlebjerg Jensen∗
Department of Earth Science and Engineering, Imperial College London, SW7 2AZ London, United Kingdom

arXiv:1508.03993v1 [] 17 Aug 2015

Abstract
We report findings related to a two dimensional viscous fingering problem solved with a timespace method and
anisotropic elements. Timespace methods have attracted interest for solution of time dependent partial differential
equations due to the implications of parallelism in the temporal dimension, but there are also attractive features in the
context of anisotropic mesh adaptation; not only are heuristics and interpolation errors avoided, but slanted elements
in timespace also correspond to long and accurate timesteps, i.e. the anisotropy in timespace can be exploited. We
show that our timespace method is restricted by a minimum timestep size, which is due to the growth of numerical
perturbations. The lower bound on the timestep is, however, quite high, which is indicative that the number of
timesteps can be reduced with several orders of magnitude for practical applications.
Keywords: anisotropic, mesh, adaptation, viscous, fingering, timespace

1. Introduction
Many engineering problems are governed by partial differential equations which can only be solved numerically.
The prospect of parallelising such simulations in the time dimension have attracted attention by the scientific community [1, 2, 3, 4], and the parareal algorithm [5] is popular and widely studied. The parallelism in the time dimension
can be achieved by leaving the distinction between time and space, in what is normally referred to as timespace
methods. The technique is often used for solving hyperbolic problems with dicontinuous Galerkin (DG) formulations
[6], as demonstrated for the full four dimensional case in the context of rotor flows [7]. For that application all the
complications of a time varying geometry disappear, when switching to a timespace method. This advantage was also
exploited for dynamic stall prediction as simulated with the Navier-Stokes equations [8]. That work utilised spacetime
adaptation and, as it is often the case, a compromise was made in the form of timeslabs, which allows for a tunable
problem size [9].
Anisotropic mesh adaptation is an established technique for ensuring computational efficiency in the context of
multiscale problems [10, 11, 12, 13]. There exists many heuristics methods for applying the technique to transient
problems [13, 14, 15], but the only way to exploit the anisotropy in the physics as seen in timespace, is by using a
timespace method. This is also a conceptually simple way to avoid the interpolation errors and heuristics associated
with combining conventional timestepping algorithms with mesh adaptation. The benefits of anisotropic mesh adaptation generally increases with dimensionality due to the fact that elements can be stretched in more dimensions, and
one could thus expect significant benefits from the use of a four dimensional anisotropic mesh adaptation algorithm,
but (to our knowledge) anisotropic mesh adaptation has only been demonstrated for two and three dimensions. We are
not aware of any studies involving the use of three dimensional anisotropic meshes for solving a timespace problem.
In this work, we focus on a miscible viscous fingering problem, which is transient and has two spatial dimensions.
The problem is known to exhibit chaotic behaviour [16], which causes some issues for the combination of a timespace
method and anisotropic mesh adaptation. We attribute this to the fact that the numerical noise can grow in time. This
∗ Corresponding

author
Email address: kristianejlebjerg@gmail.com (Kristian Ejlebjerg Jensen)

Preprint submitted to Journal of Computational Physics

January 14, 2018

could mean that a non-linear solver is unable to converge to the error level that one would expect for the computational
resources utilised. We suggest to use the method of timeslabs, where the slab thickness is a kind of timestep size. The
method is used to investigate the relation between the timeslab thickness, the computational resources used and growth
of numerical errors.
2. Anisotropic Mesh Adaptation
Mesh adaptation is the art of choosing the discretisation that minimises the error level for a given computational
cost. It is important to distinguish between p- and h-adaptation, which relate to the discretisation order and length
scale, respectively. Anisotropic mesh adaptation is the most general type of h-adaptation, because it not only considers
the local element size, but also the element orientation. The optimal size and orientation can be expressed in terms of
a metric tensor field [12], M. This has units of inverse square length, and it can be used to map elements to metric
space, where the optimal element has unit edge lengths. The metric that minimises the interpolation error [17] of a
variable with Hessian, H, is
M=

1
− 2q+d
1
det[abs(H)]
abs(H),
σ

(1)

where d is the number of dimensions, q is the error norm to be minimised, det is the determinant, σ is a scaling factor
and abs takes the absolute value in the principal frame. The metrics of several variables can be combined with the
inner ellipsoid method [13], illustrated for two dimensions in figure 1.

Figure 1: The inner ellipse method is illustrated in the case of intersection, but it is common to see one ellipse entirely within the other, such that
anisotropy is preserved.

We use 1st order polynomials for all fields, so we compute the Hessian by first performing a Galerkin projection
of the gradient onto 1st order polynomials and then repeating the process to get a nodal Hessian. The nodal metric
can then be calculated explicitly using equation (1).
We use the technique of local mesh modifications to adapt the mesh to the metric. There are four operation types:
Coarsening, swapping, refinement and smoothing, as shown in figure 2. Coarsening removes edges that are short in
metric space without regard to mesh quality, but the other modifications only allow the worst local element quality to
go up. The mesh quality is quantified by means of the Vassilevski functional [18].

Figure 2: Four local mesh modifications are illustrated: Coarsening, refinement, swapping and smoothing. Only the edge swapping operation is
illustrated in three dimensions, as it is trivial to generalise the other operations to three dimensions. Coarsening is the only operation that is allowed
to reduce the worst local element quality and therefore edges marked for refinement are not guaranteed to be split. We do not use face to edge
operations also known as 2-to-3 operations.

We use an Octave/MATLAB implementation, which is suboptimal in terms of performance, especially compared
to a similar C++ implementation [19].
2

3. Setup
We consider a simple two-phase miscible problem [16], where the viscosity, η, depends exponentially on the
saturation, φ, with φ = 0 and φ = 1 corresponding to the low and high viscosity fluid, respectively.
η(φ)

= η0 eξφ ,

(2)

where ξ determines the viscosity ratio. The velocity, v, is given by the Darcy equation,
v

= −

κ
∇p,
η(φ)

(3)

where p is the pressure and κ is the permeability. Mass conservation leads to a Poisson equation for the pressure,
!
κ
0 = ∇ · (ρv) = ∇ · ρ
∇p
η(φ)
!
κ
∇p ,
(4)
= ∇·
η(φ)
where ρ is the density. This drops out of the equations, because we choose to study the case of equal densities for the
two fluids. The saturation is convected with the velocity given in equation (3),
0

=
=

∂φ
+ v · ∇φ
∂t
∂φ
κ
−
∇p · ∇φ
∂t η(φ)

(5)

The equation is treated as a convective equation with a three dimensional velocity vector, and it is stabilised with
streamline upwind Petrov/Galerkin diffusion. The local characteristic length scale for the mesh is calculated by
computing the Steiner ellipsoid for every element and projecting this along the velocity direction.
The governing equations (4) and (5) are solved with the boundary conditions illustrated in figure 3. The use of a 13sided polygon to approximate a circle avoids extrapolation, when comparing solution fields on different meshes. The
perfect circular initial condition for the saturation makes the solution undefined, which causes the growth of numerical
perturbations to become very obvious. All equations are solved with the finite element method using FEniCS [20], an
open source finite element package. Both pressure and saturation is discretised with first order polynomials. A direct
solver (LU) is used for equations (4) and (5), while an iterative solver (CG+ILU) is used for Galerkin projections. All
computations are single threaded.
We non-dimensionalise the equations with ∆p as characteristic pressure and R as characteristic length scale,
x

=

Rx̃,

r

=

0.1R,

R2 η0
t˜, p = ∆p p̃,
ṽchar κ∆p
r1 = 0.2R, r2 = 0.3R, ξ = 2,

t=

t˜final = 1

where x̃, t˜, p̃ are dimensionless space, time and pressure, respectively. t˜final is the simulation time. ṽchar is a numerical
parameter, which determines the relative resolution of space and time. We fix it at unity. Another numerical parameters
is the norm of the interpolation error to be minimised, and we go with the popular choice of q = 2 for this.
The saturation at the initial time is imposed by means of a Dirichlet boundary condition. We initialize the nonlinear solver with φ = 0 and p̃ = p̃init , where
p̃init

r̃ − 0.1
,
0.9

=

r̃ being the radial coordinate. We use a segregated fixed point method to deal with the non-linearity of the equations. The timespace problem is solved in slabs with a thickness of ∆t (omitting the tilde), i.e.
#0 Initialise mesh, and set φ = 0, p̃ = p̃init .
3

Figure 3: The initial condition for the viscous fingering is sketched with the dark grey fluid pushing out on the light gray fluid, which has a higher
viscosity. The initial saturation is unperturbed and it is thus primarily the length scale of the initial mesh that triggers the instability. Zero pressure
is prescribed at the outer boundary, while p = ∆p and φ = 1 are prescribed at the inlet.

#1 Solve equation (4) for the pressure with fixed saturation.
#2 Solve equation (5) for the saturation by taking the viscosity to be constant using the previous saturation.
#3 If it is the first iteration at this timeslab, go to #1. Otherwise, we calculate metric fields associated with the
pressure and the saturation. These are combined with the inner ellipsoid method (see figure 1), the mesh is
updated and the fields are interpolated onto the new mesh. If it is not the first timeslab, the mesh is fixed at the
boundary matching up to the previous timeslab,
#4 If 20 iterations have been carried out with the current timeslab, we proceed to the next one by mirroring the
mesh in a plane normal to the time direction, and set φ = 0, p̃ = p̃init . We then go to #1.
The mesh is changed in every iteration, and therefore the non-linear solver cannot converge to machine precision. A
typical timeslab is shown in figure 4 together with the result of a full timespace simulation. Both figures illustrate the
value of full timespace anisotropy.
4. Results
We show result for timeslab thicknesses ∆t = 0.05, 0.1, 0.2 and 0.5 with σ = 0.01 as well as for mesh tolerances
σ = 0.04, 0.02 and 0.01 with ∆t = 0.1. The result at the end time is plotted in figures 5 and 6, respectively. The
final iteration number, time, ∆t value, σ value and final total computational time is printed above each image. As one
would expect, the shape and number of fingers clearly varies with σ, but not with ∆t.
We measure convergence in terms of the difference between saturation fields in consecutive iterations,
R
(φi − φi−1 )2 dΩ
L2i = Ω R
.
dΩ
Ω
This is plotted in figure 7(a-b). Figure 7(a) clearly shows that simulations become more accurate as ∆t is lowered,
but the convergence is different when σ is decreased, at least in the last half of the simulations: The error seems to
scale with σ, but it scales down to a certain level with ∆t. By comparing figure 7(a) and (b), we can see that this is the
level of the overall discretisation error. We attribute this effect to the fact that perturbations due to the mesh grow with
time due to the non-linear nature of the problem. The extreme case of this is illustrated in figure 8, where snapshots
of the φ = t˜ = 0.5 contour is shown for three different iterations with ∆t = 0.5. There are significant differences due
4

Figure 4: A simulation with σ = 0.02 and ∆t = 0.1 is plotted in terms of a wireframe of the final iteration with slab #5 (left) and all slabs in terms
of the φ = 0.5 isosurface (right). The time dimension is from left to right pointing out of the plane. Large elements clearly show up away from the
interface, but even at the interface the anisotropy in timespace can be exploited such that very few elements are needed. Note that the mesh is fully
tetrahedral; the quadrilaterals on the isosurface representation appear, when cutting tetrahedrons with two nodes on both sides of the threshold.

to the growth of mesh perturbations, particularly between iteration 16 and 20. This is to say, that the solution is not
well defined and therefore it has an oscillatory component arising from the changing mesh. The error associated with
this component increases with time, so it can be smaller than the overall discretisation error, if ∆t is sufficiently small.
Focusing on this case with σ = 0.01 and ξ = 2, it looks like ∆t = 0.1 is a good value.
5. Conclusion
The results show that it is possible to simulate chaotic phenomena with a timespace method. A CFL like condition
still seems to exists, so the timeslab cannot be too thick. This has to do with the fact that the chaotic nature of the
problem allows for a range of solutions and the difference between these solutions, should not exceed the overall
discretisation error. The lower limit for the timeslab thickness, is high enough to indicate that the total number of
timesteps can be reduced by orders of magnitude compared to a conventional method.
6. Suggestions for Future Work
It would be interesting to compare an established method against an optimised timespace implementation in the
context of large scale simulations and anisotropic mesh adaptation. If practical problems are to be solved, a four
dimensional anisotropic mesh adaptation algorithm will ultimately be needed.
7. Acknowledgement
This work is supported by the Villum Foundation.
8. References
References
[1] Yvon Maday and Gabriel Turinici. A parareal in time procedure for the control of partial differential equations. Comptes Rendus Mathematique, 335(4):387–392, 2002.

5

Figure 5: The saturation at the final time is plotted for four simulations with σ = 0.01 and ∆t= 0.05, 0.1, 0.2 and 0.05. The shape of the fingers
vary, but the number and length seems independent of ∆t.

Figure 6: The saturation at the final time is plotted for three simulations with ∆t = 0.1 and σ = 0.04, 0.02 and 0.01. Both the shape and number of
fingers seem highly dependent on σ.

[2] Guillaume Bal and Yvon Maday. A "parareal" time discretization for non-linear pde’s with application to the pricing of an american put. In
Recent developments in domain decomposition methods, pages 189–202. Springer, 2002.
[3] Pierre Gosselet and Christian Rey. Non-overlapping domain decomposition methods in structural mechanics. Archives of computational
methods in engineering, 13(4):515–572, 2006.
[4] Charbel Farhat and Marion Chandesris. Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and
fluid–structure applications. International Journal for Numerical Methods in Engineering, 58(9):1397–1434, 2003.
[5] Jacques-Louis Lions, Yvon Maday, and Gabriel Turinici. Résolution d’edp par un schéma en temps «pararéel». Comptes Rendus de
l’Académie des Sciences-Series I-Mathematics, 332(7):661–668, 2001.

6

Figure 7: The square root of the L2 residual is plotted as a function of normalised iterations numbers. (a) shows varying ∆t at σ = 0.01, while (b)
shows varying σ at ∆t = 0.1. The peaks occur whenever the timeslab is advanced in time.

Figure 8: The φ = t˜ = 0.5 contour is plotted for a simulation with ∆t = 0.5 and σ = 0.01 for iteration numbers 16, 18 and 20. Note how different
the fingers between iteration 16 and 20 are. We expect the mean radial position of the contours to be well defined, but not the actual fingers; they
are determined by the mesh at the initial time, which changes at every iteration.

[6] Shripad Thite. Efficient spacetime meshing with nonlocal cone constraints. arXiv preprint arXiv:0804.0946, 2008.
[7] Harmen van der Ven. An adaptive multitime multigrid algorithm for time-periodic flow simulations. Journal of Computational Physics, 227
(10):5286–5303, 2008.
[8] Christiaan M Klaij, Jaap JW van der Vegt, and Harmen van der Ven. Space–time discontinuous galerkin method for the compressible
navier–stokes equations. Journal of Computational Physics, 217(2):589–611, 2006.
[9] JJW Van der Vegt and H Van der Ven. Space–time discontinuous galerkin finite element method with dynamic grid motion for inviscid
compressible flows: I. general formulation. Journal of Computational Physics, 182(2):546–585, 2002.
[10] Adrien Loseille, Alain Dervieux, and Frédéric Alauzet. Fully anisotropic goal-oriented mesh adaptation for 3d steady euler equations. Journal
of computational physics, 229(8):2866–2897, 2010.
[11] Wagdi G Habashi, Julien Dompierre, Yves Bourgault, Djaffar Ait-Ali-Yahia, Michel Fortin, and Marie-Gabrielle Vallet. Anisotropic mesh
adaptation: Towards user-independent, mesh-independent and solver-independent cfd. part i: General principles. International Journal for
Numerical Methods in Fluids, 32(6):725–744, 2000.
[12] Adrien Loseille and Frédéric Alauzet. Continuous mesh framework part i: well-posed continuous interpolation error. SIAM Journal on
Numerical Analysis, 49(1):38–60, 2011.
[13] CC Pain, AP Umpleby, CRE De Oliveira, and AJH Goddard. Tetrahedral mesh optimisation and adaptivity for steady-state and transient
finite element calculations. Computer Methods in Applied Mechanics and Engineering, 190(29):3771–3796, 2001.
[14] Frédéric Alauzet, Pascal J Frey, Paul-Louis George, and Bijan Mohammadi. 3d transient fixed point mesh adaptation for time-dependent

7

problems: Application to cfd simulations. Journal of Computational Physics, 222(2):592–623, 2007.
[15] Xiangrong Li, Mark S Shephard, and Mark W Beall. 3d anisotropic mesh adaptation by mesh modification. Computer methods in applied
mechanics and engineering, 194(48):4915–4950, 2005.
[16] Satyajit Pramanik, GL Kulukuru, and Manoranjan Mishra. Miscible viscous fingering: Application in chromatographic columns and aquifers.
In COMSOL conference, Bangalore, 2012.
[17] Long Chen, Pengtao Sun, and Jinchao Xu. Optimal anisotropic meshes for minimizing interpolation errors in l p -norm. Mathematics of
Computation, 76(257):179–204, 2007.
[18] Yu V Vasilevski and KN Lipnikov. Error bounds for controllable adaptive algorithms based on a hessian recovery. Computational Mathematics
and Mathematical Physics, 45(8):1374–1384, 2005.
[19] Georgios Rokos, Gerard J Gorman, James Southern, and Paul HJ Kelly. A thread-parallel algorithm for anisotropic mesh adaptation. arXiv
preprint arXiv:1308.2480, 2013.
[20] Anders Logg, Kent-Andre Mardal, Garth N. Wells, et al. Automated Solution of Differential Equations by the Finite Element Method.
Springer, 2012. ISBN 978-3-642-23098-1. doi: 10.1007/978-3-642-23099-8.

8

