Élie Bretin

Roland Denis

Jacques-Olivier Lachaud

Édouard Oudet

Phase-field modelling and computing for a large number of phases

We propose a framework to represent a partition that evolves under mean curvature flows and volume constraints. Its principle follows a phase-field representation for each region of the partition, as well as classical Allen-Cahn equations for its evolution. We focus on the evolution and on the optimization of problems involving high resolution data with many regions in the partition.

In this context, standard phase-field approaches require a lot of memory (one image per region) and computation timings increase at least as fast as the number of regions. We propose a more efficient storage strategy with a dedicated multi-image representation that retains only significant phase field values at each discretization point. We show that this strategy alone is unfortunately inefficient with classical phase field models. This is due to non local terms and low convergence rate. We therefore introduce and analyze an improved phase field model that localizes each phase field around its associated region, and which fully benefits of our storage strategy.

To demonstrate the efficiency of the new multiphase field framework, we apply it to the famous 3D honeycomb problem and the conjecture of Weaire-Phelan's tiling.

Introduction

This paper is concerned by the numerical approximation of an evolving Q-partition Ω(t) following a normal velocity law given by a conserved multi-mean curvature flow. Recall that a Q-partition Ω = {Ω 1 , Ω 2 , . . . , Ω N } satisfies for i = j,

∪ N i=1 Ω i = Q, and |Ω • i ∩ Ω • j | = 0, ∀(i, j) ∈ {1, 2, . . . , N } 2 .
The evolution of the partition Ω can be obtained as the L 2 -gradient flow (with additional volume constraint) of a multiphase perimeter

P(Ω) = N i,j=1 ˆΓi,j 1dσ = N i=1 ˆ∂Ωi 1dσ,
where Γ i,j = ∂Ω i ∩ ∂Ω j . In particular, the classical physical theory [START_REF] Herring | Surface Tension as a Motivation for Sintering[END_REF][START_REF] Mullins | Two-Dimensional Motion of Idealized Grain Boundaries[END_REF] states that this evolution must follow at least two rules :

(1) At every point x ∈ Γ i,j which is not a junction point between three of more interfaces, the normal velocity V i,j of the interface Γ i,j is proportional to its mean curvature H i,j :

V i,j (x) = σ i,j H i,j (x) a.e. x ∈ Γ i,j , where σ i,j correspond to the surface tension coefficients. We assume in this study that all σ i,j are equal to one.

(2) The Herring's angle condition holds at every triple junction, e.g. if x is a junction point between phases i, j and k then σ i,j n i,j + σ j,k n j,k + σ k,i n k,i = 0, where n i,j denotes the unit normal at x to Γ i,j , pointing from Ω i to Ω j .

Additionally, we assume that the volume of each phase Ω i is conserved during the evolution, i.e.

d dt

Vol(Ω i ) = d dt ˆΩi 1dx = 0, for all i ∈ {1, 2, . . . , N }.

There is a vast literature on numerical methods for the approximation of mean curvature flows. The methods can be roughly classified into the following four categories (some of them are exhaustively reviewed and compared in [START_REF] Deckelnick | Computation of geometric partial differential equations and mean curvature flow[END_REF]): parametric methods [START_REF] Deckelnick | Computation of geometric partial differential equations and mean curvature flow[END_REF][START_REF] Barrett | On the parametric finite element approximation of evolving hypersurfaces in R 3[END_REF], level set methods [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF][START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF][START_REF] Osher | Geometric Level Set Methods in Imaging, Vision and Graphics[END_REF][START_REF] Evans | Motion of level sets by mean curvature[END_REF][START_REF] Chen | Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations[END_REF], convolution/thresholding BMO type algorithms [START_REF] Bence | Diffusion generated motion by mean curvature[END_REF][START_REF] Ishii | Threshold dynamics type approximation schemes for propagating fronts[END_REF][START_REF] Ruuth | Efficient algorithms for diffusion-generated motion by mean curvature[END_REF][START_REF] Esedoglu | Threshold dynamics for networks with arbitrary surface tensions[END_REF] and phase field approximations [START_REF] Modica | Un esempio di Γ--convergenza[END_REF][START_REF] Chen | Generation and propagation of interfaces for reaction-diffusion equations[END_REF].

The main goal of this work is to propose an efficient and accurate numerical phase field model [START_REF] Modica | Un esempio di Γ--convergenza[END_REF][START_REF] Chen | Generation and propagation of interfaces for reaction-diffusion equations[END_REF] of conserved multi-mean curvature flow in the case of a great number of phases N . More precisely, we will see that classical approaches lead to some instability issues which can be explained by using standard formal asymptotic expansion [START_REF] Caginalp | Dynamics of layered interfaces arising from phase boundaries[END_REF][START_REF] Pego | Front migration in the nonlinear Cahn-Hilliard equation[END_REF][START_REF] Bellettini | Quasi-optimal error estimates for the mean curvature flow with a forcing term[END_REF][START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF][START_REF] Garcke | On anisotropic order parameter models for multi-phase systems and their sharp interface limits[END_REF]. In particular, we will see that the classical phase field method has a precision of order 1 only. We will improve the accuracy of this numerical model by proposing a new version of order 2.

Classical phase field models

We first introduce the special case N = 2. Let us consider an evolution set Ω(t) with a normal velocity law given by V n = H.

As previously, H denotes the mean curvature of the interface and this evolution can also be obtained as the L 2 -gradient flow of the perimeter

P (Ω) = ˆ∂Ω 1dσ.
The concept of the phase field method consists to approximate the perimeter by the Cahn Hilliard's energy [START_REF] Cahn | Free energy of a nonuniform system. I. Interfacial free energy[END_REF], defined for all function u ∈ H 1 (Q, R), by

P ε (u) = ˆQ ε |∇u| 2 2 + 1 ε W (u) dx, (1)
where ε is a small approximation parameter and W is a double well potential assumed to be W (s) = 1 2 s 2 (1 -s) 2 in this work. Notice that Modica and Mortola [START_REF] Modica | Un esempio di Γ--convergenza[END_REF] proved that P ε converges to c W P in the sense of the Γ-convergence for the L 1 topology, where the constant c W = ´1 0 2W (s)ds and where

Here, |Du|(Ω) denotes the total variation of a function u ∈ BV (Ω) (see [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF]) and is defined by

|Du|(Ω) = ˆΩ |Du| = sup ˆΩ u div gdx; g ∈ C 1 0 (Ω, R n ), |g| ≤ 1 .
In particular, when Ω is a set with finite perimeter [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF][START_REF] Maggi | Sets of finite perimeter and geometric variational problems: an introduction to geometric measure theory[END_REF], the characteristic function 1 Ω of Ω satisfies P (1 Ω ) = P (Ω). More precisely, Modica and Mortola show that 1 Ω can be approximated by the sequence u ε (x) = q dist(x,Ω) ε which satisfies

lim ε→0 P ε (u ε ) = c W P (Ω).
Here, dist(x, Ω) is the signed distance function to Ω and q is the profile function associated to the potential W which is defined by q = argmin p ˆR W (p(s))|p (s)|ds; p(-∞) = 1, p(0) = 1/2, p(+∞) = 0 , where p ranges over all Lipschitz continuous functions p : R → R. It is a well-known fact that q (s) = -2W (q(s)) and q (s) = W (q(s)), for all s ∈ R, which implies that q(s) = (1 -tanh(s))/2 in the case of the standard double well potential W (s) =

1 2 s 2 (1 -s) 2 .
The classical Allen Cahn equation [START_REF] Allen | A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening[END_REF], obtained as the L 2 -gradient flow of P ε , reads

u t = ∆u - 1 ε 2 W (u).
Note that existence, uniqueness, and a comparison principle have been established for this equation (see for example chapters 14 and 15 in [START_REF] Ambrosio | Geometric evolution problems, distance function and viscosity solutions[END_REF]).

A smooth mean curvature flow Ω(t) can then be approximated by

Ω ε (t) = x ∈ R d ; u ε (x, t) ≥ 1 2
where u ε is the solution of the Allen Cahn equation with the following initial condition

u ε (x, 0) = q dist(x, Ω(0)) ε .
Notice that standard formal asymptotic expansion [START_REF] Caginalp | Dynamics of layered interfaces arising from phase boundaries[END_REF][START_REF] Pego | Front migration in the nonlinear Cahn-Hilliard equation[END_REF] of u ε near the interfaces shows that u ε is quadratically closed to the optimal profile, i.e

u ε (x, t) = q dist(x, Ω ε ) ε + O(ε 2 ),
with associated normal velocity V ε satisfying

V ε = H + O(ε 2 ).
More rigorously, the convergence of ∂Ω ε (t) to ∂Ω(t) has been proved for smooth motions by [START_REF] Chen | Generation and propagation of interfaces for reaction-diffusion equations[END_REF][START_REF] De Mottoni | Geometrical evolution of developed interfaces[END_REF][START_REF] Bellettini | Quasi-optimal error estimates for the mean curvature flow with a forcing term[END_REF] with an optimal convergence rate about O(ε 2 | log ε| 2 ).

The case of multiphase system N > 2 is very similar. The multi-perimeter P can be approximated by the multi Cahn-Hilliard energy defined for all u = (u 1 , u 2 , . . . , u N ) by

P ε (u) = 1 2 N i=1 ´Q ε |∇ui| 2 2 + 1 ε W (u i ) dx, if N i=1 u i = 1, +∞ otherwise. Notice that if u ε is a sequence satisfying lim inf P ε (u ε ) < ∞, then Modica-Mortola's result proves that u ε → u = (1 Ω1 , 1 Ω2 , . . . , 1 Ω N ), and the constraint N i=1 u i = 1 ensures that Ω = {Ω 1 , Ω 2 , . . . , Ω N } is a partition of Q.
The Γ-convergence of P ε to c W P is also established in [START_REF] Oudet | Approximation of partitions of least perimeter by Gamma-convergence: around Kelvin's conjecture[END_REF]. More general Γ-convergence results in the case of non uniform surface tension σ i,j is obtained in [START_REF] Baldo | Minimal interface criterion for phase transitions in mixtures of cahn-hilliard fluids[END_REF][START_REF] Bretin | A new phase field model for inhomogeneous minimal partitions, and applications to droplets dynamics[END_REF]. Note also that phase field models that handle the case of anisotropic surface tension are introduced and analyzed in [START_REF] Garcke | A multi phase field concept: Numerical simulations of moving phase boundaries and multiple junctions[END_REF][START_REF] Garcke | On anisotropic order parameter models for multi-phase systems and their sharp interface limits[END_REF].

The L 2 -gradient flow of P ε translates into the following Allen Cahn system

∂ t u i (x, t) = ∆u i (x, t) - 1 ε 2 W (u i ) + λ(x, t), ∀i = 1, . . . , N ,
where the Lagrange multiplier λ is defined by the partition constraints N j=1 u j = 1 and is explicitly computed by the formula

λ = 1 N N i=1 1 ε W (u i ).
More precisely, sharp interface limit derived for instance in [START_REF] Garcke | On anisotropic order parameter models for multi-phase systems and their sharp interface limits[END_REF], shows that, at least formally, the solution u ε expands near the interface Γ i,j as

       u ε i = q dist(x,Ω ε i ) ε + O(ε 2 ), u ε j = 1 -q dist(x,Ω ε i ) ε + O(ε 2 ), u ε k = O(ε 2 ), for k ∈ {1, 2, . . . , N } \ {i, j},
with associated normal velocity law V ε i,j satisfying

V ε i,j = H i,j + O(ε 2 ).
Then we obtained that this phase model if formally of order two with respect to ε.

The case with additional volume constraints is similar. Having in mind the approximation Vol(Ω i ) ´Q u i dx, the volume constraint on phase Ω i can be approximated on phase field function u i by

d dt ˆQ u i dx = 0.
Then, the constrained-L 2 -gradient flow of P ε reads as

∂ t u i (x, t) = ∆u i (x, t) - 1 ε 2 W (u i ) + λ(x, t) + µ i (t) 2W (u i ), for all i ∈ {1, 2, . . . , N }, (2)
where µ i 2W (u i ) is the Lagrange multiplier associated to the volume constraint d dt ´Q u i dx = 0.

Remark 1.1. Although it is more natural to use a Lagrange multiplier on the form µ i instead of µ i 2W (u i ), our motivation follows the case of conserved mean curvature flow where it is proved [START_REF] Brassel | A modified phase field approximation for mean curvature flow with conservation of the volume[END_REF][START_REF] Alfaro | Convergence of a mass conserving allen-cahn equation whose lagrange multiplier is nonlocal and local[END_REF] that using a Lagrange multiplier on the form 2W (u)µ(t) instead of µ(t) leads to a conservation of the volume with an error of order O(ε 2 ) instead of O(ε).

Moreover, in the case of periodic or Neumann boundary conditions on ∂Q, all these Lagrange multipliers become explicit and can be identified as

µ i = 1 ´Ω 2W (u i )dx 1 ε 2 ˆQ W (u i )dx -Λ , and
λ = 1 N N i=1 1 ε 2 W (u i ) - ´Q W (u i )dx ´Q 2W (u i )dx 2W (u i ) + 1 N N i=1 2W (u i ) ´Ω 2W (u i )dx Λ.
Here Λ = ´Q λdx is a free parameter that is the consequence of the non independence of the constraints that regions do not mix and that regions have fixed volume. In practice, we do not observe any influence on the choice of Λ with respect to the stability and the efficiency in our numerical experiments. We then take Λ = 0.

Computational complexity and limitation of standard phase field models

Recall that we have in mind applications which involve hundreds of phases on a three-dimensional spatial domain discretized with 1024 3 points. In that case, storing one phase field function u i (in double-precision floating-point format) in all the computation box requires 8 GiB of memory. Simulating the evolution of partition with more than a few phases becomes nearly impossible on an average computation server. We then need to use an adapted storage strategy. As pointed previously by many works (see e.g. [START_REF] Gruber | Sparse data structure and algorithm for the phase field method[END_REF], [START_REF] Vedantam | Efficient numerical algorithm for multiphase field simulations[END_REF] and [START_REF] Vanherpe | Bounding box algorithm for threedimensional phase-field simulations of microstructural evolution in polycrystalline materials[END_REF]), the phase-field function u i is expected of the form

u i (x) q dist(x, Ω i ) ε 1 Ωi (x), (3)
and vanishes far away from its support set Ω i . It appears then reasonable to save memory and computational time outside of this localized area.

Unfortunately, the previous classical phase field model leads to an insufficiently accurate approximation (3), which does not vanish enough far away from the interfaces. Therefore, a storage strategy based on storing non-null parts is not effective with this model. As an illustration, we plot on figure 1 a numerical experiment obtained on a two-dimensional domain with nine phases. We can then observe that the phase u 1 can be negative and that a non negligible part of its mass is located on every interface and every multiple junction. More precisely, the magnitude of this badly localized part is about 10 -3 . It corresponds to the same order as ε, that is equals to 1/256 in this numerical experiment. In particular, it suggests a phase field model of order 1 only. Other phases show the same behavior. In order to find the origin of this bad localization of the mass, we derive in section 2 an asymptotic expansion of the solution u ε . In particular, we observe that near an interface Γ i,j , u ε is in fact expected of the form

       u ε i = q dist(x,Ω ε i ) ε + O(ε), u ε j = 1 -q dist(x,Ω ε i ) ε + O(ε), u ε k = O(ε), for k ∈ {1, 2, . . . , N } \ {i, j}.
A direct consequence of this result is that a quantity of size O(ε) of the mass of u i may be localized outside of the set Ω i , which raises some instability issues when the number of the phase becomes large.

A new and more accurate phase field model

We shall now improve the accuracy of the previous phase field model. The idea is to modify the form of the Lagrange multiplier λ in order to co-localize its effect with the interface of the current phase-field function. Following the same idea as [START_REF] Brassel | A modified phase field approximation for mean curvature flow with conservation of the volume[END_REF], we propose to replace λ by λ 2W (u i ), which leads to this new slightly modified phase field model:

∂ t u i (x, t) = ∆u i (x, t) - 1 ε 2 W (u i ) + λ(x, t) 2W (u i ) + µ i (t) 2W (u i ), (4)
where λ and µ i are defined accordingly to the following constraints :

d dt ˆQ u i dx = 0, N j=1 u j = 1.
As previously, straightforward computations show that λ and µ i can be obtained as

µ i = 1 ε 2 ´Q W (u i )dx ´Q 2W (u i )dx - λ i ´Q 2W (u i )dx (5)
and

λ = N j=1 1 ε 2 W (u j ) -µ j 2W (u j ) N j=1 2W (u j ) . (6)
Here λ i = ´Q λ(x, t) 2W (u i )dx and the vector λ = (λ 1 , λ 2 , . . . , λ N ) can be computed as the solution of the linear system (Id -A)λ = b, where

A i,j = ˆQ 2W (u i ) 2W (u j )/ ´Q 2W (u j )dx N k=1 2W (u k ) dx and b i = ˆQ     1 ε 2 N k=1 W (u k ) - ´Q W (u k )dx ´Q √ 2W (u k )dx 2W (u k ) N k=1 2W (u k ) 2W (u i )     dx.
Notice that a numerical experiment similar to figure 1 is plotted in figure 2 using this new slight variant of the classical phase field model. We now clearly observe that the mass of each phase u i is well localized in the set Ω i . Moreover, the magnitude of the negative part is now about 10 -6 which corresponds to the size of ε 2 . Thus, this model presents better properties of localization which is of crucial interest from a practical point of view as it is detailed in the next sections. To estimate the order of precision of this phase field model, we also derive an asymptotic expansion of its solution. In particular, we will see, at least formally, that u ε is now expected near interface Γ i,j to be of the form

       u ε i = q dist(x,Ω ε i ) ε + O(ε 2 ), u ε j = 1 -q dist(x,Ω ε i ) ε + O(ε 2 ), u ε k = O(ε 2 ), for k ∈ {1, 2, . . . , N } \ {i, j}.
In particular, this experiment confirms the previous observations: this new model induces a quantity of size O(ε 2 ) only of the mass of u i that may lie outside of the set Ω i . With this new model, the memory storage strategy allows us to perform numerical phase field experiment with more than hundred phases in 3D and with a resolution grid of size 1024 3 on an average calculation server as it is illustrated below.

Outline of the paper

The section 2 is devoted to the formal asymptotic expansion of the two phase field models introduced above. Section 3 presents the numerical scheme while Section 4 details its implementation (compressed storage technique, algorithm details, interfaces extraction and refinement). Finally, the section 5 presents numerical experiments, where some examples in 2D and 3D involving more than a thousand of phases are computed. The classical 2D and 3D honeycomb problem is addressed numerically as a testbed for our multiphase model.

Formal asymptotic expansion of the two phase field model

This section concerns the asymptotic expansion of the solution of the two phase field models :

• The classical one

∂ t u i (x, t) = ∆u i (x, t) - 1 ε 2 W (u i ) + λ(x, t) + µ i (t) 2W (u i ).
• The modified approach

∂ t u i (x, t) = ∆u i (x, t) - 1 ε 2 W (u i ) + λ(x, t) 2W (u i ) + µ i (t) 2W (u i ).
We follow the method of matched asymptotic expansions proposed in [START_REF] Caginalp | Dynamics of layered interfaces arising from phase boundaries[END_REF][START_REF] Pego | Front migration in the nonlinear Cahn-Hilliard equation[END_REF][START_REF] Bellettini | Quasi-optimal error estimates for the mean curvature flow with a forcing term[END_REF][START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF] around the interface Γ i,j .

Outer expansion: we assume that the so-called outer expansion of u ε i , i.e. the expansion far from the front Γ i,j , is of the form

u ε i (x, t) = u 0 i (x, t) + εu 1 i (x, t) + O(ε 2 ), for all i ∈ {1, 2, . . . , N }.
In particular, analogously to [START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF], it is not difficult to see that

u 0 i (x, t) = 1 if x ∈ Ω i (t) 0 otherwise , u 0 j (x, t) = 0 if x ∈ Ω i (t) 1 otherwise
and

u 1 i = u 1 j = u 0 k = u 1 k = 0, for all k ∈ {1, 2, . . . , N } \ {i, j}.
Inner expansion: in a small neighborhood of Γ i,j , we define the stretched normal distance to the front,

z = d i (x, t) ε ,
were d i (x, t) denotes the signed distance to Ω i (t) = {x ∈ Q, u ε i (x, t) ≥ 1 2 } such that d i (x, t) < 0 in Ω i (t).
We then focus on inner expansions of u ε i (x, t), λ ε (x, t) i.e. expansions close to the front, of the form

u ε i (x, t) = U ε i (z, x, t) = U 0 i (z, x, t) + εU 1 i (z, x, t) + O(ε 2 ), for all i ∈ {1, 2, . . . , N }.
Additionally, we consider inner expansion of λ ε

λ ε (x, t) = Λ ε (z, x, t) = ε -2 Λ -2 (z, x, t) + ε -1 Λ -1 (z, x, t) + 0(1).
Moreover, we make the Ansatz that

µ ε i (t) = ε -1 M -1 i (t) + M 0 i + εM 1 i + O(ε 2 ).
Let us also define a unit normal m i to Γ i,j and the normal velocity V ε i,j to the front as

V ε i,j = -∂ t d i (x, t) = V 0 i,j + εV 1 i,j + O(ε 2 ), m i = ∇d i (x, t), x ∈ Γ i,j ,
where ∇ refers to spatial derivation only (the same holds for further derivation operators used in the sequel). Following [START_REF] Pego | Front migration in the nonlinear Cahn-Hilliard equation[END_REF][START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF] we assume that U ε i (z, x, t) does not change when x varies in the normal direction of Γ i,j with z held fixed, or equivalently ∇ x U ε i • m = 0. This is equivalent to requiring that the blow-up with respect to the parameter ε is coherent with the flow. Following [START_REF] Pego | Front migration in the nonlinear Cahn-Hilliard equation[END_REF][START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF], it is easily seen that     

∇u ε i = ∇ x U ε i + ε -1 m i ∂ z U ε i , ∆u ε i = ∆ x U ε i + ε -1 ∆d i ∂ z U ε i + ε -2 ∂ 2 zz U ε i , ∂ t u ε i = ∂ t U ε i -ε -1 V ε i,j ∂ z U ε i .
Also recall that in a sufficiently small neighborhood of Γ i,j , according to Lemma 14.17 in [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF], we have

∆d i (x, t) = n-1 k=1 κ k (π(x)) 1 + κ k (π(x))d i (x, t) = n-1 k=1 κ k (π(x)) 1 + κ k (π(x))εz
where π(x) is the projection of x on Γ i,j , and κ k are the principal curvatures on Γ i,j . In particular this implies that ∆d

i (x, t) = H i,j -εz A i,j 2 + O(ε 2 ),
where H i,j and A i,j 2 denote, respectively, the mean curvature and the squared 2-norm of the second fundamental form on Γ i,j at π(x).

Matching conditions: the matching conditions which connect outer and inner expansion (see [START_REF] Loreti | Propagation of fronts in a nonlinear fourth order equation[END_REF] for more details) imply in particular that

lim z→+∞ U 0 i (z, x, t) = 0, lim z→-∞ U 0 i (z, x, t) = 1, lim z→±∞ U 1 i (z, x, t) = 0, lim z→+∞ U 0 j (z, x, t) = 1, lim z→-∞ U 0 j (z, x, t) = 0, lim z→±∞ U 1 j (z, x, t) = 0,
and lim

z→±∞ U 0 k (z, x, t) = lim z→±∞ U 1 k (z, x, t) = 0, for all k ∈ {1, 2, . . . , N } \ {i, j}.

Classical model and inner expansion

We first consider the phase field model

∂ t u i (x, t) = ∆u i (x, t) - 1 ε 2 W (u i ) + λ(x, t) + µ i 2W (u i ).
The first order in ε -2 reads for all k ∈ {1, 2, . . . , N } as

∂ 2 zz U 0 k -W (U 0 k ) + Λ -2 = 0.
Notice that by adding the boundary conditions obtained from the matching conditions, and using U 0 i (0, x, t) = 1/2, we have

U 0 i (z, x, t) = q(z), U 0 j (z, x, t) = q(-z), U 0 k (z, x, t) = 0 and Λ -2 = 0.
Moreover, the second order in ε -1 shows that for all k ∈ {1, 2, . . . , N },

V i,j ∂ z U 0 k = ∂ 2 zz U 1 k -W (U 0 k )U 1 k + H i,j ∂ z U 0 k + M -1 k 2W (U 0 k ) + Λ -1 .
In particular, the cases k = i and k = j read respectively

V 0 i,j q (z) = ∂ 2 zz U 1 i -W (q(z))U 1 i + q (z)H i,j -M -1 i q (z) + Λ -1 ,
and

-V 0 i,j q (z) = ∂ 2 zz U 1 j -W (q(z))U 1 j -q (z)H i,j -M -1 j q (z) + Λ -1 . It follows that 2V 0 i,j q (z) = ∂ 2 zz U 1 i -W (q(z))U 1 i -∂ 2 zz U 1 j -W (q(z))U 1 j + 2q (z)H i,j -(M -1 i -M -1 j )q (z).
Multiplying the last equation by q and integrating over R leads to

V 0 i,j = H i,j - 1 2 (M -1 i -M -1 j ). Indeed, remark that ˆR ∂ 2 zz U 1 i -W (q(z))U 1 i q (z)dz = 0, as ˆR W (q(z))U 1 i q (z)dz = ˆR (W (q(z)) U 1 i dz = -ˆR W (q(z))∂ z (U 1 i )dz, = -ˆR q (z)∂ z (U 1 i )dz = ˆR q (z)∂ 2 zz (U 1 i )dz,
and where we use the matching conditions : lim z→±∞ U 1 j (z, x, t) = 0.

Moreover, the Lagrange multiplier Λ -1 satisfies

ˆR Λ -1 (z, x, t)q (z)dz = c W M -1 j + M -1 i 2 ,
where c W = ´R(q (s)) 2 ds = ´1 0 2W (s)ds. This shows that there exists a profile η such as Λ -1 equals to

Λ -1 (z, x, t) = 1 2 (M -1 j + M -1 i ) η(z).
We then obtain that each U 1 k can be defined as the solution of the system

     ∂ 2 zz U 1 i -W (q(z))U 1 i = 1 2 (M -1 j + M -1 i ) (q (z) -η(z)) , ∂ 2 zz U 1 j -W (q(z))U 1 j = 1 2 (M -1 j + M -1 i ) (q (z) -η(z)) , ∂ 2 zz U 1 k -W (U 1 k ) = 1 2 (M -1 j + M -1 i )η(z), if k ∈ {1, 2, . . . , N } \ {i, j}.
with additional Dirichlet boundary conditions at z → ±∞. Notice also that we have the additional constraint

N k=1 U 1 k = 0, which shows that U 1 k cannot vanished for all k as soon as M -1 j + M -1 i = 0.
In particular, on the neighborhood of the interface Γ i,j , the solution u of the phase field model is expected of the form

       u ε i (x, t) = q dist(x,Ωi) ε + O(ε), u ε j (x, t) = 1 -q dist(x,Ωi) ε + O(ε), u ε k (x, t) = +O(ε), for all k ∈ {1, 2, . . . , N } \ {i, j}, if M -1 j + M -1 i = 0, and        u ε i (x, t) = q dist(x,Ωi) ε + O(ε 2 ), u ε j (x, t) = 1 -q dist(x,Ωi) ε + O(ε 2 ), u ε k (x, t) = +O(ε 2 ), for all k ∈ {1, 2, . . . , N } \ {i, j}, if M -1 j + M -1 i = 0.
Finally, the classical phase field model admits an order of convergence about O(ε) only.

Modified model and inner expansion

We now consider the phase field model

∂ t u i (x, t) = ∆u i (x, t) - 1 ε 2 W (u i ) + λ(x, t) 2W (u i ) + µ i (t) 2W (u i ).
As previously, the first order in ε -2 reads for all k ∈ {1, 2, . . . , N } as

∂ 2 zz U 0 k -W (U 0 k ) + Λ -2 = 0,
which leads to (with additional boundary conditions)

U 0 i (z, x, t) = q(z), U 0 j (z, x, t) = q(-z), U 0 k (z, x, t) = 0 and Λ -2 = 0.
Moreover, the second order in ε -1 shows that for all k ∈ {1, 2, . . . , N },

V i,j ∂ z U 0 k = ∂ 2 zz U 1 k -W (U 0 k )U 1 k + H i,j ∂ z U 0 k + M -1 k 2W (U 0 k ) + Λ -1 2W (U 0 k ).
In particular, using

V 0 i,j q (z) = ∂ 2 zz U 1 i -W (q(z))U 1 i + q (z)H i,j -M -1 i q (z) -Λ -1 q (z) and -V 0 i,j q (z) = ∂ 2 zz U 1 j -W (q(z))U 1 j -q (z)H i,j -M -1 j q (z) -Λ -1 q (z) implies the following equation: 2V 0 i,j q (z) = ∂ 2 zz U 1 i -W (q(z))U 1 i -∂ 2 zz U 1 j -W (q(z))U 1 j + 2q (z)H i,j -(M -1 i -M -1 j )q (z).
Multiplying the last equation by q and integrating over R shows that

V 0 i,j = H i,j - 1 2 (M -1 i -M -1 j ). The Lagrange multiplier Λ -1 satisfies ˆR Λ -1 (z, x, t)q (z)dz = c W (M -1 j + M -1 i )/2,
and is still expected on the form

Λ -1 (z, x, t) = 1 2 (M -1 j + M -1 i ) η(z).
The functions U 1 k are then defined as solutions of the system

     ∂ 2 zz U 1 i -W (q(z))U 1 i = 1 2 (M -1 j + M -1 i ) (q (z) -η(z)) , ∂ 2 zz U 1 j -W (q(z))U 1 j = 1 2 (M -1 j + M -1 i ) (q (z) -η(z)) , ∂ 2 zz U 1 k -W (U 1 k ) = 0, if k ∈ {1, 2, . . . , N } \ {i, j}.
We can now see that U 1 k = 0 for all k ∈ {1, 2, . . . , N } \ {i, j}, and this equality holds even if

M -1 j + M -1 i = 0. Moreover, as N k=1 U 1 k = 0, we can also deduce that η(z) = q (z), U 1 i = 0 and U 1 j = 0.
Finally, the Lagrange multiplier is identified to

Λ -1 = 1 2 (M -1 j + M -1 i ) q (z)
and, on the neighborhood of the interface Γ i,j , the solution u ε of the new proposed phase field model is now expected on the form

       u ε i (x, t) = q dist(x,Ωi) ε + O(ε 2 ), u ε j (x, t) = 1 -q dist(x,Ωi) ε + O(ε 2 ), u ε k (x, t) = O(ε 2
), for all k ∈ {1, 2, . . . , N } \ {i, j}.

In conclusion, the new proposed phase field model is formally proved to have an order of convergence in O(ε 2 ).

Numerical modeling and storage

This section is devoted to the numerical scheme used to compute the solution of the modified phase field model :

∂ t u i (x, t) = ∆u i (x, t) - 1 ε 2 W (u i ) + λ(x, t) 2W (u i ) + µ i (t) 2W (u i ), for all i ∈ {1, 2, . . . , N }.
We consider its solution for times t ∈ [0, T ], in a computation box Q with periodic boundary condition and with the associated initial condition u(x, 0) = u 0 . We also assume that u 0 satisfies the partition constraint N i=1 u 0 i = 1 and we introduce the mass of u i defined by m i (t) = ´Q u i (x, t)dx. Recall that λ and µ i are free Lagrange multipliers that impose the constraints:

N i=1 u i (x, t) = 1,
and

m i (t) = m i (0) = m 0 i .
Various numerical methods have been used to solve phase field model, for instance, the finite difference method [START_REF] Chen | Convergence of numerical solutions to the allen-cahn equation[END_REF][START_REF] Li | An unconditionally stable hybrid numerical method for solving the allen-cahn equation[END_REF][START_REF] Bates | Numerical analysis for a nonlocal allen-cahn equation[END_REF] and the finite element method [START_REF] Feng | A posteriori error estimates and an adaptive finite element method for the allen-cahn equation and the mean curvature flow[END_REF][START_REF] Feng | Analysis of a fully discrete finite element method for the phase field model and approximation of its sharp interface limits[END_REF][START_REF] Zhang | Numerical studies of discrete approximations to the allen-cahn equation in the sharp interface limit[END_REF]. In this work, having in mind the periodic boundary condition, we focus on a semi-implicit splitting Fourier spectral method [START_REF] Chen | Applications of semi-implicit Fourier-spectral method to phase field equations[END_REF][START_REF] Bretin | A new phase field model for inhomogeneous minimal partitions, and applications to droplets dynamics[END_REF] in order to get high accuracy approximation in space. More precisely, our numerical approach consists in using a splitting scheme between the L 2 -gradient flow of the Cahn Hilliard system and the projection of the current solution to the partition and volume constraints. Notice that this exact discrete projection is slightly different from using a direct discretization of the explicit Lagrange multipliers ( 6) and [START_REF] Baldo | Minimal interface criterion for phase transitions in mixtures of cahn-hilliard fluids[END_REF]. Indeed, the naive discretization of the Lagrange multipliers conduces to cumulative errors, because of the truncation errors induced by our compact and sparse storage strategy.

More precisely, we introduce an approximation sequence u n of u at time nδ t defined recursively as follows :

Step 1: L 2 -gradient flow of the Cahn Hilliard energy without constraint:

Let u n+1/2 be an approximation of v(δ t ) where v = (v 1 , v 2 , . . . , v N ) is the solution of

∂ t v(x, t) = ∆v(x, t) -1 ε 2 W (v(x, t)), ∀(x, t) ∈ Q × [0, δ t ], v(x, 0) = u n (x), ∀x ∈ Q with periodic boundary conditions.
Here (W (v(x, t))

) i = W (v i (x, t)).
Step 2: Projection on the partition and volume constraints: Define u n+1 by

u n+1 i = u n+1/2 i + λ n+1 2W (u n+1/2 i ) + µ n+1 i 2W (u n+1/2 i
), for all i ∈ {1, 2, . . . , N }.

Here λ n+1 and µ n+1 i are defined by satisfying the discrete constraints

N i=1 u n+1 i = 1 and ˆQ u n+1 i = m 0 i .

3.1

Step one with an semi-implicit Fourier spectral scheme

We use a semi-implicit scheme to compute the solution u n+1/2 . More precisely, we consider the equation

Id -δ t ∆ -α/ε 2 Id u n+1/2 = u n - δ t ε 2 (W (u n ) -αu n ) ,
where α is a positive stabilization parameter. Notice that the Cahn Hilliard energy decreases unconditionally [START_REF] Eyre | Computational and mathematical models of microstructural evolution[END_REF][START_REF] Shen | Second-order convex splitting schemes for gradient flows with ehrlich-schwoebel type energy: Application to thin film epitaxy[END_REF] as soon as the explicit part s → W (s) -αs is the derivative of a concave function. This is true for a sufficiently large coefficient α and for the case of W (s) = 1 2 s 2 (1 -s) 2 , and the limit is reached at α = 2. Notice also that without the stabilization parameter (i.e. α = 0), the semi-implicit scheme is also stable under the classical condition

δ t ≤ C ε 2 , where C = s∈[0,1] |W (s)|.
Finally, recall that the equation is computed in Q with periodic boundary conditions. Consequently, the inverse of the operator Id -δ t ∆ -α/ε 2 Id is easily computed in Fourier space [START_REF] Chen | Applications of semi-implicit Fourier-spectral method to phase field equations[END_REF] using the Fast Fourier Transform.

Step two using an exact discrete constraints

We now explain how we can find a solution λ n+1 and µ n+1 i so that

N i=1 u n+1 i = 1 and ˆQ u n+1 i = m 0 i .
and where

u n+1 i = u n+1/2 i + λ n+1 2W (u n+1/2 i ) + µ n+1 i 2W (u n+1/2 i
), for all i ∈ {1, 2, . . . , N }.

Let us first introduce λ

n+1 i = ´Q 2W (u n+1/2 i
)λ n+1 dx, and remark that integrating the last equation over Q leads to

µ n+1 i = m 0 i -´Q u n+1/2 i dx -λ n+1 i ´Q 2W (u n+1/2 i
)dx and

λ n+1 = 1 -k u n+1/2 k -k µ n+1 k 2W (u n+1/2 k ) k 2W (u n+1/2 k ) .
Moreover, the coefficients λ

n+1 i satisfy λ n+1 i = ˆQ λ n+1 2W (u n+1/2 i )dx, = ˆQ 2W (u n+1/2 i ) 1 -k u n+1/2 k k 2W (u n+1/2 k ) dx - ˆQ k µ n+1 k 2W (u n+1/2 k ) 2W (u n+1/2 i ) k 2W (u n+1/2 k ) dx, = k    ˆQ 2W (u n+1/2 k ) 2W (u n+1/2 i )/ ´Q 2W (u k ) j 2W (u n+1/2 j ) dx    λ n+1 k + ˆQ       1 - k u n+1/2 k - k m 0 k -´Q u n+1/2 k dx 2W (u n+1/2 k ) ´Q 2W (u n+1/2 k )dx    2W (u n+1/2 i ) N k=1 2W (u n+1/2 k )    dx.
In particular, it implies that λ = (λ 1 , λ 2 , . . . , λ N ) is solution of the linear system

(Id -A)λ = b,
where

A i,k =    ˆQ 2W (u n+1/2 i ) 2W (u n+1/2 k )/ ´Q 2W (u n+1/2 k )dx N j=1 2W (u n+1/2 j ) dx    and b i = ˆQ       1 - k u n+1/2 k - k m 0 k -´Q u n+1/2 k dx 2W (u n+1/2 k ) ´Q 2W (u n+1/2 k )dx    2W (u n+1/2 i ) N k=1 2W (u n+1/2 k )    dx.
Notice that i A k,i = 1 and the matrix (Id -A) is not invertible. Otherwise, it is not difficult to see that b i = 0 as m 0 i = |Q|. This means that the linear system (Id -A)λ = b admits an infinity of solutions, and, in practice, we choose the solution which satisfy the following constraint i λ i = 0.

Spatial discretization and discrete Fourier transform

We recall that the P Fourier approximation of a 3D function u in a box

Q = [0, L 1 ] × [0, L 2 ] × [0, L 3 ] is given by u P (x) = P/2 p1,p2,p3=-P/2+1 c p e 2iπξp•x
where p = (p 1 , p 2 , p 3 ) and ξ p = (p 1 /L 1 , p 2 /L 2 , p 3 /L 3 ). Here c p represents the P 3 first discrete Fourier coefficients of u. Moreover, the inverse discrete Fourier transform of c p leads to u P p = IFFT[c p ] where u P p is the value of u at the points x p = (p 1 h 1 , p 2 h 2 , p 3 h 3 ) where h α = L α /P for α ∈ {1, 2, 3}. Conversely, c p can be computed by applying the discrete Fourier transform to u P p :

c p = FFT[u P p ].
We present some useful discretization operators using Fourier discretization :

Homogeneous differential operator: the operator Id -δ t ∆ -α/ε 2 Id -1 can be computed in Fourier space by using

Id -δ t ∆ -α/ε 2 Id -1 u P (x) = P/2 p1,p2,p3=-P/2+1 c p 1 + δ t (4π 2 ξ p 2 + α/ε 2 ) e 2iπξp•x .
In practice, it follows that

Id -δ t ∆ -α/ε 2 Id -1 u P (x p ) = IFFT FFT[U N n ] 1 + δ t (4π 2 ξ p 2 + α/ε 2 )
.

Computation of integrals: remark also that

ˆQ u P dx = L 1 L 2 L 3 c 0 = h 1 h 2 h 3 p u P p .
More generally, for all continuous function F : R → R, we use the following approximation ˆQ F (u

P )dx h 1 h 2 h 3 p F (u P n ).
Notice that due to periodicity, this simple quadrature formula is expected to have an order of convergence higher than any integer power of the grid step.

Optimization of the size of the box Q: using a given discrete function u P p , we also propose to optimize the size of the box Q which minimizes the Cahn Hilliard energy

J ε (u P p ) = ˆQ ε |∇u P n | 2 2 + 1 ε W (u P p )dx under the volume constraint L 1 L 2 L 3 = 1.
As ´Q W (u P p )dx depends only on the product L 1 L 2 L 3 , this optimization can done by minimizing the functional

J(L 1 , L 2 , L 3 ) = p |∂ 1 u P p | 2 + |∂ 2 u P p | 2 + |∂ 3 u P p | 2 ,
where

∂ l u P p = IFFT[2iπ p l L l FFT[u P p ]], for l ∈ {1, 2, 3}.
Then, we have

J(L 1 , L 2 , L 3 ) = A 1 L 2 1 + A 2 L 2 2 + A 3 L 2 3 ,
where

A l = p | IFFT[2iπp l FFT[u P p ]]| 2
does not depend on the length L 1 ,L 2 and L 3 . Finally, it follows that optimal values can be explicitly computed by the fomula:

L 1 = √ A 1 (A 1 A 2 A 3 ) 1/6 , L 2 = √ A 2 (A 1 A 2 A 3 ) 1/6 and L 3 = √ A 3 (A 1 A 2 A 3 ) 1/6 .
Notice that in the new box, u P p corresponds now to the value of U P at the point x p = (p 1 L 1 /P, p 2 L 2 /P, p 3 L 3 /P ).

Zero approximation and efficient storage of phase field functions

Let τ ≥ 0 be a positive tolerance, and let us introduce the approximation function ξ τ as:

ξ τ (v) = 0 if v < τ, v otherwise. ( 7
)
We propose to use this approximation function to cut-off u ε i to 0 when his value is considered insignificant, that is when it is below τ . More precisely, we propose to apply this approximation function on each of the two solving steps described in section 3, that is:

Step 1: L 2 -gradient flow of the Cahn Hilliard energy without constraint:

Let ũn+1/2 be an approximation of v(δ t ) where v = (v 1 , v 2 , . . . , v N ) is the solution of

∂ t v(x, t) = ∆v(x, t) -1 ε 2 W (v(x, t)), ∀(x, t) ∈ Q × [0, δ t ], v(x, 0) = u n (x), ∀x ∈ Q with periodic boundary conditions.
Here (W (v(x, t))

) i = W (v i (x, t)).
Let u n+1/2 be the zero-approximation of ũn+1/2 , such that:

u n+1/2 i (x) = ξ τ ũn+1/2 i (x) , for all x ∈ Q.
Step 2: Projection on the partition and volume constraints: Define ũn+1 by

ũn+1 i = u n+1/2 i + λ n+1 2W (u n+1/2 i ) + µ n+1 i 2W (u n+1/2 i
), for all i ∈ {1, 2, . . . , N }.

Here λ n+1 and µ n+1 i are defined by satisfying the discrete constraints

N i=1 ũn+1 i = 1 and ˆQ ũn+1 i = m 0 i .
Let u n+1 be the zero-approximation of ũn+1 , such that:

u n+1 i (x) = ξ τ ũn+1 i (x) , for all x ∈ Q.
Now that we have reduced the amount of significant values, we focus on an adapted storage. Since, for a given discrete space point x, only a subset of the phase field functions is significant (i.e. non-zero), we have then to store an index list I(x) of significant phase-field functions:

I(x) := {i ∈ 1, . . . , N ; u i (x) = 0},
and the associated values U i (x) for i ∈ I(x) such that

u i (x) = U i (x) if i ∈ I(x), 0 otherwise.
Thus, storing for each discrete point x only the corresponding index list I(x) and associated phase field values {U i (x) ; i ∈ I(x)} leads to a significant reduction of the memory needed to store the N phase-field functions.

The next section describes some implementation details of this strategy and section 5 illustrates their behavior on a test case.

Implementation

While some tests and the validation of the new model have been made using a Matlab implementation, the main code used on large discretization domain with many phases, along with an efficient storage strategy, has been written in modern C++.

We describe below some of the implementation choices we have made for the simulation part and for the data extraction and analysis.

Efficient storage of the phase-field functions. From a computational point of view, storing in a row the index and values for each discrete point leads to low performances. This is due to the unpredictable memory address corresponding to each discrete point and to the need of regular memory re-allocation when the size of I(x) changes.

Observing that, on average, the size of I(x) is very low (less that 3 on a typical problem of threedimensional space tiling, depending on the ratio between the total number of phases and the mesh resolution), we propose to use a fixed size initial storage for each discrete space point, which can be dynamically extended (similarly to a singly-linked list) if more phase field functions must be stored.

Consequently, we can use a contiguous block of memory for the initial value storage for all discretization points with easy predictable memory addresses. This leads to a better use of the processor caching system but also allows the C++ compiler to exploit at best vector processor instructions. Our own experiments have shown that extra storage is only rarely needed and the induced performance penalty at these few places and moments is more than counterbalanced by the global performance gain.

A first implementation of this strategy stores I(x) as a bit field where the i-th bit b i is set if and only if i ∈ I(x). The U i (x) are then stored in the same order as the elements i of I(x) (contiguously if possible, sometimes linked by a pointer if I(x) contains more than 4 elements). The bit field is finaly used to test if a given phase is actually stored and to recover its storage rank. A C++ implementation of this storage strategy is the LabelledMap class of the open-source library DGtal1 .

Another implementation that avoids to have a large bit field (e.g. when N is greather than 64 or 128), is simply to store an unordered list of (i, U i (x)) for each i ∈ I(x). The choice between these two implementations depends on the available memory, the maximum and average number of phases of the studied problem, and the read/write cost.

Discrete Fourier transform. The discrete Fourier transform is done using the FFTW library applied on a full copy of each phase-field, extracted from the sparse storage. The storage is then updated with the result of the first step on the current phase, before considering the next phase, thus leading to one temporary memory allocation for a full phase field function.

Projection step. During the second step, the λi , A i,k and b i coefficients are calculated so that no extra memory allocation is needed, even if some calculations are duplicated and may be cached otherwise.

Convergence criteria. When needed (e.g. for the experiments on Kelvin's problem), deciding when the evolution is considered at convergence is based on the evolution of the cost used by F. Morgan in [START_REF] Morgan | The hexagonal honeycomb conjecture[END_REF] for a partition of R n into countable measurable sets of unit volume.

As noted in [START_REF] Oudet | Approximation of partitions of least perimeter by Gamma-convergence: around Kelvin's conjecture[END_REF], for a finite partition Ω := (Ω i ) i=1,...,N of the computation box Q with periodic boundary condition, the Morgan's cost can be easily computed as

F (Ω) = 1 2 N i=1 P (Ω i ) N 1 n , (8)
where the perimeter P (Ω i ) of each phase is approximated by the Cahn Hilliard's energy (1) scaled by the constant c W := ´1 0 2W (s)ds.

Exported data. When exporting results, in the same way as during the computation, we avoid to store all the phase field function's values. To do so, we choose to export only two useful value fields :

• the label l(x) := arg max i (u i (x)) of the leading phase at the discrete point x,

• a coarse distance function d(x) to the partition's frontier i ∂Ω i , reconstructed from the (u i ) as

d(x) := 2ε arctanh(u a -u b ) where u a = u l(x) (x) and u b = max i =l(x) (u i (x)).
Surface extraction. To obtain a finer surface representation and a better Morgan's cost evaluation of the obtained tilings, we use the very efficient local optimization sofware Evolver (see [START_REF] Brakke | The surface evolver[END_REF]) developed by K. Brakke. To exploit it, we need to extract a mesh of the partition's frontiers and to recover the phase connectivity. This non-trivial process is done through the following steps:

• The subset {d(x) ≤ δ} of the computation grid is extracted so that it is everywhere thick and it contains the partition's frontiers, typically by choosing δ a few times greater than the space discretization step.

• This subset is transformed into a cubical complex (i.e. a set of cubes corresponding to the subset above, as well as all their lower dimensional cells); this complex is then collapsed homotopically into a thin 2-dimensional complex, that is a union of digital surfaces.2

• Digital surfaces are then converted into meshes connected with each other and their connectivity is recovered by tracking the internal surface of each cell of the partition.

These algorithms (cubical complex construction, collapsing process and surface tracking) are parts of the DGtal library.

Accurate evaluation of the Morgan's cost. The extracted mesh and connectivity informations are finally given to the Evolver software, which will then optimize the surface mesh. After many optimization steps, Evolver extracts a more reliable estimation of the Morgan's cost of our partition.

Dimension agnostic. Following the DGtal framework, the same code works for any space dimension and practical tests have been made in 2D, 3D and 4D.

5 Numerical experiments

Comparison of the two phase field models

The first experiment concerns a numerical comparison of the two phase field models: the classical one

∂ t u i (x, t) = ∆u i (x, t) - 1 ε 2 W (u i ) + λ(x, t) + µ i (t) 2W (u i ),
and the slightly modified version

∂ t u i (x, t) = ∆u i (x, t) - 1 ε 2 W (u i ) + λ(x, t) 2W (u i ) + µ i (t) 2W (u i ).
In each case, we consider an approximation of a conserved multiphase mean curvature flow associated to the same initial partition Ω(0), a Q partition with N = 9 phases. Moreover, all numerical experiments have been performed using the following numerical parameters : L = [1, 1, 1], P = 2 8 , ε = 1/P and δ t = 1/ε 2 . We plot on figure (3) the evolution of the phase field function u at different time t. The first row and the second row correspond respectively to the classical and the new phase field model. In particular, we can clearly see that the evolution of the partition seems to be very similar using the two different approaches even if the phase field vector functions u are not identically the same.

To be more precise, we plot in figure (4) the functions i iu i , u 1 , log(|u 1 |) and min(u 1 , 0) at the end of the computation. Notice that the first row and the second row still correspond respectively to the classical and the new phase field models. Similarly to our introductory example, we observe that in the case of the classical phase field model, the phase u 1 has a contribution of magnitude O(ε) on each interface Γ i,j . Concerning the new phase field model, we now observe only a negative part of magnitude 10 -6 (which is about 10 -3 for the classical phase field model) which corresponds to the size of ε 2 .

Numerical analysis of memory storing using the new phase field model

The second experiment is obtained using the new phase field model and where the initial partition Ω(0) is build as a random distribution of the N = 16 phases. Our motivation is to illustrate the quantity of memory required to obtain the numerical solutions u during the computation. As previously, we use the following numerical parameters L = [1, 1, 1], P = 2 8 , ε = 1.5/P and δ t = 1/ε 2 . Here, we additionally apply the zero approximation strategy (with τ = 10 -2 ), presented in the previous subsection, to reduce the amount of required memory. We recall that the idea is to use a sparse representation of data u i such as only non-zero data are stored.

The figure [START_REF] Baldo | Minimal interface criterion for phase transitions in mixtures of cahn-hilliard fluids[END_REF] represents the phase field partition u at iterations 1, 10, 50 and 1000. The first row corresponds to the function i iu i and we plot on the second row the number of non-zero phases u i at each pixels. In particular, starting from the moment where the set Ω i are sufficiently large with respect to the size of ε, we can see that the number of non-zero phases u i is about • 3 on each triple point junctions,

• 2 on diffuse interface Γ i,j ,

• 1 on the interior of set Ω i .

More precisely, we plot on figure [START_REF] Barrett | On the parametric finite element approximation of evolving hypersurfaces in R 3[END_REF] the evolution of the number of non-zero phases per discrete point along the computation. It shows that on average, only 1.5 phases remains on each pixel. A full storage of the phase-field functions would here use 8 MiB of memory. Storing only non-zero phases leads to a memory usage at iterations 1, 10, 50 and 1000 of about, respectively:

The Kelvin's problem

Recall that one of our first motivation in this work concerns the Kelvin's problem. The Kelvin's problem is about filling the space with cells of equal volume such that the total surface area is minimal. In the three-dimensional case, there is no proof of what is this minimal structure. The famous Kelvin structure (second picture on figure [START_REF] Bates | Numerical analysis for a nonlocal allen-cahn equation[END_REF]) is not the best partition for this problem since Weaire and Phelan found a better candidate (see [START_REF] Weaire | A counter-example to kelvin's conjecture on minimal surfaces[END_REF] and first picture on figure [START_REF] Bates | Numerical analysis for a nonlocal allen-cahn equation[END_REF]). Notice that the Weaire and Phelan structure is composed of eight cells that compose a filling of the periodic space [0, 1] 3 . As in [START_REF] Oudet | Approximation of partitions of least perimeter by Gamma-convergence: around Kelvin's conjecture[END_REF], we try to find this structure using our volumetric formulation and to understand if there exists or not better structures.

Our strategy consists then in launching a battery of tests (1655 of them are resumed on figure 8) with different numbers N of phases. To overcome the convergence issues that arise from previous cases, we start from a random initial density on a low-resolution grid (of 64 3 ) to avoid an excessive memory usage and a high ε to have a fast convergence rate. When the system is considered to be at equilibrium, we alternate decreases of ε and increases of the resolution up to 512 3 when needed. To improve our chances to find a better structure, we also apply an optimization of the size of the box Q during the computation as detailed previously.

Finally, the interface between each phase is extracted as a mesh (see section 4) and handed to the Evolver software both for a final optimization step and for evaluating more reliably the Morgan's cost of the partition.

The figure [START_REF] Bellettini | Quasi-optimal error estimates for the mean curvature flow with a forcing term[END_REF] gives a summary of our different results. Each point of this picture represents the result of one numerical experiment and for which, the abscissa and ordinate correspond respectively to the number N of phases used and the Morgan's cost of the final structure. Recall that the Morgan's cost gives a natural criteria of the optimality of a structure which is approximately equal to 2.64417 for Weaire and Phelan structure, and 2.653160 for the Kelvin structure. Unfortunately, in spite of all our numerical experiments, we did not succeed to find a better structure than the one of Weaire and Phelan. However we observed that:

• the Weaire and Phelan structure has been found with N = 8, N = 14 and N = 64 phases (see left of Figure 9),

• the Kelvin structure has often been found, as for instance with N = 128 phases (see right of Figure 9),

• intermediate structures between Kelvin and Weaire and Phelan structures have been found (see subfigures of Figure 10),

• uncommon structures with a small but non-optimal Morgan's cost have been found as for instance the ones in Figure 11.

For interested reader, all our results are available on the following web page http://math.univ-lyon1. fr/ ~denis/honeycomb/.

Gout drop experiments with a great number of phases in 2D and in 3D

Finally, we present in figures 12, 13 and 14 three numerical experiments in 2D and in 3D. The novelty is to grow up the number of phases N along the iterations and then to treat a great number of phases.

First, figure 12 shows a 2D experiment where a new phase is inserted each time the evolution is considered at equilibrium. Inserting a new phase is done by adding random amounts of mass at the exterior boundary of the cluster, proportionately to |u(1 -u)| where u is the phase field associated to the cluster's exterior. We do not expect to obtain the best configuration for each N because of the iterative process (each image is the evolution result of the previous image) and because of the described adding process. However, for N = 1, 2, 3, 4, 5, 7 and 10, we recover proved or conjectured best configurations depicted in [START_REF] Cox | Minimal perimeter for n identical bubbles in two dimensions: calculations and simulations[END_REF].

Figure 13 shows another 2D experiment with growing number of phases from N = 1 to N = 512 (of equal volume 1 512 ), with 512 discrete points along each dimension. However, in this simulation, new phases are added at a regular rhythm, after a fixed number of iterations, thus leading to non-equilibrium states when the adding process occurs.

The last figure 14 is a 3D experiment with growing number of phases from N = 1 to N = 64 with equal volume of 1 64 . Like the first experiment, each phase is added when the evolution is considered at equilibrium.

Conclusion

We have presented in this work a large scale numerical strategy to compute solutions of multiphase field model in the context of a great number of phases. In view of memory storage issues, we have noticed that the solutions to the classical phase field model are not sufficiently well localized to apply the required sparse storage strategy. We have then introduced a new phase field model which presents some better localization properties than the classical one. More precisely, a formal asymptotic analysis is also derived which shows the consistency in O(ε 2 ) of our model (instead of O(ε) for the classical one). Using this new model, we have then explained how to store the data such as the quantity of memory does not depend on the number of phases N . Finally, we applied successfully our strategy and presented some numerical experiments with more than N = 100 phases in 3D and N = 1000 in 2D, which highlight the advantages of our approach.

Figure 1 :

 1 Figure 1: Example of numerical experiment using the classical phase field model with 9 phases and ε = 1/256; At a given time t, up left: the function i iu i where each color matches to a phase ; up right: the first phase u 1 ; down left: the logarithm of the absolute value of u 1 , i.e log(|u 1 |), where a contribution of magnitude O(ε) is localized one every interface Γ i,j ; down Right: the negative part of u 1 , i.e min(u 1 , 0), which appears of magnitude O(ε).

Figure 2 :

 2 Figure 2: Example of numerical experiment using our new proposed variant of the multiphase field model (4) with 9 phases and ε = 1/256; At a given time t : Up left: the function i i u i where each color localizes each phase ; Up right: the first phase u 1 . Down left: the logarithm of the absolute value of u 1 , i.e log(abs(u 1 )). Down Right : the negative part of u 1 i.e min(u 1 , 0) which appears now of magnitude O(ε 2 )

Figure 3 :Figure 4 :

 34 Figure 3: Approximation of a conserved multiphase mean curvature flow: a numerical comparison using the classical phase field model (first row) and our slightly modified variant (second row). The initial partition contains 9 phases. Each picture represents the function i iu i at different times t.

• 8 . 5

 85 MiB, 3.1 MiB, 2.2 MiB and 1.6 MiB for the first implementation based on a bit field, with a first block of 2 values and extra blocks of 2 values, • 8.9 MiB, 2.8 MiB, 1.9 MiB and 1.3 MiB for the second implementation that stores the (i, u i ) couples, with a first block of 2 values and extra blocks of 2 values.

Figure 5 :Figure 6 :

 56 Figure 5: Numerical approximation of a conserved mean curvature flow initialised with a random initial partition. At iterations 1, 10, 50 and 1000: first row: arg max i (u i ) with the {d(x) = 0.05} contour lines (in black) of the reconstructed distance function to each domain Ω i ; second row: the number of non-zero phases u i on each pixel, superposed with the same contour lines as previously. Note that the color scale is different on each image of the second row.

Figure 7 :

 7 Figure 7: Examples of Weaire and Phelan structure (left) and Kelvin structure (right) obtained with our numerical phase field method with respectively N = 8 and N = 12, starting from a random initial state. The box size optimizer leads to a non-square calculation domain for the Kelvin structure with a stack direction that is not aligned with the box orientation.

Figure 8 :

 8 Figure 8: Summary of 1655 numerical experiments run to find optimal structures. Each cross corresponds to one numerical optimization : the number of phases N in abscissa, and the Morgan's cost of the final structure in ordinate. Each cross is surrounded by a blue dot whose opacity is linked to the number of experiments with a similar Morgan's cost.

Figure 9 :

 9 Figure 9: Examples of Weaire and Phelan structure (left) and Kelvin structure (right) obtained with our numerical phase field method and with a great number of phases, N = 64 and N = 128, starting from a random initial state. The box size optimization leads here to a near square domain, up to the third digit.

Figure 10 :

 10 Figure 10: Examples of optimal structures with a Morgan cost better than Kelvin one obtained with N = 60 (left) and N = 44 (right) phases.

Figure 11 :

 11 Figure 11: Examples of uncommon structures obtained with N = 38 (left) and N = 24 (right) phases.

Figure 12 :

 12 Figure 12: Example of a numerical experiment in 2D where a new phase of fixed volume is added each time the evolution is at equilibrium. Each picture represents the equilibrium state before adding the next phase.

Figure 13 :

 13 Figure 13: Example of a numerical experiment in 2D (512 discrete points along each dimension) where the number of phases is growing up along the iterations at a regular rhythm, up to N = 512 phases of volume 1 512 . Each picture represents the partition at a given time t.

Figure 14 :

 14 Figure 14: Example of a numerical experiment in 3D (128 discrete points along each dimension) where the number of phases is growing up along the iterations, up to N = 64 phases of volume 164 . New phases are added each time the evolution is considered to be at equilibrium. Each picture represents the partition at a given time t.

http://www.dgtal.org/

The homotopic collapse is made through elementary collapse operations, which were introduced by Whitehead[START_REF] Whitehead | Simplicial spaces, nuclei and m-groups[END_REF]. It consists in finding couples of cells (c, d), such that c has dimension one less than d, c lies in the boundary of d but not in the boundary of any other cell, d lies not in the boundary of any cell, and in removing them in sequences.

Acknowledgments The authors acknowledge the support of the French National Research Agency (ANR) under grant ANR-11-BS02-0009 (project digitalSnow) and grant ANR-15-CE40-0006 (project Comedic). This work was also supported by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).


Hoai-Nam Nguyen

Paolino Tona

Guillaume Sabiron

Dominant Wave Frequency and Amplitude Estimation for Adaptive Control of Wave Energy Converters

Adaptive control is of great interest for wave energy converters (WEC) due to the inherent time-varying nature of sea conditions. Robust and accurate estimation algorithms are required to improve the knowledge of the current sea state on a wave-to-wave basis in order to ensure power harvesting as close as possible to optimal behavior. In this paper, we present a simple but innovative approach for estimating the wave force dominant frequency and wave force dominant amplitude based on the unscented Kalman filter. The main benefits of the proposed method are illustrated by a numerical example concerning a small-scale wave energy converter prototype, with data coming from experimental recordings obtained in a test basin with irregular waves.

I. INTRODUCTION

Precise knowledge of current sea-state characteristics is crucial in advanced control methods for wave energy converters (WECs). More specifically, having real-time information at disposal about dominant frequency and dominant amplitude of the wave excitation force acting on the device, enables the controller to calculate the control action in an optimal manner.

Les us recall that the ultimate objective in the hydrodynamic control of a WEC is the maximisation, for a broad range of sea states, of the energy captured from the waves, via the force applied by the power-take off (PTO) system (or secondary converter) to the primary converter (f P T O in Fig. 1 which provides a simplified representation of a heaving-buoy wave energy converter). The reader can refer to [START_REF] Korde | Hydrodynamic Control of Wave Energy Devices[END_REF] for an extended definition of primary and secondary converters. Latching control, model predictive control (MPC) and PID control are examples of implementable WEC control strategies [START_REF] Ringwood | Control Optimisation and Parametric Design[END_REF].

For WECs with a flexible PTO, capable of both harvesting and drawing power from the grid (respectively in generator and motor modes), a commonly adopted solution, due to its simplicity, is the so called PI (or BK) controller described in [START_REF] Sanchez | Control performance assessment and design of optimal control to harvest ocean energy[END_REF]. This control law consists in computing the PTO force via a combined proportional feedback on the position and velocity of the primary converter (the float on Fig. 1) 1 . The damping controller (or P controller), where the PTO force is specified to be proportional to and oppositely directed to the velocity of the primary converter, is even more widely adopted, as it does not require the use of reactive power.

To maximise WEC performance, the P and I gains of these controllers cannot be kept constant. They have to be changed following the changes of sea state. Adaptive P or PI control has been already discussed in the literature [START_REF] Whittaker | The development of oystera shallow water surging wave energy converter[END_REF], [START_REF] Hansen | Modelling and control of the wavestar prototype[END_REF], [START_REF] Hals | A comparison of selected strategies for adaptive control of wave energy converters[END_REF], at different level of details, in the form of a gain-scheduling approach. The main idea is to compute optimal gains for a representative set of sea states. The gains are calculated offline, analytically or numerically, using a gridding approach: for each sea state, the ones leading to the best average power are picked. Then the appropriate control action f P T O is found from a look-up table whose input is the current sea state, which can be identified, for instance, in terms of dominant wave frequency and amplitude, or other characteristics of its spectrum. In most literature, only average online estimations of sea states have been proposed, with time windows of several minutes (10 minutes in [START_REF] Hansen | Modelling and control of the wavestar prototype[END_REF], 20-30 minutes in [START_REF] Hals | A comparison of selected strategies for adaptive control of wave energy converters[END_REF]). Such intermittent adaptive control laws are clearly suboptimal in terms of energy recovery, since the control gains are not continuously updated whereas the sea state is continuously time-varying.

By contrast, an efficient method to estimate online the dominant frequency and dominant amplitude of wave excitation force, would allow to adapt the control gains on a wave-towave basis. This would enable the design of a continuously adaptive (PI) control system such as the one depicted in Fig. 2, provided that the wave force excitation can be estimated online.

An estimation approach based on an Extended Kalman Filter (EKF), applied to a nonstationary, harmonic approximation of the wave excitation force, is proposed in [START_REF] Fusco | A simple and effective real-time controller for wave energy converters[END_REF] as part of the "simple-but-effective" WEC controller. The proposed approach consists in a velocity tracking controller with adaptive parameters computed using online estimates of the dominant frequency and dominant amplitude of the wave Fig. 2. Adaptive WEC control system with wave force frequency and amplitude estimation excitation force. The rationale behind the choice of such an estimation approach is that the wave force can be seen as a sum of a number of sinusoids (assuming a panchromatic sea state, [START_REF] Falnes | Ocean waves and oscillating systems: linear interactions including wave-energy extraction[END_REF]), but its dependency on frequency and amplitude is a nonlinear function, preventing the use of standard linear Kalman filter (KF). Moreover, it is quite natural to resort to robust methods to obtain accurate estimates in the presence of noise and/or uncertainties in wave force measurements or estimates.

Nonetheless, it must be recalled that the EKF provides a solution to a nonlinear estimation problem via local linearisations of the underlying model. This entails several drawbacks for this applications. In particular, when the variation of the wave force is large and/or the sampling time intervals are not sufficiently small, linearisation can produce highly unstable filters, potentially leading to divergence phenomena [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF], [START_REF] Julier | New extension of the Kalman filter to nonlinear systems[END_REF]. In practice, this jeopardises the ability to accurately track dominant wave frequency and amplitude, when the sea state changes.

To address these limitations, the Unscented Kalman Filter (UKF) is proposed in [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF]. This is a powerful nonlinear estimation technique and has been shown to be a superior alternative to the EKF in several applications including parameter estimation for time series modeling [START_REF] Wan | The unscented kalman filter for nonlinear estimation[END_REF], neural network training [START_REF] Van Der Merwe | Efficient derivative-free Kalman filters for online learning[END_REF], and state estimation for road vehicle navigation [START_REF] Julier | Navigation and parameter estimation of high speed road vehicles[END_REF]. An interesting feature of UKF is that no linearisation of the nonlinear model is required.

Although the UKF has been used in a wide range of estimation problems, to the best of the authors' knowledge there has been no attempt to use it to estimate the dominant frequency and amplitude of the wave excitation force. In this paper, we aim to fill this gap by exploring the potential benefits of the UKF in wave frequency and amplitude estimation. Real data recordings are used to analyse the validity and performance of the proposed UKF algorithm. This paper is organised as follows. Section 2 covers the problem formulation. Section 3 is dedicated to the unscented Kalman filter design. In Section 4, a state space model, that is used in the UKF, is presented. Results obtained with experimental data are evaluated in Section 5 before drawing the conclusion.

II. PROBLEM FORMULATION

In this paper we are interested in estimating in real-time the dominant frequency and dominant amplitude for a given wave excitation force acting on the WEC. The dominant frequency of the wave force is defined as the peak of its spectrum, obtained by using the Fourier transform of the corresponding time signal. For example, in Fig. 3, the dominant wave force frequency is about 1.75 rad/s. The wave force dominant amplitude is defined as the best magnitude, that the wave force can be approximated to by a sinusoid with the dominant frequency and amplitude.

III. UNSCENTED KALMAN FILTERING

The main idea behind the unscented Kalman filter is that it is easier to approximate a Gaussian distribution than it is to approximate an arbitrary nonlinear function.

Instead of linearising the nonlinear functions, the UKF uses a deterministic sampling approach to capture the mean and covariance estimates with a minimal set of sample points.

A. Unscented transformation

The idea of Unscented Transformation (UT) is:

• To choose deterministically 2n + 1 sample points (also called sigma points) with mean x and covariance matrix P xx . Here n is the dimension of x.

• To find the mean ȳ and covariance matrix P yy by propagating the sample points through a nonlinear function. Define S as Cholesky factorisation of P xx , i.e. P xx = S T S. In UT, the sample points can be obtained as

   χ 0 = x χ i = x + √ n + λS i , i = 1, 2 . . . , n χ i+n = x - √ n + λS i , i = 1, 2 . . . , n (1)
where S i is the i-th column of the matrix S, and λ = α 2 (n+ κ)n. The scaling parameter α determines the spread of the sample points and κ is a secondary scaling parameter.

Once the sample points are computed, we propagate them through the nonlinear function:

γ i = f (χ i ), i = 0, 1, . . . , 2n (2)
The next step is to calculate the mean ȳ and the covariance matrix P yy of the propagated points as

       ȳ = 2n i=0 W m i γ i , P yy = 2n i=0 W c i (γ i -ȳ)(γ i -ȳ) T (3)
where the weights W m i and

W c i are          W m 0 = λ n+λ , W m i = 1 2(n+λ) , W c 0 = λ n+λ + 1 -α 2 + β, W c i = 1 2(n+λ) (4)
Note that β is a parameter used to incorporate any prior knowledge about the distribution of x.

B. Unscented Kalman filter

The UKF is an algorithm which uses an underlying process model to make an estimate of the current state of a system and then corrects the estimate using any available sensor measurements.

Consider the following system

x(k) = f (x(k -1)) + ν(k -1), y(k) = g(x(k)) + µ(k) (5)
where x(k) ∈ R n is the state vector, y(k) ∈ R m is the measurement vector, ν(k) and µ(k) are respectively, the process noise and the measurement noise. It is assumed that ν(k) and µ(k) are Gaussian noises with zero mean and with covariance matrices Q and R, respectively. For the given state estimate x(k -1), and the given covariance matrix estimate P (k -1), there are three steps in the UKF.

1. Sigma point calculation. Using the state estimate x(k -1) and the covariance matrix estimate P (k-1) at time k-1, the sigma points are calculated, at time k, as

   χ 0 = x(k -1) χ i = x(k -1) + √ n + λS i , i = 1, 2 . . . , n χ i+n = x(k -1) - √ n + λS i , i = 1, 2 . . . , n (6)
where S is the Cholesky factorisation of P (k -1).

2. State prediction. The predicted state mean vector x(k|k -1), and the predicted covariance matrix P (k|k -1) are computed as

           x(k|k -1) = 2n i=0 W m i γ i , P (k|k -1) = 2n i=0 W c i (γ i -x(k|k -1))(γ i -x(k|k -1)) T +Q (7) where γ i = f (χ i ), i = 0, 1, . . . , 2n (8)
3. State correction.

The sigma points related to the predicted state mean vector and the covariance matrix are calculated as    ψ 0 = x(k|k -1)

ψ i = x(k|k -1) + √ n + λT i , i = 1, 2 . . . , n ψ i+n = x(k|k -1) - √ n + λT i , i = 1, 2 . . . , n (9)
where T is the Cholesky factorisation of P (k|k -1). The sigma points are then propagated through the output equation:

ζ i = g(ψ i ) (10)
The following step is to calculate the output mean vector, the output covariance matrix and the cross covariance matrix as

               ŷ(k|k -1) = 2n i=0 W m i ζ i , P yy = 2n i=0 W c i (ζ i -ŷ(k|k -1))(ζ i -ŷ(k|k -1)) T + R P xy = 2n i=0 W c i (γ i -x(k|k -1))(ζ i -ŷ(k|k -1)) T (11)
The Kalman filter gain, the state estimate and the covariance matrix at time k are given by

   K(k) = P xy P -1 yy x(k) = x(k|k -1) + K(k)(y(k) -ŷ(k|k -1)), P (k) = P (k|k -1) -K(k)P yy K(k) T (12)

C. UKF parameters

The UKF has five parameters: Q, R, α, β and κ. The matrices Q and R represent the process noise covariance and the measurement noise covariance. R is determined empirically and accounts for the uncertainty in the output measurements. Choosing the process noise covariance matrix is not straightforward. In general, there is no systematic way to calculate Q.

The scale parameter α determines the spread of sigma points and is set to be 10 -3 in this work. β is used to incorporate prior knowledge of the distribution of x; for a Gaussian distribution, β = 2 is optimal. In this work, κ = 0 is used.

IV. STATE SPACE MODEL

It is assumed that the wave excitation force can be modeled as a time varying sinusoidal signal, i.e.,

f ex (t) = A(t) sin(ω(t)t + φ(t)) + µ(t) (13)
where f ex (t) is the wave force at time t, A(t) is the timevarying dominant magnitude, ω(t) is the time-varying dominant frequency, φ(t) is the time-varying phase angle, and µ(t) is a zero mean random noise. In this model the other harmonics are considered as a random noise.

A state space model is needed for the UKF. Assume that f ex (t) is measured and/or estimated at the times kT s , where k = 0, 1, 2, . . ., and T s is the sampling time. Define the following three states :

   x 1 (k) = A sin(kT s ω + φ), x 2 (k) = A cos(kT s ω + φ), x 3 (k) = ω (14)
Here we reformulate the state equations in a recursive way to cope with UKF requirements. In order to do that, A, ω and ϕ are considered constant at first. Using the following trigonometric relations sin(α + β) = sin(α) cos(β) + cos(α) sin(β), cos(α + β) = cos(α) cos(β) -sin(α) sin(β) [START_REF] Ling | Real time estimation and prediction of wave excitation forces on a heaving body[END_REF] one has

x 1 (k) = A sin ((k -1)T s ω + T s ω + φ) = A sin((k -1)T s ω + φ) cos(T s ω)+ +A cos((k -1)T s ω + φ) sin(T s ω) = cos(T s x 3 (k -1))x 1 (k -1)+ + sin(T s x 3 (k -1))x 2 (k -1) (16)
and

x 2 (k) = A cos ((k -1)T s ω + T s ω + φ) = A cos((k -1)T s ω + φ) cos(T s ω)- -A sin((k -1)T s ω + φ) sin(T s ω) = -sin(T s x 3 (k -1))x 1 (k -1)+ + cos(T s x 3 (k -1))x 2 (k -1) (17) Denote x(k) = [x 1 (k) x 2 (k) x 3 (k)] T
Denote also 16), ( 17), and assuming ω(k) changes slowly with the time, the following state space model is obtained

A x =   cos(T s x 3 (k -1)) sin(T s x 3 (k -1)) 0 -sin(T s x 3 (k -1)) cos(T s x 3 (k -1)) 0 0 0 1   (18) Combining (
x(k) = A x x(k + 1) + ν(k) ( 19
)
where

ν(k) = [ν 1 (k) ν 2 (k) ν 3 (k)] T
Process noises ν(k -1) are included in the state equation ( 19) to take into account the time varying nature of A(k), ω(k) and ϕ(k), although ( 16), [START_REF] Nguyen | Experimental validation of a nonlinear MPC strategy for a wave energy converter prototype[END_REF] are exact. This is explained by the fact that ( 16), [START_REF] Nguyen | Experimental validation of a nonlinear MPC strategy for a wave energy converter prototype[END_REF] do not take into consideration the timevarying nature of A(t), ω(t), and φ(t). Hence ν 1 (k-1), ν 2 (k-1) are used to model these unknown uncertainties. Clearly, ν 1 (k -1), ν 2 (k -1) and ν 3 (k -1) are correlated.

For the output equation, denoting y(k) = f ex (k) and

C = [1 0 0], one obtains y(k) = Cx(k) + µ(k) (20)
where µ(k) is used to model the measurement error in the output equation and to take into account the fact that f ex (k) contains other harmonics. Combining (19), (20), one obtains,

x(k) = A x x(k -1) + ν(k -1) y(k) = Cx(k) + µ(k) (21)
The UKF presented in section III.B is then applied to the nonlinear system defined in (21) to provide estimates of x 1 (k),

x 2 (k) and x 3 (k). From the estimates, one can then derive dominant wave force frequency estimation and dominant wave force amplitude estimation:

ω(k) = x3 (k) (22)
Note that the magnitude A(k) is not a state of the system (21). However, once x(k) is available, it can be easily calculated as

Â(k) = x1 (k) 2 + x2 (k) 2 (23)
V. EXPERIMENTAL STUDY

A. Experimental setup

The experimental work was carried out in a wave basin of Aalborg University (AAU), Denmark, see Fig.

B. Wave excitation force estimation

Note that the wave excitation force f ex (t) is considered to be a known input for the UKF. However, this quantity cannot be directly measured during normal operation. Hence f ex (t) has to be estimated via the measurements of other quantities.

In [START_REF] Kracht | Performance improvement of a point absorber wave energy converter by application of an observer-based control: Results from wave tank testing[END_REF], using a combination of a bank of independent harmonic oscillators and a Luenberger observer, an estimation algorithm for wave excitation force is proposed. The strategy is tested on a real WEC system. However the reported experimental results show a relatively high phase lag in the estimated signal compared to the measured signal.

In [START_REF] Ling | Real time estimation and prediction of wave excitation forces on a heaving body[END_REF], by considering the wave excitation force as a timevarying sinusoid, an EKF approach is presented 2 . However, no experimental results are reported. In addition, the approach can clearly be effective only for very narrow-banded wave forces.

By combining several pressure sensors at discrete points on the buoy surface with the buoy heave position, and with an extended Kalman filter, another approach is proposed in [START_REF] Abdelkhalik | Estimation of excitation forces for wave energy converters control using pressure measurements[END_REF]. However, the computational complexity may be high. In addition, several extra pressure sensors are required.

In [START_REF] Nguyen | Experimental validation of a nonlinear MPC strategy for a wave energy converter prototype[END_REF], the authors have reported the successful implementation of a new estimation method developed to overcome the aforementioned drawbacks.

Its most important features, in the context studied here, are: 1) only standard WEC measurements (position, velocity, PTO force) are used by the algorithm. 2) the experimental results show that estimated wave torque values do not have any significant lag compared to "true" values, which is not the case for the approach in [START_REF] Kracht | Performance improvement of a point absorber wave energy converter by application of an observer-based control: Results from wave tank testing[END_REF]. 3) in contrast to [START_REF] Kracht | Performance improvement of a point absorber wave energy converter by application of an observer-based control: Results from wave tank testing[END_REF], no (implicit) unrealistic assumption about the time-invariant nature of the sea state is made, hence any operating condition can be efficiently dealt with. The approach is based on a linear Kalman filter and a random walk model for the variation of the wave excitation force. Since wave force estimation is not the main subject of the paper, it will not be discussed further. For more details, the interested reader is referred to [START_REF] Nguyen | Wave excitation force estimation for wave energy converters of the point absorber type[END_REF].

C. Validation on experimental data

For the present analysis, an irregular wave is considered, see Fig. 5. The wave is unidirectional, two dimensional and long-crested. The data set is measured at a sampling frequency of 1 kHz. To reduce the computational complexity, data downsampling with a factor of 100 is applied. The spectrum of the data set is presented in Fig. 6.

In Fig. 6, it can be observed that the dominant frequency is around 5 rad/s. 2 Note that the EKF approach mentioned here is applied to the estimation of wave excitation force from WEC measurement, and bears no relation to the UKF approach for the estimation of wave force dominant frequency and amplitude described in this paper, or the EKF mentionned in [START_REF] Fusco | A simple and effective real-time controller for wave energy converters[END_REF] This fact is confirmed by the dominant frequency estimate shown in Fig. 7. Finally, Fig. 8 shows the dominant magnitude estimate. An interesting feature of the UKF is that it is able to estimate the dominant frequency and dominant magnitude even in presence of significant variations of the wave force. This feature can hardly be guaranteed with the EKF.

VI. CONCLUSION

A simple solution based on the unscented Kalman filter to the problem of estimating the wave force dominant frequency and dominant amplitude has been proposed in this paper.

Various experimental results have been presented to assess the performance of the UKF. It has been shown that the UKF estimates with high accuracy and no significant lag the dominant frequency and dominant amplitude in spite of significant variations of the wave excitation force. This innovative algorithm thus paves the way to fully functional adaptive control for wave energy converters of the heavingbuoy type (and possibly of other types). It is interesting to note that the proposed algorithm is robust with respect to WEC architecture since the dynamic behaviour of the system does not interfere with the estimation approach and can thus be applied to any kind of WEC system. Future work will focus on estimating residual bias on both wave force dominant amplitude and frequency.

Fig. 1 .

 1 Fig. 1. Schematic diagram of a generic (heaving-buoy) wave energy converter

Fig. 3 .

 3 Fig. 3. Typical wave force spectrum for a fully developed sea (experimental data).

  4. It consists of four parts: the wave generation (a), the wave propagation zone (b), the point absorber laboratory model (c), and the wave termination (d). Clearly, the wave generator and the point absorber are the main interest in the context of this work.

Fig. 4 .

 4 Fig. 4. Experimental setup at AAU. The waves are generated by the wave generation system (a) and the wave elevation is measured by a series of wave gauges (b). The waves then reach the WEC laboratory model, where the wave energy absorption takes place (c). Finally, they reach the artificial beach where the termination takes place by a scattering process (d).

Fig. 5 .

 5 Fig. 5. Irregular wave.

Fig. 6 .

 6 Fig. 6. Spectrum of sample data set.

Fig. 7 .

 7 Fig. 7. Dominant frequency estimation.

Fig. 8 .

 8 Fig. 8. Dominant magnitude estimation.

This corresponds to a proportional-integral action on the velocity of the primary converter, which is the reason for it being confusingly referred to as PI controller.

ACKNOWLEDGMENTS

The authors wish to thank Morten Kramer, Associate Professor at Aalborg University, for providing data on the smallscale protototype and the experimental set-up, and support during the tests.