Implicit schemes for wave models

Implicit schemes for wave models
Mathieu Dutour Sikirić
Rudjer Bo˘sković Institute, Croatia
and Universität Rostock
November 24, 2017

Stochastic wave modelling
Oceanic models are using grids (structured or unstructured) of
size 1km ≤ d ≤ 10km to simulate the ocean
But oceanic waves have a typical wavelength 2m ≤ L ≤
100m. So, we cannot resolve waves in the ocean.
But if one uses phase averaged models and uses stochastic
assumptions then it is possible to model waves by a spectral
wave action density N(x, k)
This density satisﬁes a Wave Action Equation (WAE) which
represents advection, refraction, frequency shifting and source
terms:
∂N
∂t
+ x ((cg + uA)N) + k( ˙kN) + θ( ˙θN) = Stot
with
Stot = Sin + Snl3 + Snl4 + Sbot + Sds + Sbreak + Sbf

The WWM model
The Wind Wave Model (WWM) is a unstructured grid
spectral wave model.
It is comparable to WaveWatch III, SWAN, WAM or SWAVE.
It incorporates most existing source term formulation for wind
input and dissipation (Cycle III, Cycle IV, Ardhuin, Makin, ...)
It has been coupled to SELFE, SHYFEM, TIMOR and ROMS.
It uses Residual Distribution schemes for the horizontal
advection.
It integrates the WAE by using the Operator Splitting Method
in explicit or implicit mode.

Operator Splitting Method
A standard technique for integrating partial diﬀerential
equations is the operator splitting method.
Over the interval [t0, t1] we successively solve the equations



∂N1
∂t + θ( ˙θN1) = 0 with N1(t0) = N(t0)
∂N2
∂t + k( ˙kN2) = 0 with N2(t0) = N1(t1)
∂N3
∂t + x ((cg + uA)N3) = 0 with N3(t0) = N2(t1)
∂N4
∂t = S(t) with N4(t0) = N3(t1)
and we set N(t1) = N4(t1).
No matter what the order of the successive integration
schemes is the ﬁnal order will be 1.
It it is possible to have higher order by more complex
integration procedures (Strang splitting, iterative splitting,
etc.)

The CFL criterion
If the discretization has characteristic length l and the
physical speed is c then we have the condition
c∆t
l
≤ 1
For the integration of the frequency and directional equations
we can subdivide the integration time step if necessary
because everything is decoupled.
This is not possible for the geographical advection:
The dependency in direction/frequency is small or negligible
The problem is that the group speed is
√
gh and so the CFL
number varies with the depth and the resolution.
So, we will present an implicit scheme for integrating N1, i.e.
in order to avoid the CFL limitation for advection.
Remark: the advection scheme used in implicit mode in
WWM is the residual distribution scheme PSI.

MPI parallelization I
The parallelization of geophysical models is usually done by
using the Mesage Passing Interface MPI formalism.
The set of computational nodes of the model is thus split into
a number of diﬀerent subdomains.
in MPI the exchanges are explicit. The explicit way of doing it
is via:
CALL MPI_SEND(ArrSend,len,dest,tag,comm,ierr)
CALL MPI_RECV(ArrRecv,len,orig,tag,comm,istat,ierr)
Those operations are blocking, i.e. the program waits until all
sends and recvs have been processed.
This means that all exchanges are processed by the order in
which they are stated.
It is generally better to decrease the total number of
exchanges in order to get better performance.

MPI parallelization II
In order to avoid strictly ordained exchanges, the strategy is
to do asynchrone exchanges.
The procedure is done in the following way
DO iorig=1,nproc
CALL MPI_IRECV(U,1,type(iorig), iorig-1, tag,
comm, rqst(iorig), ierr)
END DO
CALL MPI_WAITALL(nproc,rqst, stat, ierr)
and similarly for send operations.
The idea is the following: the array type(iorig) contains
the list of positions at which the received data needs to be
put. Commands for creating such types are for example
mpi create indexed block.
By using the above the order of the exchanges is no longer
determined by the MPI program which makes it faster but
harder to debug.

II. Iterative solution
methods

Iterative solution methods
In order to resolve linear system Ax = b for typical geophysical
situation we have matrices of size N × N with N about 100000
We cannot use direct methods like Gauss elimination or LU
and so we need to use iterative methods.
For a matrix A and a vector b the Krylov space Kn(A, b) is
Kn(A, b) = Vect b, Ab, . . . , An−1
b
The Generalized Minimal Residual Method (GMRES) takes
the best solution in Kn(A, b) of Ax = b.
It is stable but it requires the storing of n vectors, which is
memory intensive.
So, in order to have a good solution strategy we need a
method with minimal storage requirement.

The conjugate gradient method
If the matrix A is positive definite, then the conjugate
gradient method can be used:
« J.W. Shewchuk, An Introduction to the Conjugate Gradient
Method Without the Agonizing Pain Edition 11
4
If the system is N dimensional then N iteration suffices.
After i iterations, the residual error ei satisfies
ei ≤ 2
√
κ − 1
√
κ + 1
i
e 0 with κ =
λmax
λmin
Operations depends on computing Ax for some vectors x.
For non-symmetric problems, the technique is to use the
biconjugate gradient stabilized (BCGS) which works similarly.

Preconditioners
The convergence of the conjugate gradient depends on κ that
is on how far A is from the identity matrix.
If κ is large, i.e. A is ill conditioned then the number of
iterations will be very large.
We may accept that but then the whole solution strategy
becomes very similar to an explicit scheme.
The idea is to ﬁnd a matrix K for which we can compute the
inverse easily.
K must similar to A, i.e. share the same property as A.
In order to apply the BCGS we need to compute Ax and
K−1x for some vectors x.
Example: Jacobi preconditioning is to take the diagonal
entries of A.

Preconditioners for advection
The essential aspect of advection is that it moves things so
Jacobi preconditioner will not work.
Instead partial factorization techniques have to be used
We write A = D + E + F with D diagonal E lower triangular
and F upper triangular.
The Successive Over Relaxation (SOR) preconditioner is to say
A = (I + ED−1
)(D + F) + R with R = −ED−1
F
The incomplete LU factorization (ILU0) is to say
Aij = (LU)ij for Aij = 0
with L and U having the same sparsity as A.
Both methods are eﬃcient because they both are of the form
K = LU.
So, when solving Kx = b we do
x = L−1
b and x = U−1
x ,
i.e. the solution propagates.

Parallelizing solvers
Suppose we have to solve Lx = b and let us assume the
diagonal is 1.



x1 = b1
l2,1x1 + x2 = b2
...
...
...
lN,1x1 + · · · + lN,N−1xN−1 + xN = bN
So, we ﬁrst determine x1 then x2 and ﬁnally xN.
Parallelization is impossible if all Lij are non-zeros because
data from one processor
What save us is sparsity because the matrices are of the
following type:
f (x)v =
v ∼v
Cv,v xv
with v ∼ v mean that v and v are adjacent nodes.

Ordering nodes
All incomplete factorizations depend on the ordering of the
nodes.
We are free to choose the ordering that suits us best and by
doing so we change the preconditioner LU.
Since the iterative solution algorithms return approximate
solutions this means that the approximate solutions depend on
the partitioning and also on the number of processors.
The situation is the following:
Colored domains The L matrix

Coloring theory
A graph G is formed by a set V of vertices and a set E of
pairs of vertices named edges.
A coloring with N colors is a function f : V → {1, . . . , N}
such that for any edge e = (a, b) we have f (a) = f (b).
The chromatic number χ(G) is the minimum number of
colors needed to color.
It is known that χ(G) ≤ 4 for G a planar graph (Appel,
Haken, 1976).
Unfortunately, the subdomains given by parmetis are not
necessarily connected and so the graph is not necessarily
planar.
But in practice we can expect that the chromatic number is
rarely above 5.

Using colorings to solve Kx = b
Suppose that we managed to color with c colors
1. The indexing is done
1.1 First index the nodes in domains of color 1 by 1, 2, . . . , n1
1.2 Then the nodes of color 2 by n1 + 1, n1 + 2, . . . , n1 + n2.
1.3 .... until color c.
2. The solution of Lx = b is then done by
2.1 Solving Lx = b on the nodes of color 1.
2.2 Nodes of color 1 send data to nodes of higher color.
2.3 Solve Lx = b on the nodes of color 2.
2.4 Continue ...
3. The solution of Uy = x is then done in reverse.

Efficiency of preconditioners
There is no general theory on the efficiency of preconditioners.
The bad news is that the ordering of the nodes has an effect
on the performance of the preconditioner.
The worst ordering for κ is the red-back ordering in finite
difference schemes. The best is the linear ordering.
1 2
3 4
5 6
7 8
9 10
11 12
Red Black
1 7
11
3 9
4 10
2 8
6 12
5
Linear Ordering
So, the best ordering for the quality of the preconditioner is
the one that is hardest to parallelize.
The ordering that we used is somewhat intermediate. It is like
red-black globally, but over individual subdomains it is linear.

Organizing the computation
1. The penalty of parallelizing come in two ways:
1.1 The preconditioner quality that decreases.
1.2 The cost of waiting for data is c.
2. If we have Nfreq frequencies and Ndir directions then this
makes Ntot = Nfreq × Ndir independent linear systems to solve.
3. The strategy is then to split Ntot into b blocks B1, . . . , Bb
3.1 After domains of color 1 have ﬁnished block B1 data is sent
and block B2 is solved.
3.2 So domains of color 2 can start working before the ones of
color 1 are ﬁnished.
4. So, by using say b = 5 we can essentially remove the second
cost.

Further work
1. The work done so far is for the SOR preconditioner.
2. We need to test the ILU0 preconditioner, it is harder to
compute but the same strategy can be applied.
3. Another possibility is to integrate implicitly the advection in
geographical, frequency and direction.
This requires an ordering of the Nnode × Nfreq × Ndir matrix
entries but by doing so we can diminish the splitting error.
4. And overall improve the speed.
THANK YOU

Implicit schemes for wave models

Recommended

More Related Content

What's hot (20)

Similar to Implicit schemes for wave models (20)

More from Mathieu Dutour Sikiric (13)

Recently uploaded (20)

Implicit schemes for wave models