# Adjoint-based optimization of two-dimensional Stefan problems

Tomas Fullana, Vincent Le Chenadec, Taraneh Sayadi

► To cite this version:

Tomas Fullana, Vincent Le Chenadec, Taraneh Sayadi. Adjoint-based optimization of two-dimensional Stefan problems. Journal of Computational Physics, 2023, 475, pp.111875. 10.1016/j.jcp.2022.111875. hal-03982348

HAL Id: hal-03982348

https://hal.science/hal-03982348

Submitted on 10 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Adjoint-based optimization of two-dimensional Stefan problems

Tomas Fullana<sup>a,*</sup>, Vincent Le Chenadec<sup>b</sup>, Taraneh Sayadi<sup>a,c</sup>

<sup>a</sup>Institut Jean le Rond d'Alembert, Sorbonne Université/CNRS, F-75005 Paris, France

<sup>b</sup>MSME, Université Gustave Eiffel, F-77454 Marne-la-Vallée, France

<sup>c</sup>Institute of Combustion Technologies (ITV), RWTH-Aachen University, Aachen, Germany

**Abstract**

A range of optimization cases of two-dimensional Stefan problems, solved using a tracking-type cost-functional, is presented. A level set method is used to capture the interface between the liquid and solid phases and an immersed boundary (cut cell) method coupled with an implicit time-advancement scheme is employed to solve the heat equation. A conservative implicit-explicit scheme is then used for solving the level set transport equation. The resulting numerical framework is validated with respect to existing analytical solutions of the forward Stefan problem. An adjoint-based algorithm is then employed to efficiently compute the gradient used in the optimisation algorithm (L-BFGS). The algorithm follows a continuous adjoint framework, where adjoint equations are formally derived using shape calculus and transport theorems. A wide range of control objectives are presented, and the results show that using parameterised boundary actuation leads to effective control strategies in order to suppress interfacial instabilities or to maintain a desired crystal shape.

**Keywords:** Stefan problem, level set, Cut Cell method, gradient-based optimization, continuous adjoint

**1. Introduction**

Stefan problem, named after J. Stefan [1] due to his substantial contribution to research on moving and free boundaries, model transport and transfer phenomena, in particular solid-liquid phase change in evaporating or chemically reacting flows. Such phenomena govern the interface motion in many engineering related problems such as dendritic solidification [2, 3], phase transformation in metallic alloys [4], and solid fuel combustion [5].

In the applications of interest to this study, the Stefan condition arises from the interaction of liquid and solid phases (both considered incompressible), resulting in a moving liquid-solid interface (freezing or melting front). The speed of the front is directly related to the jump in the conductive heat flux across the interface. In one dimension, this problem has been studied in depth using various numerical algorithms [6, 7, 8]. In higher dimensions, various methods have been used such as the level set method in Limare *et al.* [9] and

<footnote>*Corresponding author
Email address: tomas.fullana@sorbonne-universite.fr (Tomas Fullana)</footnote>
Osher *et al.* [2] and front-tracking method in Juric *et al.* [3]. One of the main challenges associated with modeling such problems in multiple dimensions is due to the unstable dentritic pattern formation [10, 11, 12]. In crystal growth, for example, under-cooling triggers an instability mechanism, causing the solid phase of the material to grow into the liquid phase in a finger-like or dendritic fashion, resulting in complex interfacial shapes, which are challenging to predict numerically. In addition, parameters such as anisotropy and surface tension or curvature effect (Gibbs-Thomson) are shown to have a large impact on the dendritic shape of the crystal, which in turn need to be modeled accurately for the simulations to remain predictive. In this work, we present a general framework for tracking and modeling crystal growth in the presence of curvature effects. This algorithm then serves as a vessel to materialise the second, and main, objective of this study which is extracting optimisation strategies to control the resulting solidification process.

The shape of the interface strongly affects the outcome and time-frame of the production process in many industrial applications involving phase change. As a result, while predicting and modeling the resulting solidification process remains at the forefront of many research areas, it is as desirable to extract efficient control strategies to manipulate the motion of the interface, for instance, by tracking a prescribed trajectory. Two major types of optimization methods in use today are (i) gradient-based, and (ii) derivative free methods. While an efficient class of generic algorithms (belonging to the class of derivative free methods) based on the surrogate management framework [13] and artificial neural networks [14] have been used for optimization in fluid mechanics, mainly in the area of aerodynamic shape optimization, they could require many function evaluations, for training purposes for example. When detailed simulations of interfacial flows are concerned, each function evaluation commands a full (potentially unsteady) CFD computation, causing gradient-based methods to be at an advantage. Most common methods in extracting the gradient information, on the other hand, are analytical or use finite differences, neither of which are suited to the configuration of interest to this work. Adjoint-based algorithms present a suitable alternative, as they allow the determination of the gradient at a cost comparable to a single function evaluation [15]. The use of adjoint methods for design and optimization has been an active area of research which started with the pioneering work of Pironneau [16] with applications in fluid mechanics, and has been extensively used in aeronautical shape optimization by Jameson and co-workers [17, 18]. Ever since these groundbreaking studies, adjoint-based methods have been widely used in fluid mechanics particularly in the areas of aero- and thermo-acoustics [19, 20]. More recently flow regimes dominated by nonlinear dynamics, such as separation and mixing have also been analysed using adjoint-based techniques [21, 22, 23, 24]. Adjoint-based methods have also been employed for the purpose of sensitivity analysis or control in flows in the presence of large gradients (flames or interfaces) [25, 26, 27, 28], showing great promise, and therefore are adopted here to carry out the optimisation procedure.

In the context of Stefan problems various control strategies have been employed to track the location of the interface. In a one-dimensional setting, for example, set-valued fixed point equations [29] or linear-quadratic defect minimization [30] have been used to control the location of the front. Adjoint-based algorithms have
also been applied to a Stefan problem caused by heterogeneous reactions on a surface of a one-dimensional solid particle [5] to extract sensitives with respect to various kinetic parameters. Alternatively, in a two-dimensional setting, adjoint-based algorithms have been utilised previously together with finite element and finite difference approaches to track and control the location of interface by imposing heat flux (or temperature) at the boundary in order to realize the desired interface motion [31, 32, 33]. In particular, Bernauer & Herzog [34], making use of shape calculus tools, derived the set of adjoint equations to extract control strategies for Stefan problems with a sharp representation of the interface.

Similar to the approach of [34] shape calculus tools have been employed in this work to extract the corresponding adjoint equations. However, contrary to the previous studies, control strategies are extracted here to suppress instabilities of the solidification process (dentritic formation). The addition of curvature effects on the interface (Gibbs-Thomson relation), and the complex shapes encountered during the growth of the crystal require dedicated numerical algorithms capable of modeling both the highly nonlinear forward problem (correct representation of the resulting interface) and the resulting adjoint equations. In addition, while previous studies mostly concentrated on actuation by imposing temperature or heat flux at the boundaries, alternative control strategies using the length of interface or surface tension coefficient are also investigated here, to identify the most relevant and effective (numerical) control strategies in the context of solidification problems.

The paper is organised as follows. Sec. 2 describes the equations governing the forward Stefan problem. The details of the numerical discretization and various schemes used to solve for the temperature field inside the two phases and to track the interface in time are described in Sec. 3, and validation cases are discussed in Sec. 4. The adjoint equations are presented in Sec. 5 and results of the optimisation performed on selected cases are offered and discussed in Sec. 6. Summary and perspectives of the work are then discussed in Sec. 7.

## 2. Governing equations in continuous form

The problem of interest is the Stefan problem in the presence of two immiscible phases with matching densities (one liquid and the other solid). The state of the problem is characterized by the temperature distributions in either phase, as well as the boundary position between them. The computational domain $\Omega$ is partitioned into the time-dependent subdomains $\Omega_1(t)$ and $\Omega_2(t)$ occupied by the liquid (1) and solid (2) phases, respectively. The external boundary of the domain, denoted $\partial\Omega$, is fixed whereas the interface separating both phases $\Gamma(t) = \bar{\Omega}_1(t) \cap \bar{\Omega}_2(t)$ evolves in time. A schematic of this configuration is shown in Fig. 1.

Let $T_i: (x, t) \to \mathbb{R}^+$ denote the temperature field in either phase ($i = 1, 2$). When the densities $\rho_1$ and $\rho_2$ are equal, convective heat transfers vanish. In addition, when the background pressure is assumed constant,
<figure><img src="image_1.png" /><figcaption>Figure 1: Schematic of the two-phase Stefan problem</figcaption></figure>

the heat transport equation simplifies to,

$$ \forall i \in [1, 2], \quad \rho_i c_i \frac{\partial T_i}{\partial t} = \nabla \cdot (k_i \nabla T_i), \quad t > 0, \quad x \in \Omega_i \qquad (2.1) $$

where $\rho_i$, $c_i$ and $k_i$ denote the density, the specific heat capacity at constant pressure, and the thermal conductivity, for each phase. Along the interface, energy balance states that

$$ V = -\lambda^{-1} \left( k_1 \frac{\partial T_1}{\partial n} - k_2 \frac{\partial T_2}{\partial n} \right), \quad x \in \Gamma \qquad (2.2) $$

where $\lambda$ denotes the latent heat of solidification and vector $n$ is the outward normal vector at the interface. This jump is taken from phase 1 to phase 2 with $\partial T_i/\partial n$ denoting the normal component of the temperature gradient in phase $i$. Eq. 2.2 is commonly referred to as the Stefan condition.

In classical Stefan problem, the temperature is set to $T_D(x, t) = T_M$ at the interface where $T_M$ is a constant equal to the melting temperature of the material. For problems involving crystal growth however [35], surface tension effects must be added to the thermodynamic boundary condition by introducing a dependence in the curvature at the front. For that purpose, we use the classical Gibbs-Thomson relation, as defined in Chen et al. [36]

$$ T_D(x, t) = T_M - \epsilon_V V - \epsilon_\kappa \kappa, \quad x \in \Gamma \qquad (2.3) $$

where $\kappa$ denotes the curvature at the interface (positive if the center of curvature lies in the solid phase), $V$ the velocity of the interface, $\epsilon_\kappa$ the surface tension coefficient, and $\epsilon_V$ the molecular kinetic coefficient. Unless stated otherwise, both $\epsilon_\kappa$ and $\epsilon_V$ are considered to be constants and the heat capacities, thermal conductivities and latent heat are all set to unity. In addition, to ease the notation, the jump in gradient of temperature is denoted as $[\nabla T]_2^1$ (also defined in Eq. 3.7).

When dealing with the numerical approximation of interfacial flows, two classes of methods are commonly used to represent the interface, namely (i) Lagrangian or "front-tracking" methods and (ii) the Eulerian or "front-capturing" methods. The former uses a parameterisation of the interface location (e.g. markers
or moving meshes), and has already been used in Stefan problems [37] but has not been adopted in this work due to the inherent difficulty of deriving the continuous adjoint equations with such methods (see Sec. 5). The latter can broadly speaking be divided into two categories: Volume-Of-Fluid (VOF) and Level Set methods. In adjoint-based optimization, the VOF method [25] may lead to numerical complications due to the piece-wise linear reconstruction of the interface. On the other hand, the Level Set method, where the interface is implicitly defined as a continuous function, has been proven to be well suited in continuous adjoint-based optimization, specifically for Stefan problems [34]. Moreover, this method has been shown to accurately reproduce dendritic pattern formation [36, 9]. An implicit signed distance strategy has therefore been used here to represent the interface.

A level set function $\phi: (x, t) \to \mathbb{R}$ is constructed, such that, at any time $t$, the front is equal to the zero level set of the function

$$ \Gamma(t) = \{x \in \Omega : \phi(x, t) = 0\}. \qquad (2.4) $$

The level set function is initially set to the signed distance function, with $d$ the distance to the front, such that

$$ \phi(x, 0) = \begin{cases} +d, & x \in \Omega_2, \\ 0, & x \in \Gamma, \\ -d, & x \in \Omega_1. \end{cases} \qquad (2.5) $$

The level set function is advected with a speed function $F$ which is extended from the normal speed $V$ at the front. The equation of motion of $\phi$ is given by the level set advection equation

$$ \frac{\partial \phi}{\partial t} + F |\nabla \phi| = 0. \qquad (2.6) $$

Equation (2.6) will move $\phi$ with the correct speed such that $\Gamma$ remains the zero level set of $\phi$. The field $\phi$ is also used to define the normal vector $n$ by

$$ n = \frac{\nabla \phi(x, t)}{|\nabla \phi(x, t)|}, \quad x \in \Gamma(t) \qquad (2.7) $$

and the curvature $\kappa$ by

$$ \kappa = \nabla \cdot \left( \frac{\nabla \phi(x, t)}{|\nabla \phi(x, t)|} \right), \quad x \in \Gamma(t). \qquad (2.8) $$

In two dimensions, the curvature $\kappa$ at the front is computed in a non-conservative form using

$$ \kappa = \frac{(\phi_y^2 \phi_{xx} - 2 \phi_x \phi_y \phi_{xy} + \phi_x^2 \phi_{yy})}{(\phi_x^2 + \phi_y^2)^{3/2}}. \qquad (2.9) $$

Considering the above mentioned equations and relations, the following generic forward two-phase Stefan problem (FP) can be recast in the level set framework given below, which is later used in the adjoint derivation.
Find a function $T: \Omega \times [0, t_f] \to \mathbb{R}$ and a function $\phi: \Omega \times [0, t_f] \to \mathbb{R}$ such that:

$$
\left\{
\begin{array}{ll}
\displaystyle \frac{\partial T_1}{\partial t} = \Delta T_1 & \text{in } \Omega_1(t) \\[1.5ex]
\displaystyle \frac{\partial T_2}{\partial t} = \Delta T_2 & \text{in } \Omega_2(t) \\[1.5ex]
T(x, 0) = T_0(x) & \text{in } \Omega \\[1.5ex]
\displaystyle \frac{\partial T(x, t)}{\partial n} = u(x, t) & \text{on } \partial\Omega \\[1.5ex]
T(x, t) = T_M - \epsilon_V V - \epsilon_\kappa \kappa & \text{on } \Gamma(t) \\[1.5ex]
\displaystyle \frac{\partial \phi}{\partial t} = -[\nabla T_i]_2^1 \cdot \nabla \phi & \text{on } \Gamma(t) \\[1.5ex]
\phi(x, 0) = \phi_0(x) & \text{in } \Omega
\end{array}
\right.
\qquad (\text{FP})
$$

Here, $T_i$ denote the restrictions of $T$ to $\Omega_i(t)$, and $T_0$ and $\phi_0$ are the initial conditions at $t=0$ of the temperature field and the level set function, respectively. The velocity coefficient in the Dirichlet boundary condition at the interface $\epsilon_V$ is set to zero.

3. Numerical framework

In this section, the methods used to solve both the forward and the adjoint Stefan problems are discussed. The solution procedure is split into six main steps: (i) solution of the heat equation in each phase, (ii) computation of the Stefan condition, (iii) extension of the velocity field normal to the interface, (iv) propagation of the level-set function, (v) reinitialisation of the level-set function (so as to satisfy the signed distance property), and finally (vi) handling of the dead and fresh cells. Each of these steps and their corresponding convergence properties are presented in the following subsections and summarised in Algorithm 1.

3.1. Temperature update

The heat equation (Eq. 2.1) is discretized on the subdomains Ω₁(t) and Ω₂(t) using a Cartesian grid (identical to the one used to solve the level set advection equation, described in detail in Sec. 3.4). Since the interface Γ(t) might not align with either of the grid axes, the discretization of the Laplacian operator is carried using a recently proposed Cut Cell Method [38]. The underlying idea is to modify the standard centered difference
**Algorithm 1:** Outline of the method

input : $T_0$, $T_D$, $\lambda$, $\kappa_1$, $\kappa_2$, $c_1$, $c_2$

output : $\phi$, $T$

Initialize the level set function $\phi_0$ and the temperature field $T_0$

repeat

1. Update the temperature field : $\tilde{T}_n \leftarrow T_{n-1}, T_D$ (Sec. 3.1)

2. Compute the Stefan condition : $V \leftarrow \lambda, \kappa_1, \kappa_2, c_1, c_2, [\nabla T]_2^1$ (Sec. 3.2)

3. Extend the velocity field : $F \leftarrow V$ (Sec. 3.3)

4. Update the level set function : $\phi_n \leftarrow \phi_{n-1}, F$ (Sec. 3.4)

5. Reinitialize the level set function : $\tilde{\phi}_n \leftarrow \phi_n$ (Sec. 3.5)

6. Clean or Initialize dead or fresh cells : $T_n \leftarrow \tilde{T}_n$ (Sec. 3.6)

until final time $t \leftarrow t_f$;

formulas for approximating the second derivatives to impose the proper boundary condition (Eq. 2.3) on the interface. Here, the procedure for imposing the Dirichlet boundary condition at the interface is briefly described (for details regarding a more general implementation of the operators, including their modification to account for Neumann boundary conditions on $\partial\Omega$, the reader is referred to [38]).

Following the Cut Cell strategy, the interface position is used to compute the areas (volume, respectively) wetted by each phase in each face (cell, respectively), henceforth referred to as face-centered surface (denoted $A_{i\pm1/2,j}$ and $A_{i,j\pm1/2}$) and cell-centered volume (denoted $V_{i,j}$) capacities. A piece-wise bi-quadratic interpolation of the level-set function values is used to determine the crossing points of the interface with the Cartesian grid, and supplied to the marching squares algorithm that infers the face and volume capacities of partially-filled faces and cells (Fig. 2). The construction of the operators is briefly described below for the $x$ contribution to the Laplacian, for Dirichlet boundary conditions.

In continuous form, the $x$ component of the temperature gradient may be integrated over a volume $\Omega$ filled by any given phase, and transformed according to Stokes' theorem to yield,

$$ \int_{\Omega} \frac{\partial T}{\partial x} dV = \oint_{\partial\Omega} T \mathbf{e}_x \cdot d\mathbf{S} \qquad (3.1) $$

where $d\mathbf{S}$ denotes the outward-pointing surface element and $\mathbf{e}_x$ the unit vector along the $x$ direction. If $\Omega = \Omega_{i+1/2,j}$ consists of the intersection of a phase domain and a $x$ face-centered computational cell $(i+1/2,j)$, its contour may be decomposed as $\partial\Omega = S_{i,j}^- \cup S_{i+1,j}^+ \cup S_{i+1/2,j-1/2}^- \cup S_{i+1/2,j+1/2}^+ \cup (\Omega \cap \Gamma)$, where the $\pm$ subscript denotes the direction the mesh-aligned surfaces point into. Inserted in Eq. 3.1, this
decomposition leads to,

$$g_{i+1/2,j} \equiv \int_{\Omega_{i+1/2,j}} \frac{\partial T}{\partial x} dV = \underbrace{\int_{S_{i+1,j}} T dS - \int_{S_{i,j}} T dS}_{\text{fluid}} + \underbrace{\int_{\Omega_{i+1/2,j} \cap \Gamma} T e_x \cdot dS}_{\text{boundary}} \quad (3.2)$$

where the right-hand side is decomposed into fluid and boundary contributions, respectively approximated as,

$$g_{i+1/2,j}^{\Omega} = \frac{A_{i+3/2,j} + A_{i+1/2,j}}{2} T_{i+1,j} - \frac{A_{i+1/2,j} + A_{i-1/2,j}}{2} T_{i,j}, \quad (3.3)$$

and,

$$g_{i+1/2,j}^{\Gamma} = \left( \frac{A_{i+3/2,j} - A_{i+1/2,j}}{2} + \frac{A_{i+1/2,j} - A_{i-1/2,j}}{2} \right) T_D. \quad (3.4)$$

To preserve symmetry, the discrete divergence is set to the negative transpose of the (fluid) gradient operator, which yields the following approximation of the temperature Laplacian,

$$\int_{\Omega_{i,j}} \Delta T dV \approx \frac{A_{i+1/2,j} + A_{i-1/2,j}}{2} \left( \frac{g_{i+1/2,j}}{V_{i+1/2,j}} - \frac{g_{i-1/2,j}}{V_{i-1/2,j}} \right) + \frac{A_{i,j+1/2} + A_{i,j-1/2}}{2} \left( \frac{g_{i,j+1/2}}{V_{i,j+1/2}} - \frac{g_{i,j-1/2}}{V_{i,j-1/2}} \right) \quad (3.5)$$

where the staggered volumes ($V_{i\pm1/2,j}$ and $V_{i,j\pm1/2}$) are interpolated from $V_{i,j}$.

<figure><img src="image_2.png" /><figcaption>Figure 2: Example of face capacities and staggered volumes necessary to compute the discrete homogeneous and inhomogeneous gradient operators in the positive $x$ direction (Equations 3.3 and 3.4).</figcaption></figure>

The values computed in the neighbouring points are projected back to the interface centroid in the partial cell using a bi-quadratic interpolation. Using the described discretisation and boundary conditions, the heat equation is then propagated in time using a Crank-Nicolson scheme.

The implementation of the Cut Cell method coupled with a Crank-Nicolson scheme is validated in the following stationary setup. A solid circle of radius $R = 0.75$ is initialized in a $2 \times 2$ domain. The initial temperature field is set to zero and a Dirichlet boundary condition is imposed at the interface $T_D = 1$. We solve only for the phase inside of the circle until a final time $t_f = 0.03125$. The simulation is performed for different
resolutions N = 32, 64, 128 with a fixed CFL = 0.5 corresponding to 16, 64, 256 iterations, respectively. In Fig. 3, the final temperature fields and the error in partial and full cells are shown. In stationary geometry, the Cut Cell method coupled with a Crank-Nicolson scheme shows an order of convergence slightly below 2 for both the partial and full cells.

<figure><img src="image_3.png" /><figcaption>Figure 3: Convergence study of the Cut Cell method in a stationary geometry with a Dirichlet boundary condition $T_D = 1$ imposed at the interface. The top figures show the position of the interface in red and the temperature field at final time for different resolutions. The bottom figure shows the convergence rate of the method in partial cells and in full cells.</figcaption></figure>

## 3.2. Stefan condition

As stated previously, the motion of the interface is solely dependent on the jump in the temperature gradient. It is therefore important to compute the normal gradient in temperature of each phase accurately, and to this end, the Johansen-Colella method [39] is used

$$ \nabla T|_{\Gamma} = \frac{1}{d_B - d_A} \left( \frac{d_B}{d_A}(T_D - T_A^*) - \frac{d_A}{d_B}(T_D - T_B^*) \right) \qquad (3.6) $$

with $T_D$ the dirichlet value imposed at the interface, $T_A^*$ and $T_B^*$ the interpolated values of the temperature field on points A and B respectively, and $d_A$ and $d_B$ the distances between the interface centroid to A and B respectively. The procedure for computing the gradient in one phase is as follows: (i) a shifted 3 × 3 stencil
<figure><img src="image_4.png" /><figcaption>Figure 4: (a) Schematic of gradient calculation. The Dirichlet value $T_d$ is imposed at the interface centroid in the partial cell and the temperatures $T_A$ and $T_B$ are determined via a quadratic interpolation from the neighboring 3 points in the vertical direction (dotted lines). (b) Convergence of the gradient computation.</figcaption></figure>

is chosen, as shown in Fig. 4(a), (ii) a line from the interface centroid is cast in the normal direction $n$, (iii) the crossing points $A$ and $B$ of this line and the vertical (or horizontal, depending on the normal orientation) segments of the neighboring 3 points are identified and the distances $d_A$ and $d_B$ are computed, and finally (iv) the values $T_A^*$ and $T_B^*$ are interpolated using $T_A^1, T_A^2, T_A^2$ and $T_B^1, T_B^2, T_B^2$, respectively. Once the normal gradient is computed in each phase, using Eq. 3.6, the jump is computed as

$$[\nabla T]_2^1 = \nabla T_1|_{\Gamma} - \nabla T_2|_{\Gamma}. \qquad (3.7)$$

The discrete velocities of the front in the partial cells are initialized with this jump and will be used as the boundary condition in the velocity extension procedure. To validate the method within our Cut Cell framework, we consider a stationary circle of radius $R = 0.5$ in a $1 \times 1$ domain and we initialize the temperature field with a similarity solution of the heat equation

$$T(r) = \begin{cases} T_{\infty} \left( 1 - \frac{F(r)}{F(R)} \right), & r > R \\ 0, & r < R \end{cases} \qquad (3.8)$$

with $T_{\infty} = -0.5$ a given under-cooling temperature, and $F(r) = E_1(1/4 r^2)$ where $E_1(t) = \int_x^\infty \frac{e^{-t}}{t} dt$. The discrete velocities at the front are computed for different resolutions $N = 32, 64, 128, 256$ and the resolution $256 \times 256$ is used as the reference solution for the convergence study. Fig. 4(b) shows that the gradient extraction procedure coupled with the interface location method results in close to second-order accuracy in both $L_2$ and $L_{\infty}$ norms.
### 3.3. Velocity extension

As highlighted in Sec. 2, the front velocity needs to be extended away from the interface. The most natural procedure is to let $V$ be a constant along the lines normal to $\Gamma$. To achieve this, the method described in Peng et al. is adopted here [40]. Using this approach, the velocity is extended in the normal direction by solving the following hyperbolic partial differential equation

$$ \left\{ \begin{array}{ll} \frac{\partial F}{\partial \tau} + S(\phi) \frac{\nabla \phi}{|\nabla \phi|} \cdot \nabla F = 0 & \text{in } \Omega \\ F(x, 0) = V & \text{on } \Gamma \end{array} \right. \qquad (3.9) $$

where $F$ is the extended velocity field equal to $V$ at the front, $\tau$ denotes a pseudo-time and $S(\phi)$ is the signature function

$$ S(\phi) = \begin{cases} -1 & \text{if } \phi < 0 \\ 0 & \text{if } \phi = 0 \\ +1 & \text{if } \phi > 0 \end{cases} \qquad (3.10) $$

Eq. 3.9 is then discretized using a first order upwind scheme and integrated in time by a forward Euler method until steady state. Taking $n$ as the normal vector defined as

$$ n = (n_x, n_y) = \left( \frac{\phi_x}{\sqrt{(\phi_x^2 + \phi_y^2)}}, \frac{\phi_y}{\sqrt{(\phi_x^2 + \phi_y^2)}} \right), \qquad (3.11) $$

the discretisation leads to

$$ F_{ij}^{n+1} = F_{ij}^{n} - \Delta\tau \left( (S_{ij}n_{ij}^{x})^{+} \frac{F_{ij} - F_{i-1j}}{h} + (S_{ij}n_{ij}^{x})^{-} \frac{F_{i+1j} - F_{ij}}{h} + (S_{ij}n_{ij}^{y})^{+} \frac{F_{ij} - F_{ij-1}}{h} + (S_{ij}n_{ij}^{y})^{-} \frac{F_{ij+1} - F_{ij}}{h} \right) \qquad (3.12) $$

where $h$ is the uniform grid spacing, $(x)^{+} = \max(0, x)$ and $(x)^{-} = \min(0, x)$, and the time step $\Delta\tau$ is chosen so that $\Delta\tau/h = 0.45$. The pseudo-time spawn in the velocity extension procedure is purely ficticious and the number of iterations in Eq. 3.12 corresponds to the width of the narrow-band (NB) around the 0-level set where the velocities are initialized. Fig. 5 shows an example of an initial velocity field for different narrow-band widths, after one iteration of the heat equation, where the temperature at the interface $T_D = \epsilon_\kappa\kappa$ depends only on the curvature. The computed discrete velocities are positive in the kinks and negative in the tips driving the initial crystal towards a circular shape.

### 3.4. Level Set advection equation

A semi-implicit scheme described in Mikula et al. [41, 42] is used here to solve the level set advection equation (Eq. 2.6). This scheme is briefly described below (for details of the scheme the reader is referred to [41]). Using this method, the equation is written in an alternative form

$$ \frac{\partial \phi}{\partial t} + F \frac{\nabla \phi}{|\nabla \phi|} \cdot \nabla \phi = 0 \qquad (3.13) $$
<figure><img src="image_5.png" /><figcaption>Figure 5: Velocity fields of a crystal-shape geometry. The red curves represent the interface location and the velocity fields correspond, from left to right, to narrow bands widths of 0, 6 and 12.</figcaption></figure>

which is then divided into conservative and non-conservative terms

$$ \frac{\partial \phi}{\partial t} + \underbrace{\nabla \cdot \left( F \phi \frac{\nabla \phi}{|\nabla \phi|} \right)}_{A} - \underbrace{\phi \nabla \cdot \left( F \frac{\nabla \phi}{|\nabla \phi|} \right)}_{B} = 0, \quad (3.14) $$

resulting in a second order partial differential equation akin to a weighted diffusion equation. The first term (A) has a diffusion coefficient $F\phi$ that depends on the solution and represents a nonlinear curvature flow whereas in the second term (B) the solution is multiplied by the curvature of its level-sets. The main idea behind Mikula et al. [41] scheme is to distinguish two cases: if the product $F\phi$ is negative (positive, respectively) then A represents a forward (backward, respectively) diffusion and B represents a backward (forward, respectively) diffusion. The forward diffusion is treated implicitly while the backward diffusion is treated explicitly leading to a semi-implicit scheme with a diffusive CFL number.

This scheme allows us to relax the CFL condition, usually present in most of level set methods. To check the robustness of the method for a $\Delta t/h^2$ ratio exceeding the usual CFL condition, we consider an initial level set function $\phi_0$ in a $1 \times 1$ domain, given by

$$ \phi_0(x, y) = \sqrt{x^2 + y^2} - R, \quad (3.15) $$

with $R = 0.8$ the radius of the 0-level set. The velocity field is set to $F = -1$ in the whole domain and we run the simulation until the final time $t_f = 0.3625$ for different resolutions $N = 16, 32, 64$ and different CFL numbers ranging from 1 to 16. The $L_2$ norm of the error with respect to the analytical solution is computed in the whole domain. Results in Fig. 6 show a second-order accuracy for any given CFL number for the retracting circle case.
<figure><img src="image_6.png" /><figcaption>Figure 6: The top figures represent the normalized error field for different CFL = 1, 4, 16 and fixed N = 64. The red (white, respectively) curve represents the initial (final, respectively) 0-level set. In the bottom figure the error in L₂ norm is plotted for different resolutions and CFL numbers. The vertical blue lines correspond to a fixed number of points for varying CFL numbers.</figcaption></figure>

### 3.5. Reinitialisation step

Depending on the case, usually after one or more iterations of the time advancement scheme, the level set function will cease to be an exact signed distance function, necessitating a reinitialisation step to enforce the level set function $\phi$ to remain an exact distance function from the front $\Gamma(t)$. Given a function $\phi_0$, which is not a signed distance function, it can be evolved into a function $\phi$ by solving the Eikonal equation

$$ \left\{ \begin{array}{ll} \frac{\partial \phi}{\partial \tau} = S(\phi_0)(1 - |\nabla \phi|) & \text{in } \Omega \\ \phi(x, \tau) = 0 & \text{on } \Gamma(t) \end{array} \right. \qquad (3.16) $$

where $\tau$ denotes a pseudo-time, $S(\phi_0)$ denotes a sign function, and $\phi$ a signed distance function (once steady-state is reached). There exist many numerical methods for solving the Eikonal equation and, here we have adopted that of Min et al. [43], also recently used in the work of Limare et al. [9]. This approach relies on a second order ENO spatial discretization with sub-cell resolution near the interface. The one-sided ENO
finite difference (in the x direction only) yields

$$ 
\begin{aligned}
D_x^+ \phi_{ij} &= \frac{\phi_{i+1,j} - \phi_{ij}}{h} - \frac{h}{2} \minmod (D_{xx}\phi_{ij}, D_{xx}\phi_{i+1,j}), \\
D_x^- \phi_{ij} &= \frac{\phi_{i,j} - \phi_{i-1,j}}{h} + \frac{h}{2} \minmod (D_{xx}\phi_{ij}, D_{xx}\phi_{i-1,j})
\end{aligned}
\qquad (3.17) $$

where $D_{xx}\phi_{ij} = (\phi_{i-1,j} - 2\phi_{ij} + \phi_{i+1,j})/h^2$ is the second order derivative of $\phi_{ij}$ and the “minmod” limiter is zero when both arguments have opposite signs, and takes the argument of smallest absolute value otherwise. The numerical Hamiltonian $|\phi|$ is computed as follows

$$ |\nabla \phi|_{ij} \approx H_G (D_x^+ \phi_{ij}, D_x^- \phi_{ij}, D_y^+ \phi_{ij}, D_y^- \phi_{ij}) \qquad (3.18) $$

where function $H_g$ is given by

$$ H_G(a, b, c, d) = \begin{cases} \sqrt{\max((a^-)^2, (b^+)^2) + \max((c^-)^2, (d^+)^2)} & \text{when } S(\phi^0) \ge 0, \\ \sqrt{\max((a^+)^2, (b^-)^2) + \max((c^+)^2, (d^-)^2)} & \text{when } S(\phi^0) < 0. \end{cases} \qquad (3.19) $$

Near the interface, the finite differences need to be modified in order to impose $\phi = 0$ where $\phi_0 = 0$. To this end, a quadratic ENO polynomial interpolation gives

$$ D_x^+ \phi_{ij} = \frac{0 - \phi_{ij}}{h^+} - \frac{h^+}{2} \minmod (D_{xx}\phi_{ij}, D_{xx}\phi_{i+1,j}) \qquad (3.20) $$

where

$$ \Delta h^+ = \begin{cases} h \left( \frac{1}{2} + \frac{\phi_{ij}^0 - \phi_{i+1,j}^0 - S(\phi_{ij}^0 - \phi_{i+1,j}^0) \sqrt{D}}{\phi_{xx}^0} \right) & \text{where } |\phi_{xx}^0| > \epsilon, \\ h \frac{\phi_{ij}^0 - \phi_{i+1,j}^0}{\phi_{ij}^0 - \phi_{i-1,j}^0} & \text{elsewhere} \end{cases} \qquad (3.21) $$

with $\phi_{xx}^0 = \minmod (\phi_{i-1,j}^0 - 2\phi_{ij}^0 + \phi_{i+1,j}^0, \phi_{ij}^0 - 2\phi_{i-1,j}^0 + \phi_{i-2,j}^0)$ and $D = (\phi_{xx}^0/2 - \phi_{ij}^0 - \phi_{i-1,j}^0)^2 - 4\phi_{ij}^0\phi_{i-1,j}^0$. The negative one-sided ENO difference $D_x^-$ is obtained similarly. A forward Euler scheme is then used for time advancement. This Hamiltonian extraction will also be used for check-pointing in the adjoint problem (Sec. 5) when needed. We validate the method by initializing a perturbed solution, similar to the test case found in [43], in a $4 \times 4$ domain

$$ \phi_0(x, y) = ((x-1)^2 + (y-1)^2 + 0.1) \times (\sqrt{(x \pm a)^2 + y^2} - R) \qquad (3.22) $$

where $R = 1$ and $a = 0.7$. Fig. 7 shows the initial perturbed level set function in a $128 \times 128$ grid converging towards a signed distance function with increasing number of iterations in time.

### 3.6. Dead and fresh cells

The last step of the method presented in this study is the treatment of fresh and dead cells. As the interface moves through the Cartesian grid, the temperature field needs to be initialized according to the front position.
<figure><img src="image_7.png" /><figcaption>Figure 7: The -2 to 1 level sets of the function $\phi$ are shown with a 0.1 step. The red curve represents the 0-level set. After 90 iterations the level set function is a true signed distance function.</figcaption></figure>

Two cases can be distinguished: (i) a full or partial cell becomes an empty cell, and (ii) an empty cell becomes a partial cell. In the first case, the previous temperature value is simply set to 0. In the second case, the temperature previously non-existing needs to be initialized. The procedure is similar to that of the Stefan condition, Sec. 3.2: (i) a shifted 3 × 3 stencil, as shown in Fig. 8, is chosen, (ii) a line from the interface centroid is cast in the opposite normal direction $-n$, (iii) the crossing points A and B of this line and the vertical (or horizontal depending on the normal orientation) segments of the neighboring 3 points, are identified, (iv) the values $T_A^*$ and $T_B^*$ are interpolated using $T_A^1, T_A^2, T_A^2$ and $T_B^1, T_B^2, T_B^2$ respectively, and (v) the coordinates $\mathbf{o} = (x_{new}, y_{new})$ of the barycentre of the cell to initialize is located. This last step is done by using the discrete face capacities defined in Sec. 3.1. An orthonormal coordinate system $\mathcal{R} = (\mathbf{o}, (x', n))$ similar to the parabola-fitted curvature found in [44] is then defined. Finally, the new temperature value, $T_{new}$, is linearly extrapolated by solving for $a_0$ and $a_1$ in the resulting system of equations

$$ \begin{cases} T_A^* = a_0x'_A + a_1 \\ T_B^* = a_0x'_B + a_1 \end{cases} \qquad (3.23) $$

where $a_1 = T_{new}$.

## 4. Validation and assessment

In this section, the numerical scheme previously described is validated and the effect of grid refinement, as well as, various physical parameters of the system, such as surface tension at the interface, on the resulting crystal growth are demonstrated and discussed.
<figure><img src="image_8.png" /><figcaption>Figure 8: Example of empty cells becoming partial cells from the point of view of the solid phase. The interface at time $t^n$ moves in the normal direction with speed $V$. At time $t^{n+1}$, the newly initialized value $T_{new}$ located at the partial cell centroid is extrapolated from $T_A$ and $T_B$.</figcaption></figure>

### 4.1. Frank's sphere

The growth of an ice sphere surrounded by an under-cooled liquid was initially studied by Frank [45], and the initial sphere radius was shown to evolve in a self-similar manner as the square-root of time. This problem can therefore be used to validate the accuracy of the numerical scheme implemented here.

In this configuration, the temperature field is given by

$$ T(r,t) = T(s) = \begin{cases} T_{\infty} \left( 1 - \frac{F(s)}{F(S)} \right), & s > S \\ 0, & s < S \end{cases} \qquad (4.1) $$

where $r = \sqrt{x^2 + y^2}$, $s = r/t^{1/2}$, $T_{\infty}$ is a given under-cooling temperature, and $F(s)$ denotes the similarity solution of the heat equation $F(s) = E_1(1/4 s^2)$. Numerical errors may lead to an unwanted alteration of the initial shape due to the Mullins-Sekerka instability [11], where a perturbed solution can lead to unstable dendritic growth in the case of zero melting temperature ($T_M = 0$). We therefore test the robustness of our method using the initial parameters recommended in Almgren et al. [46]. The initial level set function is set to a circle of radius $R_0 = 1.56$ in a $8 \times 8$ domain surrounded by an initial temperature field (Eq. 4.1) with $T_{\infty} = -0.5$. The initial time is set to 1 and the simulation is advanced until a final time $t_f = 2$. The results of this study are summarized in Fig.9. The intrinsic regularization of our method allows the level set function to retain its initial circular shape, avoiding spurious interface oscillations. The following test cases will focus
<figure><img src="image_9.png" /><figcaption>Figure 9: Convergence of the growing Frank’s sphere for different resolutions N = 32, 64, 128. The top left figure shows the convergence of the radius towards the analytical solution as a function of time. The bottom left figure represents the error of the final radius as a function of N. The right figure shows the position of the interface at times t = 0, 0.5, 1 for the different grids. The dotted lines correspond to the initial radius R₀ = 1.56 and the analytical radius at final time R(t = 1) = 2.206.</figcaption></figure>

instead on the dendritic growth of an initial crystal shape in a uniform under-cooled temperature field.

## 4.2. Dendritic growth

Crystal growth is an unstable phenomenon that occurs spontaneously in nature. Its appearance is the result of a competition between the natural growth due to the Stefan condition (Eq. 2.2), which tends to stretch the tips, and the surface tension effect present in the Gibbs-Thomson relation (Eq. 2.3), which tends to restore the flatness of the interface (see [35] for further details). Numerical reproduction of such patterns is a challenge, since the dendritic growth depends on the smallest resolved length scale. In order to validate our method, a case similar to that of Chen et al. [36] is considered, where an initial solid crystal is surrounded by an under-cooled liquid. The level set function is initialized in a 2 × 2 domain as

$$ \phi_0(x, y) = \sqrt{x^2 + y^2} \times (0.1 + 0.02 \cos(4\alpha) - 0.01), \quad (4.2) $$

where α is the angle of the interface with respect to the x axis, and the under-cooled temperature is $T_∞ = -0.5$.

At the first stage, the effect of the grid resolution on the solution (Fig. 10) is assessed. Surface tension coefficient is fixed at $\epsilon_κ = 0.0004$, and velocity coefficient is set to $\epsilon_V = 0$. The simulations are advanced with
<figure><img src="image_10.png" /><figcaption>Figure 10: Effect of grid refinement on the dendrite growing tips. The initial condition is a solid crystal at temperature $T = 0$ surrounded by an under-cooling temperature $T_{\infty} = -0.5$. The color map represents the final temperature field. The final time is $t_f = 0.5$ and the interface is plotted with a time step of 0.025.</figcaption></figure>

insulated boundary conditions until a final time $t_f = 0.5$ for different grid resolutions $N = 50, 100, 150, 200$. Fig. 10 shows that as the grid size decreases, the tip-splitting appears earlier and the length of the dendrites increases. The final shape converges towards an 8-fold symmetric crystal shape.

In the second case, the surface tension coefficient is varied from 0.0004 to 0.001 for a fixed $N = 200$. Fig. 11
clearly demonstrates the stabilizing effect of the Gibbs-Thomson relation. The tip-splitting disappears as the surface tension increases. This behaviour is explained by a stronger regularization due to higher surface tension coefficient, reducing the growth rate of the instability, and is in agreement with the observations in similar regimes [36].

<figure><img src="image_11.png" /><figcaption>Figure 11: Effect of the surface tension coefficient on the dendrite growing tips. The initial condition is a solid crystal at temperature $T = 0$ surrounded by an under-cooling temperature $T_\infty = -0.5$. The color map represents the final temperature field. The final time is $t_f = 0.6$ and the interface is plotted with a time step of 0.03.</figcaption></figure>
### 4.3. Anisotropy effects

In order to accurately reproduce the crystal shapes produced by dendritic growth, we introduce a variable surface tension coefficient, similar to [47] and [9], in the Gibbs-Thomson relation (Eq. 2.3) to account for anisotropy effects,

$$ \bar{\epsilon}_{\kappa}(\alpha) = \epsilon_{\kappa} \left( 1 + A \left[ \frac{8}{3} \sin^4 \left( \frac{1}{2} M (\alpha - \alpha_0) \right) - 1 \right] \right), \quad (4.3) $$

where $\epsilon_κ$ is the surface tension coefficient, A represents the weight of the anisotropy effect, M the mode number, α the angle of the interface with respect to the x-axis and $\alpha_0$ the prescribed angle of symmetry on which the dendrites will grow. To validate these effects, we initialize a six-fold crystal in an under-cooled liquid with $T_∞ = -0.8$ and a fixed velocity coefficient $\epsilon_V = 0$. The anisotropic weight is fixed at $A = 0.4$ and mode number $M = 6$ and the simulations are performed for two different prescribed angles $\alpha_0 = \pi/2, \pi/4$ until a final time $t_f = 0.09$ in a 2 × 2 domain with $N = 300$. Fig. 12 shows the crystals growing in the direction of the prescribed angle as well as secondary dendrites appearing from the main branches. When $\alpha_0$ is different from the initial crystal shape (Fig. 12(b)), a rotation of the primary dendrites towards that prescribed direction can be observed.

<figure><img src="image_12.png" /><figcaption>Figure 12: Anisotropy effects on the crystal growth for different prescribed angles $\alpha_0 = \pi/2, \pi/4$. The initial condition is a six-fold solid at temperature $T = 0$ surrounded by an under-cooling temperature $T_∞ = -0.8$. The color map represent the final temperature field. The final time is $t_f = 0.09$ and the interface is plotted with a time step of 0.0045.</figcaption></figure>
## 5. Optimisation algorithm

In this section, the optimization problem corresponding to the forward problem (FP) presented in Sec. 2 is presented, together with the respective cost functional and the corresponding adjoint problem, where the control variable is a Neumann boundary condition that acts on the domain boundary. The resulting optimisation algorithm, used throughout the various cases considered in this study, is also briefly described. The adjoint heat equation and the transport theorem used for the derivation are presented in detail in Appendix A. For the adjoint level set derivation, we refer the reader to Bernauer et al. [34].

We first define $T_{tf}$ and $\phi_{tf}$ as the desired temperature field and level set function respectively. The following tracking-type cost functional provides a mathematical description of the control goals stated above:

$$ \mathcal{J}(T, \phi, u) = \frac{\beta_1}{2} \int_{\Omega} |T(t_f) - T_{tf}|^2 + \frac{\beta_2}{2} \int_{\Omega} |\phi(t_f) - \phi_{tf}|^2 + \frac{\beta_3}{2} \int_{\Gamma} dS + \frac{\beta_4}{2} \int_{0}^{t_f} \int_{\partial\Omega} |u|^2. \quad (5.1) $$

The first two terms monitor the final temperature distribution and interface position and are used to initialize the adjoint fields. The third term is used in cases where the instabilities induce an increase in interface length while the final term models and penalises the control cost. The constants $\beta_1$ to $\beta_4$ act as weights for each term. We assume that each choice of the control $u$ leads to unique states $T(u)$ and $\phi(u)$. Therefore, the minimization problem (MP) reads,

$$ \begin{array}{ll} \min_{u} & \mathcal{J}(T(u), \phi(u), u) \\ \text{subject to} & (\text{FP}) \end{array} \qquad (\text{MP}) $$

Let $\theta$ be the adjoint temperature and $\psi$ the adjoint level set function. To compute the gradient of $\mathcal{J}(u)$, we introduce the Lagrange functional $\mathcal{L}$ (we omit $dx$, $dt$ and $ds$ for the sake of brevity),

$$ \begin{aligned} \mathcal{L}(T, \phi, u, \theta, \theta_D, \theta_I, \psi) = & \mathcal{J}(T, \phi, u) \\ & - \int_0^{t_f} \int_{\Omega_1} \left( \frac{\partial T_1}{\partial t} - \Delta T_1 \right) \theta_1 - \int_0^{t_f} \int_{\Omega_2} \left( \frac{\partial T_2}{\partial t} - \Delta T_2 \right) \theta_2 \\ & - \int_0^{t_f} \int_{\partial \Omega} (T - u) \theta_D \\ & - \int_0^{t_f} \int_{\Gamma} (T - (T_M - \epsilon_V V - \epsilon_\kappa \kappa)) \theta_I \\ & - \int_0^{t_f} \int_{\Gamma} \left( \frac{\partial \phi}{\partial t} + [\nabla T_i]_2^1 \cdot \nabla \phi \right) \psi \end{aligned} \quad (5.2) $$

The adjoint system is obtained by setting $\mathcal{L}_T(\cdot) = \mathcal{L}_\phi(\cdot) = 0$. The resulting adjoint Stefan problem (AP) reads,
Find a function $\theta: \Omega \times [t_f, 0] \to \mathbb{R}$ and a function $\psi: \Omega \times [t_f, 0] \to \mathbb{R}$ such that:

$$
\left\{
\begin{array}{ll}
-\displaystyle\frac{\partial\theta_1}{\partial t} = \Delta\theta_1 & \text{in } \Omega_1(t) \text{ (AP.a)} \\
\\
-\displaystyle\frac{\partial\theta_2}{\partial t} = \Delta\theta_2 & \text{in } \Omega_2(t) \text{ (AP.b)} \\
\\
\theta(x, t_f) = \beta_1(T(t_f) - T_{t_f}) & \text{in } \Omega \text{ (AP.c)} \\
\\
\displaystyle\frac{\partial\theta(x,t)}{\partial n} = 0 & \text{on } \partial\Omega \text{ (AP.d)} \qquad (\text{AP}) \\
\\
\theta(x,t) = \psi|\nabla\phi| & \text{on } \Gamma(t) \text{ (AP.e)} \\
\\
\displaystyle\frac{\partial\psi}{\partial t} + \text{div}(\psi\vec{V}) = \frac{1}{|\nabla\phi|} \frac{\partial T}{\partial n} [\nabla\theta_i]_2^1 \cdot n & \text{on } \Gamma(t) \text{ (AP.f)} \\
\\
\psi(x, t_f) = -\frac{\beta_2}{2} \left( \frac{\partial}{\partial n} |\phi_{t_f}|^2 + \kappa \left( |\phi_{t_f}|^2 + \frac{\beta_3}{\beta_2} \right) \right) & \text{in } \Omega \text{ (AP.g)} \\
\\
0 = \beta_4 u + \theta & \text{on } \partial\Omega \text{ (AP.h)}
\end{array}
\right.
$$

Equations AP.a and AP.b correspond to the heat equations with reverse time direction while Equation AP.c is used to initialize the adjoint temperature field in both phases. In Equation AP.e, the in-homogeneous Neumann boundary condition in $\partial\Omega$ existing in the forward problem (FP) is mapped to an homogeneous alternative. The Dirichlet boundary condition at the interface is defined in Equation AP.f where the adjoint temperature is equal to the adjoint level set value $\psi$ augmented by the Hamiltonian of the level set function $|\phi|$. In Equation AP.f, the adjoint level set advection equation is defined. The source term corresponds to the adjoint Stefan condition and takes into account the normal jump in gradient of adjoint temperature across the interface. The adjoint level set function is initialized using Equation (AP.g), and the final Equation (AP.h) represents the optimality condition, where the left-hand side represents the gradient used to update the control variable $u$.

Even thought it is known that gradient-based methods exhibit slow convergence, they are the method of choice to prove that the optimal control approach is reasonable due to their straight forward implementation. Therefore, we chose to solve the minimization problem (MP) by using the limited memory BFGS (L-BFGS) method, a quasi-Newton method originally described in [48]. The main characteristic of this method is that it determines the descent direction by preconditioning the gradient with an approximation of the Hessian matrix. This information is obtained using past approximations (the number of approximations is determined
by the memory length parameter which is set to m = 10 in all our optimization test cases) as well as the gradient. As an initial guess for the initial Hessian, we use the scaled identity matrix as described in [49]. The following algorithm summarizes the L-BFGS method used in our numerical examples,

<caption>**Algorithm 2:** Optimization procedure using L-BFGS method</caption>

**input** : $u_0$, $m = 10$
**output** : $u$, $T$, $\phi$, $\theta$, $\psi$
$k \leftarrow 0, l \leftarrow 0$
**while** not converged **do**
  Solve the forward Stefan problem (FP) for $T^k$ and $\phi^k$
  Solve the adjoint Stefan problem (AP) for $\theta^k$ and $\psi^k$
  Compute the gradient:
  $$ \nabla \mathcal{J}^k = \beta_4 u^k + \theta^k $$
  **if** $k \ge 1$ **then**
    $s^{k-1} = \sigma^{k-1} d^{k-1}$, $g^{k-1} = \nabla \mathcal{J}^k - \nabla \mathcal{J}^{k-1}$
    **if** $(s^{k-1})^T g^{k-1} \le 0$ **then**
      $l \leftarrow 0$
    **else if** $(s^{k-1})^T g^{k-1} > 0$ **then**
      $l \leftarrow l+1$
      **if** $l > m$ **then**
        $\mid$ Remove $\{s^{l-m}, g^{l-m}\}$
      **end**
      Add $\{s^{l-m}, s^{l-m}\}$
    **end**
  **end**
  Choose an initial approximation to the inverse of the Hessian $H_0^k$
  Construct the direction $d^k = -H^k \nabla \mathcal{J}^k$
  Determine $\sigma^k$ using a Line Search algorithm with backtracking where $\sigma^k = \argmin \mathcal{J}(u^k + \sigma^k d^k)$
  Update $u^{k+1} = u^k + \sigma^k d^k$
  $k \leftarrow k+1$
**end**
# 6. Results

A range of optimization cases for the Stefan problem with varying complexity are presented here. The results and performance of the algorithm MP is then analysed for these different setups. In all of the cases the desired temperature field, $T_{tf}$, and level set, $\phi_{tf}$, are computed beforehand and are used to drive the control variable $u$. A summary of the optimization parameters and optimization procedure results for each case is shown in Table 1 and the histories of the relative cost functionals $J/J^0$ is presented in Fig. 16.

*Implementation:* The code is written in Julia [50] a high level scientific programming language. The package *Flower.jl*, developed by the authors of this study, contains all the methods and test cases presented. We use *Iterative.jl* to solve the linear systems and *Optim.jl* [51] to build the optimization procedure.

## 6.1. Melting circle

In this case, an initial circle of radius $R = 0.75$ in a $2 \times 2$ domain is considered. The surface tension coefficient is constant and set to $\epsilon_k = 0.002$ while $\epsilon_V = 0$. The non-zero surface tension coefficient is added to further regularize the level set function. The number of points per dimension is $N = 64$ and the final time is set to $t_f = 0.1$. The control $u$ acts on the whole domain boundary and is of Dirichlet type. The adjoint boundary condition, on the other hand, remains a homogeneous Neumann boundary condition. The actuator $u$ is parameterized as

$$u = \sum_{p=1}^{2} a_p \cos(p\pi x) + \sum_{p=1}^{2} b_p \sin(p\pi x) \quad (6.1)$$

where $x = [-1, 1]$ corresponds to the bounds of the domain. Through the optimisation process the amplitude of each basis function, $[a_1, a_2, b_1, b_2]$ is prescribed using the gradient equation (AP.h). By opting for a parameterized distribution we ensure the smoothness of the actuation function. Due to the high sensitivity of the cost functional to the basis considered, with too many parameters creating multiple local minima, the number of parameters are kept at a low enough value to ensure the convexity of the problem while allowing spatial variation of the actuation function. Same strategy will be adopted in the following examples when parametrising the actuation functional.

The initial guess (ie. the starting point in our cost functional space) is set to $u = 0$ in order not to add a bias regarding the initial descent. Moreover, in this case, the term controlling the interface length $\beta_3$ is set to zero.

Fig. 13 shows the level set of the forward problem evolving towards the desired shape after a few iterations. The algorithm is able to recover both the retraction and the expansion of the initial circle.
<figure><img src="image_13.png" /><figcaption>Figure 13: Iterations 0, 3, 10 and 16 of the optimization procedure for the circle test case. The blue curve represents the desired shape and the red one the final position of the interface at a given iteration. The final time is $t_f = 0.1$ and the interface is plotted with a time step of 0.01. The color map corresponds to the error in temperature field. The inset shows the actuator u at a given iteration (red) and the desired one (blue).</figcaption></figure>

## 6.2. Mullins-Sekerka instability

In this second case, we consider a perturbed initial planar shape in a 2 × 2 domain such that the initial level set function is given by,

$$ \phi_0(x, y) = y + 0.6 + A \cos(2\pi x) \quad (6.2) $$

and the temperature field by,

$$ T_0(x, y) = \begin{cases} -1 + e^{-T_\infty \phi_0(x, y)}, & \phi_0 > 0 \\ 0, & \phi_0 < 0 \end{cases} \quad (6.3) $$
with the amplitude $A = 0.05$ and the under-cooling temperature $T_\infty = 1.2$. It is well known that an initial perturbation such as that prescribed in this problem leads to a Mullins-Sekerka type instability [11, 36] characterised by unstable dendritic growth. The surface tension and velocity coefficients are set to zero in this case. The purpose of this optimization test case is to reduce the instability by imposing an optimal actuation on the top boundary. The control variable $u$ is of Neumann type and we use $u = 0$ as the initial guess. Both surface tension and velocity coefficients are set to zero and we set the final time $t_f = 0.5$ with $N = 64$. Similar to the previous case the actuation is parametrised as

$$u = \sum_{p=1}^{4} a_p \cos(p\pi x) + \sum_{p=1}^{4} b_p \sin(p\pi x) \qquad (6.4)$$

In this case, we can add four extra basis function in the actuator parametrization and still recover the global minimum. The extra term in the cost functional (Eq. 5.1) which controls the length of the interface is also included and $\beta_3 = 0.1$. In Fig. 14, we can observe how the initial tip splitting is reduced as we go through the optimization procedure.

### 6.3. Anisotropic growing crystals with topology changes

In this last case, we consider three crystals asymmetrically disposed in an under-cooled liquid. The crystals will grow and eventually merge. The objective of this optimization is to drive the final shape towards the desired one by acting on the boundaries of the whole domain, and thereby to suppress the anisotropy effects. The surface tension coefficient is set to $\epsilon_\kappa = 0.0005$ and the velocity coefficient to $\epsilon_V = 0.002$. We choose a relatively small surface tension coefficient with respect to the velocity one in order to allow for strong dendritic formation and to examine the optimisation algorithm in a challenging case where the topology remains complex. Moreover, we add anisotropy effects by setting $\alpha_0 = \pi/2$ and $M = 4$ in equation4.3. The heat flux induced through actuation will have to compete with these effects in order to drive the interface towards the desired one. The simulations are run in a $4 \times 4$ domain with $N = 100$ until a final time $t_f = 0.45$. The under-cooling initial temperature is set to $T_\infty = -0.6$. The control variable $u$ is of Neumann type, with an initial guess $u = 0$. The following parametrisation is used,

$$u(x, p) = p_1((1 + \cos(\pi/8x))/2)^4 + p_2((1 + \sin(\pi/8x))/2)^4 \qquad (6.5)$$

In this final case, we need to restrict the basis used for the actuator (by fitting only two parameters) in order to descend towards the global minimum. Fig. 15 shows the evolution of the surface. It can be seen that initially the crystal is driven towards the domain corners by the anisotropy parameters. Through the optimisation procedure however, the final crystal shape tends towards the desired one. The topology changes are implicitly taken into account by our level set method.
<figure><img src="image_14.png" /><figcaption>Figure 14: Iterations 0, 2, 3 and 9 of the optimization procedure for the Mullins-Sekerka test case. The blue curve represents the desired shape and the red one the final position of the interface at a given iteration. The final time is $t_f = 0.5$ and the interface is plotted with a time step of 0.05. The color map corresponds to the error in temperature field. The inset shows the actuator $u$ at a given iteration (red) and the desired one (blue).</figcaption></figure>

## 7. Conclusion

In this study, we proposed a simple level set method coupled with a novel Cut Cell approach to tackle two-phase Stefan problems and an adjoint-based optimization procedure. The key features of the level set related algorithms are, (i) an implicit-explicit scheme to solve the level set advection allowing us to relax the usual CFL condition, (ii) the high order Johansen-Colella method used to compute the normal gradient across the interface, and (iii) the sub-cell resolution reinitialization procedure to retain the signed distance function property as the interface moves. Moreover, we use a novel Cut Cell method coupled with a Crank-Nicolson time integrator that allows us to solve the two-phase problem for any given geometry. The adjoint-based
<figure><img src="image_15.png" /><figcaption>Figure 15: Iterations 0, 15, 17 and 20 of the optimization procedure for the growing crystals test case. The blue curve represents the desired shape and the red one the final position of the interface at a given iteration. The final time is $t_f = 0.45$ and the interface is plotted with a time step of 0.045. The color map corresponds to the error in the temperature field. The inset shows the actuator $u$ at a given iteration (red) and the desired one (blue).</figcaption></figure>

optimization procedure is shown to be a robust algorithm to control the shape of a melting or solidification front, even in the presence of dendritic instabilities and anisotropic effects. As a future work, one could consider to solve the Navier-Stokes equations in the fluid phase, allowing for Rayleigh-Bénard instabilities in the case of the unidirectional solidification process.
<table><caption>Table 1: Optimization parameters and final results of the different considered cases. The columns J and ∇J correspond to the number of FP and AP calls respectively.</caption><thead><tr><th rowspan="2">Case</th><th colspan="4">Optimization parameters</th><th colspan="3">Results</th></tr><tr><th>β<sub>1</sub></th><th>β<sub>2</sub></th><th>β<sub>3</sub></th><th>β<sub>4</sub></th><th>J calls</th><th>∇J calls</th><th>J<sub>final</sub>/J<sub>0</sub></th></tr></thead><tbody><tr><td>Melting circle</td><td>1</td><td>10</td><td>0</td><td>10<sup>-3</sup></td><td>66</td><td>16</td><td>1.81 × 10<sup>-3</sup></td></tr><tr><td>Mullins-Sekerka</td><td>1</td><td>10</td><td>0.1</td><td>10<sup>-4</sup></td><td>54</td><td>10</td><td>1.89 × 10<sup>-3</sup></td></tr><tr><td>Growing crystals</td><td>1</td><td>1</td><td>0</td><td>10<sup>-2</sup></td><td>25</td><td>20</td><td>3.21 × 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>Figure 16: Evolution of the cost functional for the different cases.</figcaption></figure>

## Appendix A. Adjoint temperature derivation

We detail the derivations from the Lagrange function $\mathcal{L}$ (Eq. 5.2) to the adjoint Stefan problem (AP). We have the initial conditions,

$$ 
\begin{aligned}
T(x, 0) &= T_0(x), \\
\phi(x, 0) &= \phi_0(x)
\end{aligned}
\qquad (A.1) 
$$

meaning that when calculating the derivatives in the direction $\delta T$ and $\delta \phi$, we require: $\delta T(x, 0) = \delta \phi(x, 0) = 0$. We consider the adjoint temperature by setting $\mathcal{L}_T(\cdot) = 0$,

$$ 
\begin{aligned}
\mathcal{L}_T \delta T = & - \int_{\Omega} \beta_1 (T(t_f) - T_{t_f}) \delta T(t_f) \\
& - \int_0^{t_f} \int_{\Omega_1} \left( \frac{\partial \delta T_1}{\partial t} - \Delta \delta T_1 \right) \theta_1 - \int_0^{t_f} \int_{\Omega_2} \left( \frac{\partial \delta T_2}{\partial t} - \Delta \delta T_2 \right) \theta_2 \\
& - \int_0^{t_f} \int_{\partial \Omega} \delta T \theta_D - \int_0^{t_f} \int_{\Gamma} \delta T \theta_I \\
& - \int_0^{t_f} \int_{\Gamma} \psi [\nabla \delta T_i]_2^1 \cdot \nabla \phi
\end{aligned}
\qquad (A.2) 
$$

We move the spatial and temporal derivatives towards the adjoint state $\theta$. We apply integration by parts, once with respect to time (Corollary 1) and twice with respect to space (using Green's formula),
$$
\begin{equation}
\begin{aligned}
\mathcal{L}_T \delta T = & - \int_{\Omega} \beta_1 (T(t_f) - T_{t_f}) \delta T(t_f) \\
& - \int_{\Omega_1(t_f)} \delta T(t_f) \theta_1(t_f) + \int_{\Omega_1(0)} \delta T(0) \theta_1(0) \\
& - \int_0^{t_f} \int_{\Omega_1} \frac{\partial \theta_1}{\partial t} \delta T + \int_0^{t_f} \int_{\partial \Omega_1} \delta T \theta_1 V_1 \\
& + \int_0^{t_f} \int_{\partial \Omega_1} \frac{\partial \delta T}{\partial n} \theta_1 - \int_0^{t_f} \int_{\Omega_1} \nabla \delta T \nabla \theta_1 \\
& - \int_{\Omega_2(t_f)} \delta T(t_f) \theta_2(t_f) + \int_{\Omega_2(0)} \delta T(0) \theta_2(0) \\
& - \int_0^{t_f} \int_{\Omega_2} \frac{\partial \theta_2}{\partial t} \delta T + \int_0^{t_f} \int_{\partial \Omega_2} \delta T \theta_2 V_2 \\
& + \int_0^{t_f} \int_{\partial \Omega_2} \frac{\partial \delta T}{\partial n} \theta_2 - \int_0^{t_f} \int_{\Omega_2} \nabla \delta T \nabla \theta_2 \\
& - \int_0^{t_f} \int_{\partial \Omega} \delta T \theta_D - \int_0^{t_f} \int_{\Gamma} \delta T \theta_I \\
& - \int_0^{t_f} \int_{\Gamma} \psi [\nabla \delta T]_2^1 \cdot \nabla \phi
\end{aligned}
\tag{A.3}
\end{equation}
$$

with $V_1$ and $V_2$ the velocities of the control volumes $\Omega_1$ and $\Omega_2$ respectively with $V_1 = \vec{V} \cdot n$ and $V_2 = -\vec{V} \cdot n$.
These two terms are only non-zero on $\Gamma$. Sorting the terms by their domain of integration and setting
$\delta T(0) = 0$, we obtain,

$$
\begin{equation}
\begin{aligned}
\mathcal{L}_T \delta T = 0 & = \int_0^{t_f} \int_{\Omega_1} \left( \frac{\partial \theta_1}{\partial t} + \Delta \theta_1 \right) \delta T \\
& \quad + \int_0^{t_f} \int_{\Omega_2} \left( \frac{\partial \theta_2}{\partial t} + \Delta \theta_2 \right) \delta T \\
& \quad - \int_{\Omega_1} (\theta_1(t_f) - \beta_1 (T(t_f) - T_{t_f})) \delta T(t_f) \\
& \quad - \int_{\Omega_2} (\theta_2(t_f) - \beta_1 (T(t_f) - T_{t_f})) \delta T(t_f) \\
& \quad + \int_0^{t_f} \int_{\partial \Omega} -\frac{\partial \theta}{\partial n} \delta T + \frac{\partial \delta T}{\partial n} (\theta - \theta_D) \\
& \quad + \int_0^{t_f} \int_{\Gamma} \left( \theta_1 \vec{V} \cdot n - \theta_2 \vec{V} \cdot n - \frac{\partial \theta_1}{\partial n} + \frac{\partial \theta_2}{\partial n} - \theta_I \right) \delta T \\
& \quad + \int_0^{t_f} \int_{\Gamma} (\theta - \psi |\nabla \phi|) [\nabla \delta T]_2^1 \cdot n
\end{aligned}
\tag{A.4}
\end{equation}
$$

with the second to last term simplifying to, $\theta_I = -[\nabla\theta]_2^1 \cdot n$ as $\theta_1 = \theta_2$ on $\Gamma$.

**Theorem 1 (Reynolds transport theorem).** The derivative of the quantity

$$
F(t) := \int_{\Omega(t)} f(x, t) dx
$$
is given by

$$
\begin{align*}
\frac{dF}{dt}(t) &= \int_{\Omega(t)} \frac{\partial f}{\partial t} + \operatorname{div}(fV)dx \\
&= \int_{\Omega(t)} \frac{df}{dt} + f \operatorname{div}(V)dx \\
&= \int_{\Omega(t)} \frac{\partial f}{\partial t}dx + \int_{\partial\Omega(t)} fV \cdot nds
\end{align*}
$$

where V is the velocity field in which the control volume Ω(t) moves.

Corollary 1 (Integration by Parts in Time in Moving Domains). For g = g(x, t) and h = h(x, t),

we have

$$
\int_{0}^{t_f} \int_{\Omega(t)} g h_t dx dt = \int_{\Omega(t_f)} g(x, t_f) h(x, t_f) dx - \int_{\Omega(0)} g(x, 0) h(x, 0) dx \\
- \int_{0}^{t_f} \int_{\Omega(t)} g_t h dx dt - \int_{0}^{t_f} \int_{\partial \Omega(t)} g h V \cdot nds dt
$$

References

[1] B. Šarler, Stefan's work on solid-liquid phase changes, Engineering Analysis with Boundary Elements 16 (2) (1995) 83-92.
doi:10.1016/0955-7997(95)00047-X.

[2] S. Osher, J. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations, J Comput. Phys. 79 (1) (1988) 12-49. doi:10.1016/0021-9991(88)90002-2.

[3] D. Juric, G. Tryggvason, A front-tracking method for dendritic solidification, J Comput. Phys. 123 (1996) 127-148.

[4] G. Segal, C. Vuik, F. Vermolen, A conserving discretization for the free boundary in a two-dimensional stefan problem, J Comput. Phys. 141 (1998) 1-21.

[5] A. Hassan, T. Sayadi, V. LeChenadec, H. Pitsch, A. Attili, Adjoint-based sensitivity analysis of steady char burn out, Com. Theo. Modeling 25 (1) (2021) 96-120. doi:10.1080/13647830.2020.1838614.

[6] E. Javierre, C. Vuik, F. Vermolen, S. van der Zwaag, A comparison of numerical models for one-dimensional Stefan problems, J Comp. Applied Mathematics 192 (2) (2006) 445-459. doi:10.1016/j.cam.2005.04.062.

[7] M. Rose, An enthalpy scheme for stefan problems in several dimensions, App. Numerical Math. 12 (1993) 229.

[8] K. Brattkus, D. Meiron, Numerical simulations of unsteady crystal growth, SIAM J. Appl. Math. 52 (1992) 1303.

[9] A. Limare, S. Popinet, C. Josserand, Z. Xue, A. Ghigo, A hybrid level-set / embedded boundary method applied to solidification-melt problems, arXiv (2022). doi:10.1063/5.482020.

[10] J. S. Langer, Instabilities and pattern formation in crystal growth, Rev. Mod. Phys. 52 (1) (1980) 1-28. doi:10.1103/RevModPhys.52.1.

[11] W. W. Mullins, R. F. Sekerka, Stability of a planar interface during solidification of a dilute binary alloy, Journal of Applied Physics 35 (2) (1964) 444-451, publisher: American Institute of Physics. doi:10.1063/1.1713333.
URL https://aip.scitation.org/doi/10.1063/1.1713333

[12] A. Woods, Melting and dissolving, J. Fluid Mech. 239 (1992) 429. doi:10.1017/S0022112092004476.

[13] A. L. Marsden, J. A. Feinstein, C. A. Taylor, A computational framework for derivative-free optimization of cardiovascular geometries, Computer Methods in Applied Mechanics and Engineering 197 (21-24) (2008) 1890-1905.

[14] S. Pierret, R. F. Coelho, H. Kato, Multidisciplinary and multiple operating points shape optimization of three-dimensional compressor blades, Structural and Multidisciplinary Optimization 33 (1) (2007) 61-70.

[15] M. B. Giles, N. A. Pierce, An Introduction to the Adjoint Approach to Design, Flow Turb. Combus. 65 (3) (2000) 393-415.
[16] O. Pironneau, On optimum design in fluid mechanics, J Fluid Mech. 64 (1) (1974) 97-110.

[17] A. Jameson, Aerodynamic design via control theory, J Scientific Computing 3 (3) (1988) 233-260.

[18] A. Jameson, L. Martinelli, N. A. Pierce, Optimum Aerodynamic Design Using the Navier-Stokes Equations, Theor. Comp. Fluid Dyn. 10 (1) (1998) 213-237.

[19] M. Juniper, Triggering in the horizontal Rijke tube: non-normality, transient growth and bypass transition, J. Fluid Mech. 667 (2010) 272-308.

[20] M. Lemke, J. Reiss, J. Sesterhenn, Adjoint-based analysis of thermoacoustic coupling, ICNAAM (2013) 2163-2166.

[21] S. Schmidt, C. Ilic, V. Schulz, N. R. Gauger, Three-dimensional large-scale aerodynamic shape optimization based on shape calculus, AIAA J 51 (11) (2013) 2615-2627.

[22] S. Rabin, C. Caulfield, R. Kerswell, Designing a more nonlinearly stable laminar flow via boundary manipulation, J Fluid Mech. 738 (2014) R1:1-12.

[23] D. Foures, C. Caulfield, P. Schmid, Optimal mixing in two-dimensional plane Poiseuille flow at finite Péclet number, J Fluid Mech. 748 (2014) 241-277.

[24] K. Duraisamy, J. Alonso, Adjoint-based techniques for uncertainty quantification in turbulent flows with combustion, 42nd AIAA Fluid Dynamics Conference and Exhibit (2012) 25-28.

[25] A. Fikl, V. Le Chenadec, T. Sayadi, Control and optimization of interfacial flows using adjoint-based techniques, Fluids 5 (3) (2020). doi:10.3390/fluids5030156

[26] K. Ou, A. Jameson, Unsteady adjoint method for the optimal control of advection and burger's equations using high-order spectral difference method, American Institute of Aeronautics and Astronautics (2011) 1-18 doi:10.2514/6.2011-24
URL https://arc.aiaa.org/doi/10.2514/6.2011-24

[27] K. Braman, T. Oliver, V. Raman, Adjoint-based sensitivity analysis of flames, Comb. Theo. Modelling 19 (1) (2015) 29-56.

[28] M. Lemke, L. Cai, J. Reiss, H. Pitsch, J. Sesterhehn, Adjoint-based sensitivity analysis of quantities of interest of complex combustion models, Comb. Theo. Modelling 23 (1) (2019) 180-196.

[29] K. Hoffmann, J. Sprekels, Real-time control of the free boundary in a two-phase stefan problem, Numerical Functional Analysis and Optimization 5 (1) (1982) 47-76.

[30] P. Knabner, Control of stefan problems by means of linear-quadratic defect minimization, Numerische Mathematik 46 (3) (1985) 429-442.

[31] S. Kang, N. Zabaras, Control of the freezing interface motion in two-dimensional solidification processes using the adjoint method, International J Numerical Methods in Engineering 38 (1995) 63-80.

[32] Z. Yang, The adjoint method for the inverse design of solidification processes with convection, PhD thesis, Cornell University (1997).

[33] M. Hinze, S. Ziegenbalg, Optimal control of the free boundary in a two-phase stefan problem, J Comput. Phys. 223 (2) (2007) 657-684.

[34] M. Bernauer, R. Herzog, Optimal Control of the Classical Two-Phase Stefan Problem in Level Set Formulation, SIAM J. Sci. Comput. 33 (1) (2011) 342-363. doi:10.1137/100783327

[35] J. S. Langer, Instabilities and pattern formation in crystal growth, Reviews of Modern Physics (1980).

[36] S. Chen, B. Merriman, S. Osher, P. Smereka, A Simple Level Set Method for Solving Stefan Problems, J Comp. Phys. 135 (1) (1997) 8-29. doi:10.1006/jcph.1997.5721

[37] D. Juric, G. Tryggvason, A Front-Tracking Method for Dendritic Solidification, J Comp. Phys. 123 (1) (1996) 127-148. doi:10.1006/jcph.1996.0011

[38] T. Fullana, A. Q. Rodríguez, V. L. Chenadec, T. Sayadi, A cut cell method for the solution of the incompressible navier-stokes equations around stationary geometries, Journal of Computational Physics (2022).

[39] H. Johansen, P. Colella, A cartesian grid embedded boundary method for poisson's equation on irregular domains, J Comp.
Phys. 147 (1) (1998) 60-85. doi:10.1006/jcph.1998.5965.

[40] D. Peng, B. Merriman, S. Osher, H. Zhao, M. Kang, A PDE-Based Fast Local Level Set Method, J Comp. Phys. 155 (2) (1999) 410-438. doi:10.1006/jcph.1999.6345.

[41] K. Mikula, M. Ohlberger, J. Urbán, Inflow-implicit/outflow-explicit finite volume methods for solving advection equations, Applied Numerical Mathematics 85 (2014) 16-37. doi:10.1016/j.apnum.2014.06.002.

[42] K. Mikula, M. Ohlberger, A New Level Set Method for Motion in Normal Direction Based on a Semi-Implicit Forward-Backward Diffusion Approach, SIAM J. Sci. Comput. 32 (3) (2010) 1527-1544. doi:10.1137/09075946X.

[43] C. Min, On reinitializing level set functions, J Comp. Phys. 229 (8) (2010) 2764-2772. doi:10.1016/j.jcp.2009.12.032.

[44] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, Journal of Computational Physics 228 (16) (2009) 5838-5866. doi:10.1016/j.jcp.2009.04.042.
URL https://linkinghub.elsevier.com/retrieve/pii/S002199910900240X

[45] F. C. Frank, Radially symmetric phase growth controlled by diffusion, Proc. R. Soc. Lond (1950). doi:http://doi.org/10.1098/rspa.1950.0080.

[46] R. Almgren, Variational algorithms and pattern formation in dendritic solidification, J Comp. Phys. 106 (2) (1993) 337-354. doi:https://doi.org/10.1016/S0021-9991(83)71112-5.

[47] L. Tan, N. Zabaras, A level set simulation of dendritic solidification with combined features of front-tracking and fixed-domain methods, Journal of Computational Physics 211 (1) (2006) 36-63. doi:10.1016/j.jcp.2005.05.013.

[48] D. C. Liu., J. Nocedal, On the limited memory method for large scale optimization, Mathematical Programming 3 (45) (1989) 503-528.

[49] S. J. Wright, J. Nocedal, Numerical optimization, Springer, 2006.

[50] J. Bezanson, A. Edelman, S. Karpinski, V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Review 59 (1) (2017) 65-98. doi:10.1137/141000671.

[51] P. K. Mogensen, A. N. Riseth, Optim: A mathematical optimization package for Julia, Journal of Open Source Software 3 (24) (2018) 615. doi:10.21105/joss.00615.