A numerical level set method for the Stefan problem
with a crystalline Gibbs-Thomson law
Norbert Požár
Institute of Science and Engineering
Kanazawa University
Abstract
We present a simple numerical method for the Stefan problem with the GibbsThomson law based on a coupling of a total variation minimization algorithm andthe finite element method for the heat equation, and discuss a few numerical results.

1

Introduction

In this note we present a simple numerical method for solving a Stefan-type problem with
the Gibbs-Thomson law to model a growth of small crystals. Our goal is to explore thefeasibility of coupling the algorithm due to Oberman, Osher, Takei and Tsai [OOTT]
to solve the crystalline mean curvature flow problem without any regularization witha finite element method for the heat equation. Our main motivation are the recent
advances in the theory of the crystalline mean curvature flow in an arbitrary dimension[CMP, CMNP, GP1, GP2].
In a dimension n ∈ N, the problem is to find a relatively open set E ⊂ Rn × [0, ∞)and a function u : Rn × [0, ∞) → R that satisfy
in Rn × (0, ∞) \ ∂E, ut = div(k∇u)(1.1)
LV = (−k∇u)out · ν − (−k∇u)in · ν on ∂E,
V = β(ν)(ακσ − u)on ∂E,
with appropriate initial data u0 and E0 . The set Et := {x : (x, t) ∈ E} represents theshape of an evolving crystal at the time t ≥ 0, and V and ν are respectively the outernormal velocity and the spatial outer unit normal vector of its boundary ∂Et . u represents
the temperature, and the heat conductivity k is a given function that can depend on u as
well. (·)out and (·)in denote the limits from the outside and the inside of E respectively.
Note that a negative temperature on the crystal surface will force the crystal to grow,which will in turn release the latent heat of phase transition. The crystalline mean curvature κσ of the surface ∂Et will be explained below. The function β : S n−1 → (0, ∞) is agiven mobility on the unit sphere, and α, L > 0 are fixed constants.
Problem (1.1) is a model of a phase transition with the so-called kinetic undercooling
represented by the Gibbs-Thomson law in the third equation. The second equation states
the conservation of energy during the phase transition, where the amount of energy needed

to change the phase per unit volume is given by the specific latent heat of phase transition
L, and this energy is delivered by the heat flux −k∇u. This problem appears in models
of dendritic growth and solidification. For a more detailed discussion, see [CMOS, BGN]
and the survey article [Ga] and the references therein.

2

Crystalline mean curvature

Following [AG, T], the crystalline mean curvature can be introduced as the first variationof the anisotropic surface energy functional defined for any suﬃciently smooth set U ⊂ Rn
as∫
F(U ) :=
σ(ν) dS,∂U

where σ : S n−1 → (0, ∞) is a given anisotropy. For the detailed overview of this topic,
see the survey [B] or [Gu, GG, GP1, CMP]. In what follows we extend σ positively onehomogeneously to whole Rn as{p), p ̸= 0,|p|σ( |p|σ(p) =
0,
p = 0.
We will assume that the extended σ is convex.
The anisotropy σ determines the optimal crystal shape (Wulﬀ shape), that is, theshape that minimizes the anisotropic surface energy among shapes with the same volume
asW := {x : x · p ≤ σ(p), p ∈ Rn },(2.1)

see [T].
If σ is smooth on Rn \{0} and {p : σ(p) ≤ 1} is a strictly convex set, and ∂U is smooth,the first variation of F with respect to a change of volume is div∂U (∇σ(ν)), where div∂U
is the surface divergence, see [B]. We defineκσ := − div∂U (∇σ(ν)).

(2.2)

We are, however, specifically interested in anisotropies that are not smooth, in particular whose level set {σ ≤ 1} is a convex polytope and W above is the dual convex polytope[R]. In other words, when σ : Rn → [0, ∞) is a convex, piece-wise linear function. We
call such anisotropies crystalline. In this case, the definition of κσ becomes much more
involved since ∇σ(ν) is no longer defined pointwise and might be discontinuous even on
a smooth surface. A natural generalization of (2.2), consistent with the general abstract
theory of monotone operators due to Kōmura and Brézis [K, Br], is to replace the firstvariation of F with a subdiﬀerential on an appropriate Hilbert space. This leads to thedefinition [B, CMP, GP1]
κσ := − div∂U zmin ,where zmin is the element of the set{}z ∈ L∞ (∂U ) : z ∈ ∂σ(ν), div z ∈ L2 (∂U )(2.3)

that minimizes ∥div z∥L2 (∂U ) , if the set is nonempty. Here ∂σ is the subdiﬀerential of σ,∂σ(p) := {ξ ∈ Rn : σ(p + h) − σ(p) ≥ ξ · h, h ∈ Rn }.
While zmin might not be unique, div zmin is. Unfortunately, even for smooth ∂U the aboveset (2.3) might be empty, in which case κσ is not defined.
We call such κσ the crystalline mean curvature of ∂U . This quantity enjoys a numberof interesting properties, including the comparison principle. Furthermore, it is a singular,
nonlocal quantity on the surface ∂U . In particular, if ∂U has a flat part, called a facet,with normal ν such that σ is not diﬀerentiable at ν, κσ then depends on the shape of the
facet. For example, in dimension n = 2, facets of U are line segments, and on each ofthem κσ is a constant that is inversely proportional to their length.
The well-posedness of (1.1) with a crystalline σ seems to be open. A general theoryof solutions for the crystalline mean curvature flow (3.1) has been available only recently[CMP, CMNP, GP1, GP2].

3

Numerical algorithm

There have been many very successful algorithms developed to solve (1.1) with a smooth
σ numerically and listing them all is beyond the scope of this note. One of the diﬃculties
is the unstable nature of the dendritic growth and an extra care must be taken to avoid
the discretization and the choice of a mesh to produce unwanted artifacts. See [CMOS,
BGN, Ga] and the references therein. In this work to treat the crystalline mean curvature
directly without any regularization, we instead use the algorithm proposed by [OOTT]
that relies on the split Bregman iteration scheme introduced in [GO] to eﬃciently findminimizers of anisotropic total variation functionals. This approach builds on the level set
formulation due to Chambolle [C] of the minimizing movements time semidiscretizationof the anisotropic mean curvature flow
V = β(ν)(κσ + f )

(3.1)

on ∂E

by Almgren, Taylor and Wang [ATW].
In this algorithm, we fix a computational domain Ω ⊂ Rn suﬃciently large so thatEt ⊂ Ω for all t ∈ [0, T ] for some fixed T > 0, and choose a time step h > 0. Then
we approximate Etm , tm = mh, m = 1, 2, · · · , by a sequence of sets Em ⊂ Ω, where
Em := {x : vm (x) < 0}, and vm ∈ L2 (Ω) ∩ BV (Ω) is the minimizer of the functional∫
12
Jm (v) :=∥v − wm−1 ∥L2 + σ(∇v) dx − ⟨v, fm ⟩L2 .
2hΩ
Here wm−1 is the signed distance function to the set Em−1 with respect to the anisotropyβ,wm−1 (x) = signdistβ Em−1 (x) :=where β ◦ is the polar of β [R],◦β (x) := sup{infy∈Em−1β ◦ (x − y) −

}
x·ν
n−1
:ν∈S
.
β(ν)

infc

y∈Em−1β ◦ (y − x),The minimizer of Jm can be found eﬃciently by the split Bregman minimization
method of [GO] as observed in [OOTT]. We proceed by introducing a new variable d :Rn → Rn and add a constraint d = ∇v that is then relaxed by a L2 -penalization term.
In other words, we choose λ > 0 and in place of Jm we minimize
∫
1
λ
2
˜
Jm (v, d) :=∥v − wm−1 ∥L2 + σ(∇v) dx − ⟨v, fm ⟩L2 + ∥d − ∇v∥2L2 .
2h2ΩSince all the terms are convex, we can attempt to find the minimizer by alternatively
minimizing with respect to v and d. This produces a sequence converging to the uniqueminimizer of J˜m . However, this minimizer does not generally satisfy d = ∇v and therefore
is not a minimizer of Jm .
This is addressed by introducing a third variable b : Rn → Rn that accumulates the“error” ∇v − d during the iteration. The full algorithm can be stated as follows: Set
bm,0 = dm,0 = 0 and for k = 0, 1, . . . iterate1λ∥v − wm−1 ∥2 − ⟨v, fm ⟩ + ∥dm,k − ∇v − bm,k ∥2 ,
2h
2
v
∫
λ
dm,k+1 ← arg min σ(d) dx + ∥d − ∇vm,k − bm,k ∥2 ,2dΩ
bm,k+1 ← bm,k + ∇vm,k+1 − dm,k+1 ,vm,k+1 ← arg min(3.2a)(3.2b)(3.2c)

until convergence is achieved, for instance, until ∥vm,k+1 − vm,k ∥ is suﬃciently small.
This algorithm converges to a stationary point, see [GO]. When that happens, we
deduce from the third step (3.2c) that the limit (vm , dm , bm ) satisfies ∇vm = dm andtherefore vm is the unique minimizer of Jm . In particular, it does not depend on the
choice of λ. However, λ influences the diﬃculty to solve the minimization problem (3.2a)and the overall speed of convergence.
Note that solving (3.2a) is equivalent to solving a linear elliptic partial diﬀerential
equation for v, while solving (3.2b) after discretization yields a completely decoupled
minimization problem at each node, which only requires evaluating the projection on the
optimal shape W in (2.1). See [GO, OOTT] for more details.
To couple this scheme with the heat equation in (1.1), we interpret the second equationin (1.1) as an energy source concentrated on the surface ∂E and add it to the heat equation:
ut = div(k∇u) + LV Hn−1 ⌊∂Et .
After time discretization, for each m = 1, 2, . . ., we
1. solve for Em using the above algorithm with fm = −α−1 um−1 (in the TV minimization, we take β̃ = αβ with β from (1.1)), and then
2. solve for the solution um of
um − h div(k∇um ) = um−1 + hLVm Hn−1 ⌊∂Em ,(3.3)

where the normal velocity Vm is estimated using the value of the signed distancefunction signdist Em−1 with respect to the usual Euclidean metric on the set ∂EmashVm = signdist Em−1 .

We discretize both (3.2) and (3.3) by the standard piece-wise linear finite elements
on a regular triangular grid in two dimensions or a regular tetrahedral grid in three
dimensions [P]. The discrete signed distance function is computed using the fast sweeping
algorithm [Z] with the initialization scheme proposed in [P].
Recomputing the signed distance function necessarily moves the level set and introduces an error. However, we observe that vm is still very close to a distance function of
Em near its boundary and therefore we usually set wm = vm . We recompute the actual
signed distance function using the fast sweeping algorithm only after a prescribed numberof steps.
It would be more appropriate to set the value of fm using the extension of the valueof −um−1 |∂Em−1 via a transport equation as in [AS] to increase the accuracy and reduce
the need to reinitialize the distance function. This will be explored in the future.
It is natural to consider an adaptive mesh or perform the computation only near theboundary of the evolving domain. This is a subject of ongoing work.

4

Numerical results

In this section we present results of a few simple numerical experiments using the above
algorithm, see Figures 1–4.
We fix the domain Ω = (− 12 , 21 )2 for the} minimization and the fast
{ total2 variation
1
sweeping method, and the domain Ωheat = x ∈ R : |x| < 2 for the heat equation. The
boundary data for the heat equation is set to −10−3 . For the mobility, we always use
β(p) = 100|p|2 . The conductivity is taken to be 1 outside the crystal and 0 inside.
We discretize the domain by subdividing it uniformly into 2562 squares and split eachof those in half to produce the triangular mesh. The time step is fixed at τ = 2.5 × 10−4 .
Larger time steps lead to instability, most likely due to the explicit nature of the couplingof the heat equation with the anisotropic mean curvature flow.
The optimal choice of the parameter λ is an interesting problem. In the original paper[GO] for total variation denoising applications, the value λ = h2 was used. However, in the1
was found using numerical experiments to provide
current computations the value λ = 10h
a significant speed-up of the convergence. With this choice, only about 3–5 iterations of
the algorithm (3.2) per time step are needed, and each of the presented results takes 15–30
minutes to compute on an Intel Core i7-4770K processor.
The simple numerical results show that the method can reproduce facet breaking due
to the uneven cooling of the facets, and that neighboring crystals do not merge because
the diﬀusion is limited inside the crystal. Due to the technical limitations of the current
implementation, the fixed boundary for the heat equation is relatively close to the growing
crystal, and this leads to a strong forcing at the exposed vertices. The facets there are
therefore short and the crystal surface is not fully faceted. On the other hand, thisillustrates that the method is capable of handling not fully faceted crystals, which isobserved in practice in such strongly forced situations.
In summary, the initial implementation produces promising results, but more work isneeded to make it competitive with other available algorithms. Adaptive meshes, a moreaccurate treatment of the boundary forcing term in the level set method and various otherimprovements are being considered. In particular, a three dimensional simulation is inpreparation.

Figure 1: The evolution of a single square crystal of initial side length 0.05 growing with
α = 5 × 10−7 , plotted at 2 second intervals aligned with the mesh (left) and rotated by
0.2 radians (right).

Figure 2: The evolution of a single square crystal of initial side length 0.05 growing with
α = 10−8 , plotted at 2 second intervals aligned with the mesh (left) and rotated by 0.2
radians (right). Since the curvature eﬀect is significantly weaker than in Figure 1, only
short facets appear in regions with a relatively weak forcing. The eﬀect of the mesh ismore pronounced: the dendrites splitting in half along the diagonal directions are almost
certainly mesh artifacts.

Figure 3: The evolution of a single hexagonal crystal of initial side length 0.05 growing
with α = 2.5 × 10−7 (left) and α = 10−8 (right), plotted at 2 second intervals. The
horizontal dendrites grow faster, indicating an anisotropy caused by the mesh.

Figure 4: The evolution of three initially hexagonal or square crystals growing with α =5 × 10−7 , plotted at 2 second intervals. Note that the crystals do not merge.

Acknowledgments
The work of the author is partially supported by JSPS KAKENHI Grant No. 26800068(Wakate B).

References[AS] D. Adalsteinsson and J. A. Sethian, The fast construction of extension velocities in level set
methods, J. Comput. Phys. 148 (1999), no. 1, 2–22, DOI 10.1006/jcph.1998.6090.
[ATW] F. Almgren, J. E. Taylor, and L. Wang, Curvature-driven flows: a variational approach, SIAM
J. Control Optim. 31 (1993), no. 2, 387–438, DOI 10.1137/0331020.
[AG] S. Angenent and M. E. Gurtin, Multiphase thermomechanics with interfacial structure. II. Evolution of an isothermal interface, Arch. Rational Mech. Anal. 108 (1989), no. 4, 323–391, DOI
10.1007/BF01041068.
[BGN] J. W. Barrett, H. Garcke, and R. Nürnberg, Finite-element approximation of one-sided Stefan problems with anisotropic, approximately crystalline, Gibbs-Thomson law, Adv. Diﬀerential
Equations 18 (2013), no. 3-4, 383–432.
[B] G. Bellettini, An introduction to anisotropic and crystalline mean curvature flow, Hokkaido
Univ. Tech. Rep. Ser. in Math. 145 (2010), 102–162.
[Br] H. Brézis, Monotonicity methods in Hilbert spaces and some applications to nonlinear partial
diﬀerential equations, Contributions to nonlinear functional analysis (Proc. Sympos., Math. Res.
Center, Univ. Wisconsin, Madison, Wis., 1971), Academic Press, New York, 1971, pp. 101–156.
[C] A. Chambolle, An algorithm for mean curvature motion, Interfaces Free Bound. 6 (2004), no. 2,
195–218, DOI 10.4171/IFB/97.
[CMNP] A. Chambolle, M. Morini, M. Novaga, and M. Ponsiglione, Existence and uniqueness foranisotropic and crystalline mean curvature flows, preprint, available at https://arxiv.org/
abs/1702.03094.
[CMP] A. Chambolle, M. Morini, and M. Ponsiglione, Existence and Uniqueness for a Crystalline Mean
Curvature Flow, Comm. Pure Appl. Math. 70 (2017), no. 6, DOI 10.1002/cpa.21668.
[CMOS] S. Chen, B. Merriman, S. Osher, and P. Smereka, A simple level set method for solving Stefan
problems, J. Comput. Phys. 135 (1997), no. 1, 8–29, DOI 10.1006/jcph.1997.5721. MR1461705[Ga] H. Garcke, Curvature driven interface evolution, Jahresber. Dtsch. Math.-Ver. 115 (2013), no. 2,
63–100, DOI 10.1365/s13291-013-0066-2.
[GG] M.-H. Giga and Y. Giga, Evolving graphs by singular weighted curvature, Arch. Rational Mech.
Anal. 141 (1998), no. 2, 117–198.
[GP1] Y. Giga and N. Požár, A level set crystalline mean curvature flow of surfaces, Adv. Diﬀerential
Equations 21 (2016), no. 7-8, 631–698.
[GP2], Approximation of general facets by admissible facets for anisotropic total variationenergies and its application to the crystalline mean curvature flow, Comm. Pure Appl. Math.
(to appear).

[GO] T. Goldstein and S. Osher, The split Bregman method for L1-regularized problems, SIAM J.
Imaging Sci. 2 (2009), no. 2, 323–343, DOI 10.1137/080725891.
[Gu] M. E. Gurtin, Thermomechanics of evolving phase boundaries in the plane, Oxford Mathematical
Monographs, The Clarendon Press, Oxford University Press, New York, 1993.
[K] Y. Kōmura, Nonlinear semi-groups in Hilbert space, J. Math. Soc. Japan 19 (1967), 493–507.
[P] N. Požár, On the self-similar solutions of the crystalline mean curvature flow in three dimensions(in preparation).

[OOTT] A. Oberman, S. Osher, R. Takei, and R. Tsai, Numerical methods for anisotropic mean curvature
flow based on a discrete time variational formulation, Commun. Math. Sci. 9 (2011), no. 3, 637–
662, DOI 10.4310/CMS.2011.v9.n3.a1.
[R] R. T. Rockafellar, Convex analysis, Princeton Mathematical Series, No. 28, Princeton University
Press, Princeton, N.J., 1970. MR0274683 (43 #445)[T] J. E. Taylor, Constructions and conjectures in crystalline nondiﬀerential geometry, Diﬀerential
geometry, Proceedings of the Conference on Diﬀerential Geometry, Rio de Janeiro (B. Lawson
and K. Tanenblat, eds.), Pitman Monogr. Surveys Pure Appl. Math., vol. 52, Longman Sci.
Tech., Harlow, 1991, pp. 321–336, DOI 10.1111/j.1439-0388.1991.tb00191.x.
[Z] H. Zhao, A fast sweeping method for eikonal equations, Math. Comp. 74 (2005), no. 250, 603–
627, DOI 10.1090/S0025-5718-04-01678-3.

Faculty of Mathematics and Physics
Institute of Science and Engineering
Kanazawa University
Kakuma town, Kanazawa
Ishikawa 920-1192
JAPANE-mail address: npozar@se.kanazawa-u.ac.jp