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

,


Introduction
This paper is concerned with 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, The evolution of the partition Ω can be obtained as the L 2 -gradient flow (with additional volume constraint) of a multiphase perimeter where Γ i,j = ∂Ω i ∩ ∂Ω j .In particular, the classical physical theory [32,39] 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.
The main goal of this work is to propose an efficient and accurate numerical phase field model [16,37] 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 [9,14,28,35,45].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 The concept of the phase field method consists to approximate the perimeter by the Cahn Hilliard's energy [15], defined for all function u ∈ H 1 (Q, R), by 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 [37] 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 [4]) and is defined by In particular, when Ω is a set with finite perimeter [4,36], 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 Here, dist(x, Ω) is the signed distance function to Ω and q is the profile function associated to the potential W which is defined by 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 [2], obtained as the L 2 -gradient flow of P ε , reads Note that existence, uniqueness, and a comparison principle have been established for this equation (see for example chapters 14 and 15 in [3]).
A smooth mean curvature flow Ω(t) can then be approximated by where u ε is the solution of the Allen Cahn equation with the following initial condition Notice that standard formal asymptotic expansion [14,45] of u ε near the interfaces shows that u ε is quadratically closed to the optimal profile, i.e.
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 The Γ-convergence of P ε to c W P is also established in [44].More general Γ-convergence results in the case of non uniform surface tension σ i,j is obtained in [5,13].Note also that phase field models that handle the case of anisotropic surface tension are introduced and analyzed in [28,29].
The L 2 -gradient flow of P ε translates into the following Allen Cahn system where the Lagrange multiplier λ is defined by the partition constraints N j=1 u j = 1 and is explicitly computed by the formula More precisely, sharp interface limit derived for instance in [28], shows that, at least formally, the solution u ε expands near the interface Γ i,j as with associated normal velocity law V ε i,j satisfying Then we obtained that this phase model is 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 Then, the constrained L 2 -gradient flow of P ε reads as 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 [1,12] 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 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.[31,48,49]), the phase-field function u i is expected of the form 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 (1.3), which does not vanish enough far away from the interfaces.Therefore, a storage strategy based on storing nonnull 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 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 first explanation of the linear convergence is the consequence of the Lagrange multiplier µ i .Indeed, if µ i differs from µ j , the Lagrange multiplier λ becomes active on the interface Γ i,j and modifies the one-order term of the asymptotic expansion of u i and u j .Additionally, without volume conservation, the Lagrange multiplier λ is active only at multiple point junctions and not on the interface Γ i,j which explain why we observe a model of order 2 in that case.

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 [12], we propose to replace λ by λ 2W (u i ), which leads to this new slightly modified phase field model: where λ and µ i are defined accordingly to the following constraints: As previously, straightforward computations show that λ and µ i can be obtained as 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 The only difference between the classical phase-field model and the new version is on the form and treatment of the Lagrange multiplier associated to the constraint N j=1 u j = 1.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 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 (1.2) and the modified approach (1.4).
We follow the method of matched asymptotic expansions proposed in [9,14,35,45] 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 ), for all i ∈ {1, 2, . . ., N }.In particular, analogously to [35], it is not difficult to see that and Inner expansion: in a small neighborhood of Γ i,j , we define the stretched normal distance to the front, We then focus on inner expansions of u ε i (x, t), λ ε (x, t) i.e. expansions close to the front, of the form Additionally, we consider inner expansion of λ ε Moreover, we make the Ansatz that Let us also define a unit normal m i to Γ i,j and the normal velocity V ε i,j to the front as where ∇ refers to spatial derivation only (the same holds for further derivation operators used in the sequel).Following [35,45] 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 [35,45], it is easily seen that Also recall that in a sufficiently small neighborhood of Γ i,j , according to Lemma 14.17 in [30], we have 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 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 [35] for more details) imply in particular that lim z→+∞ and

Classical model and inner expansion
We first consider the phase field model The first order in ε −2 reads for all k ∈ {1, 2, . . ., N } as 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 }, In particular, the cases k = i and k = j read respectively and Multiplying the last equation by q and integrating over R leads to and where we use the matching conditions: lim z→±∞ U 1 j (z, x, t) = 0.
Moreover, the Lagrange multiplier Λ −1 satisfies 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 We then obtain that each U 1 k can be defined as the solution of the system 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 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 (1.4).As previously, the first order in ε −2 reads for all k ∈ {1, 2, . . ., N } as 0, which leads to (with additional boundary conditions) Moreover, the second order in ε −1 shows that for all k ∈ {1, 2, . . ., N }, In particular, using i satisfies the equation (2.1).Multiplying the last equation by q and integrating over R shows that V 0 i,j satisfies (2.2).Consequently, the Lagrange multiplier Λ −1 satisfies (2.3) and is still expected on the form (2.4).The functions U 1 k are then defined as solutions of the system We can now see that U 1 k = 0 for all k ∈ {1, 2, . . ., N } \ {i, j}, and this equality holds even if Finally, the Lagrange multiplier is identified to 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 ), 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 (1.4).
We consider its solution for times t ∈ [0, T ], in a computation box Q with periodic boundary conditions 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: Various numerical methods have been used to solve phase field model, for instance, the finite difference method [8,19,34] and the finite element method [26,27,52].In this work, having in mind the periodic boundary condition, we focus on a semi-implicit splitting Fourier spectral method [13,17] 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 Allen-Cahn 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 (1.6) and (1.5).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 Here (W (v(x, t)) Step 2: Projection on the partition and volume constraints: Define u n+1 by ) , for all i ∈ {1, 2, . . ., N }.
Here λ n+1 and µ n+1 i are defined by satisfying the discrete constraints Remark 3.1.This approximation can be view as a Lie splitting between Cahn-Hilliard energy and the partition and volume constraints, which is known to be a method of order 1 in the best case.
Remark 3.2.In the projections step, the Lagrange multipliers λ n+1 and µ n+1 i are defined so that the discrete constraints are exactly satisfied.In practice, this approach seems to be better than to use a discretization of the analytic expressions of λ and µ i given in 1.6 and 1.5.

Step one with a 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 where α is a positive stabilization parameter.Notice that the Cahn Hilliard energy decreases unconditionally [25,47] 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 [17] 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 , for all i ∈ {1, 2, . . ., N }.

Let us first introduce
λ n+1 dx, and remark that integrating the last equation over In particular, it implies that λ = (λ 1 , λ 2 , . . ., λ N ) is solution of the linear system where 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 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 : 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 In practice, it follows that More generally, for all continuous function F : R → R, we use the following approximation 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.
Here λ n+1 and µ n+1 i are defined by satisfying the discrete constraints Let u n+1 be the zero-approximation of ũn+1 , such that: 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: and the associated values U i (x) for i ∈ I(x) such that 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 three-dimensional 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 ith 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 finally 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 DGtal 5 .
Another implementation that avoids to have a large bit field (e.g. when N is greater 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 Morgan in [38] for a partition of R n into countable measurable sets of unit volume.
As noted in [44], 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 where the perimeter P (Ω i ) of each phase is approximated by the Cahn Hilliard's energy (1.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 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 [11]) 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 2dimensional complex, that is a union of digital surfaces6 .
-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.

Comparison and validation of the two phase field models in the case of two and three bubbles
The first experiment concerns a numerical comparison of the two phase field models: the classical one (1.2) and the slightly modified version (1.4).
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 = 3 and 4 phases.Moreover, all numerical experiments have been performed using the following numerical parameters: L = [1, 1], P = 2 8 , ε = 2/P and δ t = 1/ε 2 .
In Figure 3, we plot the approximation of the evolution of two bubbles obtained with the two different phase field models.We can first observe that the two evolutions are very similar.Moreover, we represent on the left picture of Figure 4, a comparison between the 1/2 level-set of each phase field solution observed at the final time and with the theoretical optimal shape of the two bubbles.This shows a good agreement between the simulation and the expected optimal shape.
We also plot on the right picture of Figure 4, the evolution of the Cahn-Hilliard energy associated to the solutions of the two phase field models.As expected, for each model, the energy decreases in a similar fashion.
Similarly, we make the same experiments with three bubbles plotted on Figures 5 and 6.Finally, all these experiments show that the two phase-field models are equivalent for a small number of phases with quantitative comparisons with the theoretical solutions, which validates our numerical scheme at least in these two cases.

Comparison of the two phase field models with a large number of phases
We now consider an approximation of a conserved multiphase mean curvature flow associated to the same initial partition Ω(0), a Q partition with N = 9 phases.All numerical experiments have been performed using the following numerical parameters: We plot on Figure 7 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 8 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], P = 2 8 , ε = 1.5/P and δ t = 1/ε 2 .Here, we additionally apply the zero approximation strategy presented in the previous subsection, to reduce the amount of required memory, with

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 threedimensional case, there is no proof of what is this minimal structure.The famous Kelvin structure (third and fourth pictures on Fig. 11) is not the best partition for this problem since Weaire and Phelan found a better candidate (see [50] and first and second pictures on Fig. 11).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 [44], 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 Fig. 12) 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.and N = 128 (fourth picture).All experiments have been obtained with our numerical phase field method, 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.
Finally, the interface between each phase is extracted as a mesh (see Sect. 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 12 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.644 17 for Weaire and Phelan structure, and 2.653 160 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 Fig. 11), -the Kelvin structure has often been found, as for instance with N = 128 phases (see Fig. 11),  -intermediate structures between Kelvin and Weaire and Phelan structures have been found (see subfigures of Fig. 13), -uncommon structures with a small but non-optimal Morgan's cost have been found as for instance the ones in Figure 13.
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 14-16 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 14 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 [20].
Figure 15 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 16 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. .New phases are added each time the evolution is considered to be at equilibrium.Each picture represents the partition at a given time t.

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 .
Figure 1.Example of numerical experiment using the classical phase field model with 9 phases and ε = 1/256; At a given time t, from left to right: the function i iu i where each color matches to a phase; the first phase u 1 ; 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 ; the negative part of u 1 , i.e. min(u 1 , 0), which appears of magnitude O(ε).

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

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

Figure 4 .
Figure 4. Multiphase mean curvature flow of two bubbles in 2D.Left panel: comparison of optimal shapes obtained with the two phase field models and the theoretical one.Right panel: evolution of the multi-phase Cahn-Hilliard energy associated to the solutions of the two phase field models.

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

Figure 6 .
Figure 6.Multiphase mean curvature flow of three bubbles in 2D.Left panel: comparison of optimal shapes obtained with the two phase field models and the theoretical one.Right panel: evolution of the multi-phase Cahn-Hilliard energy associated to the solutions of the two phase field models.

Figure 9 .
Figure 9. Numerical approximation of a conserved mean curvature flow initialized 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 10 .
Figure 10.Numerical approximation of a conserved mean curvature flow starting from a random initial partition.Evolution of the number (minimum, maximum and mean values) of non-zero phases per point: from iteration 1 to 50 on the left, from iteration 50 to 1000 on the right.

Figure 11 .
Figure 11.Examples of Weaire and Phelan structure obtained with N = 8 (first picture) and N = 64 (second picture).Examples of Kelvin structure obtained with N = 12 (third picture)and N = 128 (fourth picture).All experiments have been obtained with our numerical phase field method, 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 12 .
Figure 12.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 13 .
Figure 13.Examples of optimal structures with a Morgan cost better than Kelvin one obtained with N = 60 (first picture) and N = 44 (second picture) phases.Examples of uncommon structures obtained with N = 38 (third picture) and N = 24 (fourth picture) phases.

Figure 14 .
Figure 14.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 15 .
Figure 15.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 16 .
Figure 16.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.