Open AccessArticle

Efficient Discretization of the Laplacian: Application to Moving Boundary Problems

Sebastian-Josue Castillo

^†

and

Ferenc Izsák

^*,†

Department of Applied Analysis and Computational Mathematics, Eötvös Loránd University, 1117 Budapest, Hungary

Author to whom correspondence should be addressed.

^†

These authors contributed equally to this work.

Algorithms 2024, 17(12), 541; https://doi.org/10.3390/a17120541

Submission received: 22 October 2024 / Revised: 19 November 2024 / Accepted: 27 November 2024 / Published: 29 November 2024

(This article belongs to the Special Issue Numerical Optimization in Honor of the 60th Birthday of Marko M. Mäkelä)

Download

Browse Figures

Figure 1
A uniform (left) and a general (right) quadrilateral grid for finite difference approximations. A uniform rectangular grid (left) related to the approximation (<a href="#FD4-algorithms-17-00541" class="html-disp-formula">4</a>) and a non-uniform quadrilateral grid (right) related to (<a href="#FD5-algorithms-17-00541" class="html-disp-formula">5</a>). "> Figure 2
Location of non-zero elements in the discretization matrix <math display="inline"><semantics> <msub> <mi>D</mi> <mrow> <mi>h</mi> <mo>,</mo> <mi>FD</mi> </mrow> </msub> </semantics></math> corresponding to (<a href="#FD4-algorithms-17-00541" class="html-disp-formula">4</a>) (left) and (<a href="#FD5-algorithms-17-00541" class="html-disp-formula">5</a>) (right), respectively. The total matrix size <math display="inline"><semantics> <mrow> <mn>36</mn> <mo>×</mo> <mn>36</mn> </mrow> </semantics></math> corresponds to the number of the interior grid points. The notation <tt>nz</tt> yields the number of non-zeros. "> Figure 3
Spatial discretization using structured grids in two cases. Number of uniformly distributed grid points: 20 (left) and 40 (right) in the horizontal direction. Number of grid points: 16 (left) and 33 (right) in the vertical direction. "> Figure 4
Motion of the segments <math display="inline"><semantics> <mrow> <msub> <mi mathvariant="bold">x</mi> <mrow> <mi>i</mi> <mo>−</mo> <mn>1</mn> </mrow> </msub> <msub> <mi mathvariant="bold">x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi mathvariant="bold">x</mi> <mi>i</mi> </msub> <msub> <mi mathvariant="bold">x</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>,</mo> </mrow> </semantics></math> and <math display="inline"><semantics> <mrow> <msub> <mi mathvariant="bold">x</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <msub> <mi mathvariant="bold">x</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>2</mn> </mrow> </msub> </mrow> </semantics></math> of <math display="inline"><semantics> <msup> <mo>Γ</mo> <mi>n</mi> </msup> </semantics></math> in step (ii) (a). <math display="inline"><semantics> <msub> <mi>d</mi> <mi>i</mi> </msub> </semantics></math> denotes the shift vectors obtained from the jump term in (<a href="#FD26-algorithms-17-00541" class="html-disp-formula">26</a>). Dashed lines show the shifting of the segments of <math display="inline"><semantics> <msup> <mo>Γ</mo> <mi>n</mi> </msup> </semantics></math>. Their average on the vertical lines gives the new position of the vertices of <math display="inline"><semantics> <msup> <mo>Γ</mo> <mrow> <mi>n</mi> <mo>+</mo> <mn>2</mn> </mrow> </msup> </semantics></math>. "> Figure 5
The steps of shifting the interface <math display="inline"><semantics> <msup> <mo>Γ</mo> <mi>n</mi> </msup> </semantics></math> to obtain the new one <math display="inline"><semantics> <msup> <mo>Γ</mo> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> </semantics></math>. "> Figure 6
Initial condition (left, with a heat map) and the corresponding initial grid (right) for the simulation of the problem in (<a href="#FD11-algorithms-17-00541" class="html-disp-formula">11</a>) and (<a href="#FD12-algorithms-17-00541" class="html-disp-formula">12</a>). Bottom region: <math display="inline"><semantics> <msub> <mo>Ω</mo> <mrow> <mn>1</mn> <mo>,</mo> <mi>t</mi> </mrow> </msub> </semantics></math>, depicted with red color (left) and red grid points (right); top region: <math display="inline"><semantics> <msub> <mo>Ω</mo> <mrow> <mn>0</mn> <mo>,</mo> <mi>t</mi> </mrow> </msub> </semantics></math>, depicted with blue color (left) and blue grid points (right). The common interface and the corresponding grid points are colored with green. The horizontal distribution of grid points is uniform, while the vertical distribution becomes more dense towards the common interface. "> Figure 7
The result of the simulation procedure for the Stefan problem after 7 time steps (left) and 17 time steps (right), respectively. The corresponding parameters are given in the list below. Again, blue color represents negative temperature, while red indicates positive temperature. ">

Review Reports Versions Notes

Abstract

An efficient approximation is developed for the Laplacian operator by merging the advances of finite difference and finite element approximations. This approach is applicable to a general quadrilateral grid. The optimal coefficients for the approximation are computed using a pointwise optimization process. In this process, an overdetermined system is solved in the least-square sense using weighted polynomial approximation. The proposed algorithm is a vectorized procedure, keeping the computational time at a low level. The performance of this method is demonstrated on a model problem involving the numerical solution of a Poisson problem. Its true potential is evident when applied to moving boundary problems, which typically require a dynamic grid for efficient simulation. Within the framework of the proposed algorithm, we can compute the spatial discretization on the new grid quickly. This procedure is tested in the Stefan problem. For this, we give the simulation algorithm in detail utilizing the quadrilateral grid geometry. The performance is again demonstrated in a series of numerical experiments.

Keywords:

boundary value problems; Laplacian; moving boundary problem; finite difference discretization; optimization; non-uniform mesh; Stefan problem

1. Introduction

The numerical solution of PDEs is one of the most important problems in computational mathematics. In real-life cases, these equations serve as a mathematical model of different phenomena. To simulate them, the corresponding equations have to be solved. For this, usually, we have to perform a numerical solution procedure.

From a modeling perspective, the primary challenge in these simulations lies in defining and incorporating realistic boundary conditions. At the same time, the basis of any numerical solution is an appropriate discretization. Typically, this is applied to both spatial and temporal variables, provided that we deal with time-dependent problems.

Our focus will be spatial discretization. The computational methods include mainly the finite element [1], finite difference (FD) [2], and spectral or Fourier-type methods [3]. Alternatively, in the case of simple linear differential operators, a family of mesh-free methods (e.g., the method of fundamental solutions [4,5]) can also be applied. However, for complex spatial domains, which is common in realistic scenarios, both Fourier-type and traditional finite difference discretizations often fall short.

Accordingly, several attempts were made to extend finite difference methods for non-uniform grids, see [6,7].

At the same time, in the engineering practice, for the numerical solution of PDEs, usually finite element-based discretizations are used. More precisely, one should call them Galerkin-type discretizations, which include the non-conforming or the discontinuous Galerkin-type discretizations. However, a significant drawback of this approach is its time-consuming nature. Constructing a finite element discretization involves generating a comprehensive data structure that includes a hierarchical organization of vertices, edges, and faces, along with their interrelations. Moreover, within such a framework, one should assemble all data into a large algebraic equation. In the case of (local) linear differential operators, we obtain a large linear algebraic system with a sparse matrix. Here, the full discretization (i.e., the construction of the above matrix) can take more time compared to the solution of the resulting system. This can extremely slow down the simulation procedure, especially when the computational domain changes dynamically over time, such as in moving boundary problems. Addressing these issues remains a significant challenge in computational mathematics.

A main source of the enhanced computing time is that the above finite element assembling procedure can hardly be vectorized. Accordingly, to speed up the assembling, several procedures were proposed, which are also called sequencing. In [8], the authors introduced a vectorized algorithm for this purpose, which was then developed further in [9] in another way to reduce computing time. Another conventional way to speed up finite element assembling is to perform parallel processing and computations on GPUs. Regarding this, several attempts were also made, see, e.g., [10,11]. Whether all of these results led indeed to a significant speed-up or not, they all need some extra effort. In each case, one has to build up an extra structure to support the algorithm. In concrete terms, usually, the authors employ a kind of graph coloring algorithm to partition elements into non-overlapping sets. In this way, the elemental matrices can be assembled into a global finite element matrix in parallel without conflicts. In such procedures, a suitable sparse format of the global matrix is needed [11].

The aim of the present work is to propose a vectorized algorithm for assembling the discretization matrix of the Laplace operator in general domains.

The efficient (accurate and fast) discretization of the Laplacian operator has great importance in a number of physical models. For the distribution of electric potential, for the viscosity term in fluid dynamics, and for the modeling of irrotational flows using velocity potentials (e.g., in the case of linear wave propagation), one always makes use of this operator. Moreover, in many source-free and steady-state models, the Laplacian becomes zero. Our algorithm optimizes finite difference methods for general geometric configurations around a given grid point, enabling its application in diverse domains using only point data. We will demonstrate the performance of the proposed algorithm first for the solution of Laplace equations. Note that in a preceding work, instead of pointwise optimization, in the case of a lower-order FD approximation, we made use of a neural network [12] to automatize the optimization. The algorithm is then integrated into the numerical solution of a Stefan problem, serving as a critical component. This classic moving boundary problem, first proposed in 1891 to model the melting of ice, has since been applied to modeling crystallization and simulating the dynamics of two-phase materials [13]. In all cases, the evolution of the interface between the different materials or phases is determined by a jump relation between the normal components of the gradient. This becomes the key source of the dynamics, and can be derived as a conservation principle [14,15]. A number of studies first investigated the corresponding one-dimensional case [16].

Regarding the realistic multidimensional model, a few simulation procedures were proposed, beginning in the 1980s [17], and developed later [18]. With advancements in computational capacity and renewed interest in modeling melting ice surfaces, several new simulation algorithms have been developed [15]. Our approach is related to these by completing them with an efficient discretization procedure.

The outline of our contribution is the following. After this short introduction, we set the mathematical problem to solve. Then, the details of our pointwise optimization algorithm are given. The performance is demonstrated and analyzed in a case where the analytic solution of a Poisson problem is known. This method is then incorporated as a building block into the numerical solution of a two-dimensional two-phase Stefan problem. Here, a detailed algorithm is also given for moving the mesh in each time step. On this new mesh, the spatial discretization of the Laplacian has to be computed. Finally, the full simulation process is shown in a model problem illustrating the evolution of the grids and the interface between the two phases. Additionally, we summarize the results in a short discussion.

2. Materials and Methods

The proposed algorithm will be performed on the two-dimensional bounded domain

Ω

to solve numerically the problem

\{\begin{matrix} Δ u (x) = f (x) x \in Ω \\ u (x) = g (x) x \in \partial Ω, \end{matrix}

(1)

where

f \in L_{2} (Ω)

and

g \in C (\partial Ω)

are given. A generic grid point in our finite difference method will be denoted by

x_{0}

In the following, vector quantities along with their components will be denotes with boldface characters.

2.1. Problem Statement and Mathematical Background

The Laplacian of a function

u \in C^{2} (Ω)

will be approximated in

x_{0}

based on the corresponding function values in its neighboring grid points. Our structured mesh points can be recognized as the vertices of a generic quadrilateral mesh. Using such a mesh, one discretize

\bar{Ω}

and a number of software tools are available to perform the approximation.

For a classical survey on these, we refer to [19].

Here, the neighbors

x_{1}, x_{2}, x_{3}

and

x_{4}

x_{0}

are determined by the edges of all quadrilaterals containing

x_{0}

. Similarly,

x_{5}, x_{6}, x_{7}

and

x_{8}

are determined by the diagonals containing

x_{0}

. For a visualization of these, see Figure 1. At the same time, we do not need the entire structure of such a mesh, including the relation of quadrilaterals, edges and vertices.

We note that within the framework of a classical finite element method, we have to perform the following steps.

(i)

In the domain, we have to define a mesh, numbering its elements, construct a whole mesh structure, numbering the elements (now the quadrangles), and generate the corresponding adjacency matrices:

(a): for each node, the number of each subdomain (element) containing this node;
(b): for each element, the list of its nodes;
(c): for each element, the list of its neighboring elements.

(ii)

By mapping a reference element to each of the physical ones, one has to evaluate the element stiffness matrices, computing Jacobians for this mapping.

(iii)

In the assembling procedure, applying a loop, for each node, we have to collect the above contributions from the adjacent elements (given by the matrix in (a)).

To explain our approach, we consider two different discretizations of

Δ

in the grid point

x_{0}

. If only the data at

x_{0}, x_{1}, x_{2}, x_{3}

and

x_{4}

are used, we employ the approximation

Δ u (x_{0}) \approx \sum_{j = 0}^{4} a_{j} u (x_{j}) .

(2)

Otherwise, if all the data in

x_{0}, x_{1}, \dots, x_{8}

are used, we change (2) accordingly to obtain

Δ u (x_{0}) \approx \sum_{j = 0}^{8} a_{j} u (x_{j}) .

(3)

Note that in the the presence of a uniform rectangular grid with the parameters

h_{x}

and

h_{y}

, as shown in Figure 1, left, the optimal approximation in (2) can be given as

Δ u (x_{0}) \approx - (\frac{2}{h_{x}^{2}} + \frac{2}{h_{y}^{2}}) u (x_{0}) + \frac{1}{h_{x}^{2}} u (x_{1}) + \frac{1}{h_{x}^{2}} u (x_{3}) + \frac{1}{h_{y}^{2}} u (x_{2}) + \frac{1}{h_{y}^{2}} u (x_{4}),

(4)

so that the error becomes

O (h_{x}^{2}) + O (h_{y}^{2})

, assuming that u is four times differentiable. By changing the geometry in Figure 1, left, in any way, this accuracy is lost.

In a similar way, using a square grid with

h = h_{x} = h_{y}

and all the nine-point data in (3), the optimal approximation in (3) becomes

\begin{matrix} Δ u (x_{0}) & \approx - \frac{20}{6 h^{2}} u (x_{0}) + \frac{4}{6 h^{2}} \cdot (u (x_{1}) + u (x_{2}) + u (x_{3}) + u (x_{4})) \\ + \frac{1}{6 h^{2}} \cdot (u (x_{5}) + u (x_{6}) + u (x_{7}) + u (x_{8})) . \end{matrix}

(5)

In this case, provided that

Δ u = 0

, the error becomes

O (h^{4})

. For further details and extensions, we refer to [20]. Again, for a general geometry, such as for the one shown on the right of Figure 1, this approximation order is no longer valid.

In any case, we collect the coefficients in (4) or (5) into a matrix

D_{h, FD}

, which is called the finite difference discretization matrix. Here, the entry

D_{h, FD} [j, k]

is the coefficient of

u (x_{k})

in the approximation of

Δ u (x_{j})

, either in (4) or in (5). In this way,

D_{h, FD}

becomes a sparse matrix with either at most 5 or 9 non-zero elements per row; see Figure 2. In the case of regular geometries and homogeneous boundary conditions,

D_{h, FD}

is a positive definite (or semidefinite) symmetric one.

To compare the above cases, in Figure 2, we also give the position of the non-zero elements and the total number of non-zeros (denoted with nz) in the matrices

D_{h, FD}

corresponding to the discretizations in (4) and (5), respectively.

The conventional way to find optimal parameters in (2) and (3) is to use Taylor expansion and determine the coefficients to obtain the desired approximation order. In a general case, this approach completely fails since then the accuracy in (4) and (5) are not available and for lower-order accuracy, multiple solutions exist. In this way, we somehow have to optimize the coefficients in (2) and (3).

In a previous study [12], instead of pointwise optimization procedures, a neural network was utilized, which “learns” this process and can quickly compute the corresponding coefficients. That study focused on the five-point approximation in (4). With the matrix

D_{h, FD}

at hand, the final form of the FD discretization of the problem in (1) is the linear system

D_{h, FD} u_{h, FD} = f_{h},

(6)

where

f_{h}

denotes the discretization of the source term f with

f_{h} [j] = f (x_{j})

. The solution

u_{h, FD}

is then called a finite difference approximation of u in (1).

2.2. The Pointwise Optimization

To achieve possibly the highest-order approximation in (2) and (3), these should hold with equality for some polynomials until a certain order.

For example, if we have a uniform mesh with the grid size h and the coefficients are scaled as

h^{- 2}

, then to obtain an approximation order

N - 1

, the equality

Δ p (x_{0}) = a_{0} p (x_{0}) + a_{1} p (x_{1}) + \dots + a_{8} p (x_{8}) .

(7)

should hold for all polynomials p up to order N. Henceforth, we focus to this approximation with nine grid points.

At the same time, for general mesh geometries, the approximation order is decreased and, accordingly, (7) is not valid any more for all polynomials up to order N. The principal problem here is to determine them, which should be taken into account in (3). Also, we demand that the right-hand side of (3) is zero for any constant function. This immediately gives the constraint

a_{0} = - (a_{1} + a_{2} + \dots + a_{8})

for the coefficients in (3). To satisfy (7), we will first consider the set of polynomials

\begin{matrix} p_{1} (x, y) = 100 x, p_{2} (x, y) = 100 y, p_{3} (x, y) = 10 x y, \\ p_{4} (x, y) = 10 x^{2}, p_{5} (x, y) = 10 (x^{2} - y^{2}), \end{matrix}

(8)

which span the space of first- and second-order ones. Additionally, we take

\begin{matrix} p_{6} (x, y) = x^{3} - 3 x^{2} y, p_{7} (x, y) = y^{3} - 3 x y^{2}, \\ p_{8} (x, y) = x^{3} y - x y^{3}, p_{9} (x, y) = x^{4} + 6 x^{2} y^{2} - y^{4}, \\ p_{10} (x, y) = 5 x^{4} y - 10 x^{2} y^{3} + y^{5}, p_{11} (x, y) = x^{5} - 10 x^{3} y^{2} + 5 x y^{4}, \end{matrix}

which are the harmonic polynomials up to order five. To motivate this choice, we recall that in a number of real-life cases, we have to compute with harmonic functions, such that

Δ u = 0

. In other words, we can only have a certain accuracy for approximating the solution of

Δ u = f

, if we can perform it in the case of

f = 0

. In this way, an approximation in (5) has the form

a_{1} p_{j} (x_{1}) + a_{2} p_{j} (x_{2}) + \dots + a_{8} p_{j} (x_{8}) \approx Δ p_{j} (x_{0}), j = 1, 2, \dots, 11 .

(9)

Performing it for all polynomials

p_{1}, p_{2}, \dots, p_{11}

above, for the geometrical setup given with

x_{0}, x_{1}, \dots, x_{8}

, we obtain the following optimization task:

Find the vector of coefficients $a = {[a_{1}, a_{2}, \dots, a_{8}]}^{T}$ in (9) such that $P \cdot a \approx Δ p$ , where
-
$P \in R^{11 \times 8}$ with $P_{j k} = p_{j} (x_{k})$ ;
-
$Δ p = {[0, 0, 0, 20, 0, 0, 0, 0, 0, 0, 0]}^{T} \in R^{11}$ , so that $Δ p_{j} = p_{j} (x_{0})$ .

Using the coefficient 100 in the first-order polynomials

p_{1}

and

p_{2}

and the coefficient 10 in the second-order ones

p_{3}, p_{4}

and

p_{5}

in (8), we force more accuracy in the case of the lower-order polynomials. This procedure significantly improves the approximation in (9).

The optimization task in (9) is solved in the least-square sense for each point

x_{0}

using the np.linalg.lstsq subroutine in Python. This can be performed for all grid points independently. In this way, we can easily vectorize the procedure in the following way:

We create a matrix, where each row contains the geometry around a fixed grid point.
The optimization procedure is applied (independently) to the rows of this matrix using the np.apply_along_axis subroutine in Python.

Accordingly, this final step could also be executed with parallel processing.

2.3. Numerical Experiments

We have first tested the performance of the above optimization task on a fixed domain.

2.3.1. The Model Problem

We investigate the numerical solution of the test problem

\{\begin{matrix} Δ u (x, y) = - 16 x cos 4 y (x, y) \in Ω \\ u (x, y) = u_{g} (x, y) (x, y) \in \partial Ω, \end{matrix}

(10)

where the computational domain is given with

Ω = {(x, y) \in R^{2} : 0 < x < \frac{π}{2}, 0 < y < 1 + \frac{cos x}{4}}

such that

u_{g} (0, y) = y

on the left,

u_{g} (\frac{π}{2}, y) = y + \frac{π}{2} \cdot cos 4 y

on the right,

u_{g} (x, 0) = x

on the bottom and

u (x, 1 + \frac{cos x}{4}) = 1 + \frac{cos x}{4} + x \cdot cos 4 (1 + \frac{cos x}{4})

on the top of

\partial Ω

We will compare the numerical solution with the analytic one

u (x, y) = y + x \cdot cos 4 y

2.3.2. Spatial Discretization

We utilize a scale of non-uniform structured grids in

Ω

. In each case, the vertical positions of the grid points will be uniformly distributed, while in the horizontal direction, we fix only the number of grid points. Here, the grid spacing is refined in the vicinity of the top boundary. This is a usual situation in the case of modeling water waves, as the dynamics are primary influenced by the motion of the free surface. Two different grids are shown in Figure 3.

2.3.3. Numerical Results

To validate our method for approximating the Laplacian, we evaluate the accuracy of the numerical solution of (10). Here, we measure the accuracy in the discrete

L_{2}

-norm depending on the number of grid points. The results, displayed in Table 1, clearly demonstrate the efficiency of the proposed vectorized discretization method. For a relatively fine grid with 2420 grid points, the construction of the FD discretization matrix

D_{h, FD}

took approximately one-sixth of a second.

2.4. Application to the Numerical Solution of the Stefan Problem

We introduce the Stefan problem and specify the boundary and initial conditions we apply. Then, the simulation algorithm is given in details. Finally, we present and discuss the results of the numerical experiments.

2.4.1. General Setting

Stefan problems have been used as a mathematical model of different physical phenomena. In a number of applications, this describes the interaction of different phases or materials. This can be a dissolution of a chemical compound into an alloy [18], the solidification of a binary alloy [21], a chrystallization process [22] or just the melting and freezing of ice [23].

Taking a minimal model, the unknown function

u : Ω \to R

in the Stefan problem corresponds to the temperature in a melting system [24].

At each time

t \in [0, T]

, the domain

Ω

is separated into the two disjoint ones

Ω_{0, t} = {x \in Ω : u (t, x) < 0} and Ω_{1, t} = {x \in Ω : u (t, x) > 0},

so that

\bar{Ω} = \bar{Ω_{0, t}} \cup \bar{Ω_{1, t}}

. Here,

Ω_{0, t}

represents the domain filled with ice, while

Ω_{1, t}

represents the domain filled with water. Accordingly, the set

Γ_{t} = {x = (x (t), y (t)) \in \bar{Ω} : u (t, x) = 0}

gives their common boundary.

We also introduce the notation

u_{0} = u {|_{[0, T] \times Ω_{0, t}} and u_{1} = u |}_{[0, T] \times Ω_{1, t}},

for the solutions on the subdomains above, where, for the sake of simplicity, we omit the time parameter t. With these, the governing equations can be given as

\{\begin{matrix} \partial_{t} u_{0} (t, x) = D_{0} \cdot Δ u_{0} (t, x), x \in Ω_{0, t} \\ \partial_{t} u_{1} (t, x) = D_{1} \cdot Δ u_{1} (t, x), x \in Ω_{1, t} . \end{matrix}

(11)

Here,

D_{0}

and

D_{1}

are the diffusion coefficients in the corresponding domains, where usually

D_{0} > D_{1}

. The evolution of the common boundary is given by

\partial_{t} Γ_{t} (x) = - L_{H} \cdot [[\partial_{ν} u (t, x (t), y (t))]],

(12)

where

Γ_{t}

is considered as a level set parametrized with t and

L_{H}

denotes the latent heat of solidification [25]. Physically, this gives the speed of the interface between the ice and liquid water. The symbol

[[\cdot]]

is for the jump operator, which is defined on the common boundary

Γ_{t}

with

[[\partial_{ν} u (t, x)]] = D_{0} \cdot ν_{0} (x) \cdot \nabla u_{0} (t, x) + D_{1} \cdot ν_{1} (x) \cdot \nabla u_{1} (t, x),

(13)

where

ν_{0}

and

ν_{1}

denote the corresponding (opposite) outward normals.

Note again that depending on the model, the unknown could be the concentration of some chemical compound or the degree of crystallinity. Also, in a more involved model, multiple phases (even for ice [23]) or a separate nucleation process [22] may occur.

2.4.2. Specific Initial and Boundary Conditions

The simulation will be initiated using the initial condition

u (0, x, y) = 1 + \frac{cos x}{4} - y .

During the time steps, the top boundary condition is set to be

u (t, x, 2 + \frac{cos x}{2}) = - 1 + (\frac{t}{4} - 1) \frac{cos x}{4}, x \in [- 1, 1],

(14)

such that, initially, it is given by

u (0, x, 2 + \frac{cos x}{2}) = - 1 - \frac{cos x}{4} .

We also set up the bottom boundary condition as

u (t, x, 0) = 1 + (t + 1) \cdot \frac{cos x}{4}, x \in [- 1, 1] .

(15)

which, initially, takes the value

u (0, x, 0) = 1 + \frac{cos x}{4} .

The right and left boundary conditions are homogeneous Neumann-type, as

\partial_{x} u (0, 1, y) = 0 and \partial_{x} u (0, - 1, y) = 0 y \in (0, 2 + \frac{cos (1)}{4}) .

Note that beyond to this specific case, the algorithm will be given in case of general Neumann-type boundary data.

2.4.3. Construction of the Grid

In order to track the evolution of the computational domains

Ω_{0}

and

Ω_{1}

, we apply a structured quadrilateral grid. Since the mesh is moving only vertically, we fix the x coordinates

x_{1}, x_{2}, \dots, x_{N x}

of the grid points. Their number—in the interior—is denoted with

N x

such that the horizontal grid spacing is

h x = 2 / (N x + 1)

. To incorporate the Neumann boundary conditions on the left and right, we also include a ghost grid point near to both the left and the right boundaries.

We also fix the number of interior grid points in the vertical direction; this is

N y / 2

both in

Ω_{0}

and

Ω_{1}

. Additionally, we place grid points on the fixed top and bottom boundaries as well as to the moving interface. According to (2) and (3), we apply a row-wise numbering for the grid points. To visualize the notations above, the grid point structure and the numbering are shown in Figure 4.

During the simulations, the interface evolves, requiring adjustments to the vertical positions of the grid points. For a fixed horizontal position, we choose a uniform vertical mesh. In this way, regarding the computing time, generating a new grid is a rather cheap operation. At the same time, the vertical distribution of the grid points is non-uniform. In the vicinity of the interface

Γ_{t}

, we apply a finer grid. This refinement allows for a more accurate approximation of the evolution of

Γ_{t}

. In concrete terms, for computing the jump terms (12), a fine grid will make our interpolation procedure more accurate. Note that a similar procedure can be used if free surface water waves are simulated [26].

2.5. Principles of the Algorithm for Solving (11)–(13)

The construction of the main building blocks of our algorithm follows the classical setting to solve moving boundary problems. Accordingly, the sketch of our algorithm is the following, where the main steps are repeated to simulate the evolution of the heat distribution.

(i)

We have a diffusion problem both on the top domain

Ω_{0}

and on the bottom domain

Ω_{1}

with some initial data.

(a): We construct the discretization of the Laplacian using the pointwise optimization algorithm in Section 2.2.
(b): We incorporate the Dirichlet- and Neumann-type boundary conditions in Section 2.4.2.
(c): We perform a time step using discrete problems on the two subdomains.

(ii)

We have the new Neumann boundary data on the moving surface.

(a): We move the surface according to Equation (13).
(b): We interpolate the solutions on the two new domains.

2.5.1. STEP (i) (b)

We now discuss step (i) (b). According to the implicit Euler time discretization, we have the general scheme

u^{n + 1} = u^{n} + δ \cdot (Δ_{h} u^{n + 1} + f^{n + 1}),

(16)

where the matrix

Δ_{h}

contains data from the virtual outward points. To incorporate them into the scheme, we will use the given Neumann boundary data. To simplify the notations, the algorithm will be presented on the left side of the computational domain.

Here, a generic interior node has the index K and its left virtual neighbors have the indices

K - N x - 3, K - 1

and

K + N x + 1

. Yielding the components of matrices and vectors in the subscript, and using the notation

\begin{matrix} J = { & K - N x - 3, K - N x - 2, K - N x - 1, K - 1, K, \\ K + 1, K + N x + 1, K + N x + 2, K + N x + 3}, \end{matrix}

we can express the component K of the product

Δ_{h} u^{n + 1}

{(Δ_{h} u^{n + 1})}_{K} = \sum_{j \in J} {(Δ_{h})}_{K, j} u_{j}^{n + 1} .

(17)

Here, we have to eliminate the sum

\begin{matrix} {(Δ_{h})}_{K, K - N x - 3} u_{K - N x - 3}^{n + 1} + {(Δ_{h})}_{K, K - 1} u_{K - 1}^{n + 1} + {(Δ_{h})}_{K, K + N x + 1} u_{K + N x + 1}^{n + 1} . \end{matrix}

(18)

Using the second-order approximations

\begin{matrix} u_{K - N x - 3}^{n + 1} - u_{K - N x - 2}^{n + 1} \approx \partial_{ν} u ((n + 1) δ, 0, y_{K - 1}) \\ u_{K - 1}^{n + 1} - u_{K}^{n + 1} \approx \partial_{ν} u ((n + 1) δ, 0, y_{K}) \\ u_{K + N x + 1}^{n + 1} - u_{K + N x + 2}^{n + 1} \approx \partial_{ν} u ((n + 1) δ, 0, y_{K + 1}), \end{matrix}

we can rewrite (17) as

\begin{matrix} {(Δ_{h} u^{n + 1})}_{K} = {(Δ_{h})}_{K, K - N x - 3} (u_{K - N x - 3}^{n + 1} - u_{K - N x - 2}^{n + 1}) + {(Δ_{h})}_{K, K - 1} (u_{K - 1}^{n + 1} - u_{K}^{n + 1}) \\ + {(Δ_{h})}_{K, K + N x + 1} (u_{K + N x + 1}^{n + 1} - u_{K + N x + 2}^{n + 1}) \\ + ({(Δ_{h})}_{K, K - N x - 3} + {(Δ_{h})}_{K, K - N x - 2}) u_{K - N x - 2}^{n + 1} + ({(Δ_{h})}_{K, K - 1} + {(Δ_{h})}_{K, K}) u_{K}^{n + 1} \\ + ({(Δ_{h})}_{K, K + N x + 2} + {(Δ_{h})}_{K, K + N x + 1}) u_{K + N x + 2}^{n + 1} \\ + {(Δ_{h})}_{K, K - N x - 1} u_{K - N x - 1}^{n + 1} + {(Δ_{h})}_{K, K + 1} u_{K + 1}^{n + 1} + {(Δ_{h})}_{K, K + N x + 3} u_{K + N x + 3}^{n + 1} \\ \approx {(Δ_{h})}_{K, K - N x - 3} \partial_{ν} u ((n + 1) δ, - 1, y_{K - 1}) + {(Δ_{h})}_{K, K - 1} \partial_{ν} u ((n + 1) δ, - 1, y_{K}) \\ + {(Δ_{h})}_{K, K + N x + 1} \partial_{ν} u ((n + 1) δ, - 1, y_{K + 1}) \\ + ({(Δ_{h})}_{K, K - N x - 3} + {(Δ_{h})}_{K, K - N x - 2}) u_{K - N x - 2}^{n + 1} + ({(Δ_{h})}_{K, K - 1} + {(Δ_{h})}_{K, K}) u_{K}^{n + 1} \\ + ({(Δ_{h})}_{K, K + N x + 2} + {(Δ_{h})}_{K, K + N x + 1}) u_{K + N x + 2}^{n + 1} \\ + {(Δ_{h})}_{K, K - N x - 1} u_{K - N x - 1}^{n + 1} + {(Δ_{h})}_{K, K + 1} u_{K + 1}^{n + 1} + {(Δ_{h})}_{K, K + N x + 3} u_{K + N x + 3}^{n + 1} \end{matrix}

(19)

such that we can eliminate the sum in (18). Accordingly, we perform the following:

on the right-hand side of (16), we have to add $δ$ times the first three terms in (19) modifying $f^{n + 1}$ ;
we have to modify the matrix elements ${(Δ_{h})}_{K, K - N x - 2}, {(Δ_{h})}_{K, K}$ and ${(Δ_{h})}_{K, K - N x + 2}$ by adding the entries ${(Δ_{h})}_{K, K - N x - 3}, {(Δ_{h})}_{K, K - 1}$ and ${(Δ_{h})}_{K, K - N x + 1}$ , respectively, corresponding to the virtual grid points as shown in the next three terms in (19).

Additionally, we have to incorporate the Dirichlet boundary conditions. Regarding the common boundary, the situation is easy: we do not have to change anything except for leaving the virtual points located here. On the top of the top domain, we have to collect the Dirichlet data

u_{K + N x + 1}^{n + 1}, u_{K + N x + 2}^{n + 1}, u_{K + N x + 3}^{n + 1}

associated with the interior point with index K. Then, on the right-hand side of (16), we have to keep

δ \cdot ({(Δ_{h})}_{K, K + N x + 1} u_{K + N x + 1}^{n + 1} + {(Δ_{h})}_{K, K + N x + 2} u_{K + N x + 2}^{n + 1} + {(Δ_{h})}_{K, K + N x + 3} u_{K + N x + 3}^{n + 1})

modifying with this the source term

f^{n + 1}

again, while the rest has to be moved to the left-hand side. The new source term is denoted with

{\tilde{f}}^{n + 1}

Performing a similar process on the bottom face of the bottom domain, we finally obtain the linear system

(I - δ \cdot ({\tilde{Δ}}_{h})) u^{n + 1} = u^{n} + δ \cdot {\tilde{f}}^{n + 1},

(20)

where

${\tilde{Δ}}_{h}$ contains only entries corresponding to interior elements, which are modified according to the Neumann boundary data in (17).
Similarly, the new source term, according to the incorporation of Neumann and Dirichlet data, is denoted with ${\tilde{f}}^{n + 1}$ .

2.5.2. STEP (i) (c)

For performing a time step, we only have to compute the numerical solution of the system in (20). Here, we simply use the built-in linear solver of Python. It is crucial that we apply implicit time stepping. In this way, we can maintain the stability of the method.

2.5.3. STEP (ii) (a)

We denote by

Γ^{n}

the approximation of

Γ_{n δ}

. As we have a quadrilateral mesh,

Γ^{n}

becomes a broken line. It is determined by a number of

N x

points such that their horizontal coordinates are fixed and uniformly distributed.

According to (13), we can approximate the evolution of

Γ_{t}

at each point

(x (t), y (t)) \in Γ_{t}

using the first-order approximation of (12) as

\begin{matrix} (x ((n + 1) δ), y ((n + 1) δ)) \approx (x (n δ), y (n δ)) + δ \cdot \partial_{t} Γ_{n δ} (x (n δ), y (n δ)) \\ = (x (n δ), y (n δ)) + δ \cdot [[\partial_{ν} u (x (n δ), y (n δ))]], \end{matrix}

(21)

which also defines a time stepping. Observe that using only (21) for the given grid points results in general points on

Γ^{n + 1}

, which do not necessarily have x-coordinates fitting to the mesh, Figure 4.

Before dealing with this problem, we should first compute

[[\partial_{ν} u^{n}]]

, for which we need the gradients

\nabla u_{i}^{n}

. With these, for the jump term in (21), we have

[[\partial_{ν} u (x (n δ), y (n δ))]] = (ν \cdot (\nabla u_{0}^{n} - \nabla u_{1}^{n})) ν .

Another chief problem is that function values are given in the grid points, where the outward normal of

Γ^{n}

is not defined.

Instead, to move each broken line segment of

Γ^{n}

, we compute the normal vector

ν

in the midpoint of these sections. These are given as

\frac{x_{i} + x_{i + 1}}{2}, i = 1, 2, \dots, N x

, see Figure 4. We have to approximate the gradients here, which is performed in the following way.

We first develop a bilinear interpolation of the numerical solution on the quadrilateral with the vertices

x_{i}, x_{i + 1}, x_{i + N x + 3}, x_{(i + 1) + N x + 3}

lying on the top of the zero level set as follows.

u_{0}^{n} (x, y) \approx a_{0} + a_{1} x + a_{2} y + a_{3} x y,

(22)

where the coefficients

a_{j}

are obtained by solving a linear system. Shifting the origin to

x_{i} = (x_{j}, y_{k})

and evaluating (22) here, we first observe that

a_{0} = 0

. Similarly, relating point values and the right-hand side of (22) at the points

x_{i}, x_{i + 1}, x_{i + N x + 3}, x_{(i + 1) + N x + 3}

, we obtain the equalities

\{\begin{matrix} 0 = a_{1} h_{x} + a_{2} (y_{k + 1} - y_{k}) + a_{3} h_{x} (y_{k + 1} - y_{k}) \\ u_{0}^{n + 1} (x_{j + 1 + N x + 3}) = a_{1} h_{x} + a_{2} (y_{k + 1 + N x + 3} - y_{k}) + a_{3} h_{x} (y_{k + 1 + N x + 3} - y_{k}) \\ u_{0}^{n} (x_{j + N x + 3}) = a_{2} (y_{k + N x + 3} - y_{k}), \end{matrix}

where the number in the bracket of the right-hand side yields the vector component. Similarly, we interpolate the numerical solution on the quadrilaterals lying on the bottom of

Γ^{n}

to obtain

u_{1}^{n} (x, y) \approx b_{0} + b_{1} x + b_{2} y + b_{3} x y .

(23)

Corresponding to the previous case, we obtain

\{\begin{matrix} u_{1}^{n} (x_{i - (N x + 3)}) = b_{2} (y_{k - (N x + 3)} - y_{k}) \\ u_{1}^{n} (x_{i + 1 - (N x + 3)}) = b_{1} h_{x} + b_{2} (y_{k + 1 - (N x + 3)} - y_{k}) + b_{3} h_{x} (y_{k + 1 - (N x + 3)} - y_{k}) \\ 0 = b_{1} h_{x} + b_{2} (y_{k + 1} - y_{k}) . \end{matrix}

Then, by (22) and (23), we obtain

\nabla u_{0}^{n} (x_{j + 1 / 2}, y_{k + 1 / 2}) = (a_{1} + a_{3} y_{k + 1 / 2}, a_{2} + a_{3} x_{j + 1 / 2}),

(24)

\nabla u_{1}^{n} (x_{j + 1 / 2}, y_{k + 1 / 2}) = (b_{1} + b_{3} y_{k + 1 / 2}, b_{2} + b_{3} x_{j + 1 / 2}) .

(25)

Next, we move the segment according to (21). That is, we compute

d_{i} = δ \cdot (ν_{i} \cdot (\nabla u_{0}^{n} - \nabla u_{1}^{n})) ν_{i}

(26)

using (24), (25) and shift the segment between

x_{i}

and

x_{i + 1}

with

d_{i}

in (26).

In this way, applying it on two neighboring segments between

x_{i}

and

x_{i + 1}

and between

x_{i}

and

x_{i + 1}

, we have two new vertices of

Γ^{n + 1}

on the fixed vertical lines. Therefore, instead, we take their average, which gives the location of the interface

Γ^{n + 1}

on the corresponding vertical line, as shown in Figure 4. Also, the main steps of the procedure in step (ii) (a) are shown in the block diagram in Figure 5.

2.5.4. STEP (ii) (b)

Finally, we have the new grid points of

Γ^{n + 1}

and the fixed ones on the bottom and the top of

Ω

. Between these, we again generate the new grid points according to Section 2.4.3 such that the consecutive vertical distances follow a geometric series.

On the top domain

Ω_{0, t}

, on a fixed vertical line, we then have the new grid point

x_{i}^{new}

of the interface

Ω^{n + 1}

(see Figure 4) and the generated new ones

x_{i + N x + 3}^{new} \in (x_{i}^{new}, x_{i + N x + 3}), x_{i + 2 N x + 6}^{new} \in (x_{i + N x + 3}, x_{i + 2 N x + 6}), \dots .

Also, we have the approximations

\begin{matrix} u ((n + 1) δ, x_{i}^{new}) = 0, u ((n + 1) δ, x_{i + N x + 3}) = u_{i + N x + 3}^{n + 1}, \\ u ((n + 1) δ, x_{i + 2 N x + 6}) = u_{i + 2 N x + 6}^{n + 1}, \dots . \end{matrix}

(27)

Accordingly, in the new points

x_{i + N x + 3}^{new}, x_{i + 2 N x + 6}^{new}, \dots

, we use the linear interpolation between the values in (27), replacing the values

u_{i + N x + 3}^{n + 1}, u_{i + 2 N x + 6}^{n + 1}, \dots

, which constitutes the starting point of our new time step from

(n + 1) δ

(n + 1) δ

2.5.5. Numerical Results

Our model problem represents the physical scenario of a melting volume with heat applied from the bottom. This is modeled using positive Dirichlet boundary conditions on the fixed bottom of

Ω_{1, t}

. Additionally, negative Dirichlet boundary conditions are applied to the top surface of

Ω_{0, t}

to account for the outward heat flow due to a negative temperature gradient.

The initial condition and the corresponding initial grid points are shown in Figure 6. The regions and grid point corresponding to positive temperature are depicted with red and those with negative temperature with blue color. The discretization of the interface

Γ_{t}

is defined by the green grid points.

In practical scenarios involving the modeling of ice layer melting, the characteristic horizontal dimension is often significantly larger than their average thickness. In this way, their interface becomes flat and horizontal [23], corresponding to the numerical experiments in this manuscript.

Elapsing the time, the material in

Ω

is melting such that the zero temperature level is moving upwards. This is also shown in the numerical simulations in Figure 7. Since we applied homogeneous Neumann boundary conditions, starting from the left and right boundary, the front will be flattened. Also, by the warm-up both from the bottom and the top boundary, corresponding to the boundary conditions in (14) and (15), the “cold” region shrinks significantly.

The simulation data are summarized as follows.

Governing equations and the corresponding initial and boundary data are given in Section 2.4.1 and Section 2.4.2, respectively.
Diffusion coefficients:
-
$D_{0} = 2.02 \cdot 10^{- 1}$ on the top domain;
-
$D_{1} = 0.143 \cdot 10^{- 2}$ on the bottom domain.
Number of interior grid points:
-
$N x = 41$ —horizontal;
-
$N y = 48$ —vertical: 24 + 24 interior points taking twice the common interface.
Time step: $δ = 0.08$ .
Number of time steps: 30.
Simulation time of one step: 1.05 s.
Software and hardware: Python code on a Laptop with Intel i7-6700 processor with 2.6 GHz.

3. Discussion and Conclusions

An optimization-based FD discretization of the two-dimensional Laplace operator was developed. The core of the optimization procedure was the least-square solution of an overdetermined system, where a special weighting of the corresponding equations was utilized. The pointwise execution of this procedure provides a vectorized alternative to the conventional assembly procedure.

The use and benefits of this algorithm are demonstrated on the Stefan problem, which is a mathematical model of the dynamics of two-phase materials with a moving interface. With this, the computational domains also evolve with time. Consequently, at each time step of the numerical simulations, a new spatial discretization must be performed, making the speed of this process critical to the overall efficiency of the numerical method. The proposed algorithm is notably fast, relying solely on neighboring relationships between grid points. Furthermore, it is effective on non-uniform and non-rectangular grids, combining the strengths of traditional finite difference and finite element methods.

The primary challenge in such simulations is that the interface often becomes unstable. This instability can lead to spurious oscillations on the free surface,

Γ_{t}

, which may grow over time and ultimately disrupt the simulation. To mitigate this issue, several stabilization techniques have been proposed, see [18]. To validate our approach in this context, the choice of a flat, common interface proves to be highly suitable.

Also, the algorithm proposed here is rather quick and uses only neighboring relations between the grid points. At the same time, it also works on non-uniform and non-rectangular grids, merging the advances of classical finite difference and finite element methods. Note that using classical finite element on such grids implies only first-order convergence [27].

The main application field of this new algorithm is the family of problems with moving boundaries or interfaces, including the more complex Stefan problems mentioned at the beginning of Section 2.4.1. Another very important family of such problems arises by modeling free-surface water waves. As the waves evolve, the domain filled with water will change step-by-step. When forecasting tsunamis [28] in real time, it is highly important to reduce the computing time to an affordable level for these 3-dimensional problems.

The disadvantage of the method is that the optimization procedure highly depends on the weighting, which was applied for the harmonic polynomials to obtain the optimal coefficients. Another related problem is that in the proposed algorithm, we cannot guarantee a fixed spatial convergence (or consistency) order. For irregular geometries, even the optimal coefficients will not deliver a very accurate approximation of the differential operator.

A natural extension of this work is to apply the method for approximating the surface Laplacian (or Laplace–Beltrami operator), a significant topic in current research. Also, the principle in Section 2.2 is applicable for any non-linear differential operator D. For this, we have to modify the overdetermined system in the optimization step, where on the left-hand side of (7), we have to apply D.

Author Contributions

Conceptualization, F.I.; methodology, F.I.; software, S.-J.C.; validation, S.-J.C.; formal analysis, F.I.; investigation, S.-J.C.; resources, S.-J.C.; data curation, S.-J.C.; writing—original draft preparation, F.I.; writing—review and editing, S.-J.C.; visualization, S.-J.C.; supervision, F.I.; project administration, F.I.; funding acquisition, F.I. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Research, Development and Innovation Office within the framework of the Thematic Excellence Program 2021—National Research Sub programme: “Artificial intelligence, large networks, data security: mathematical foundation and applications”.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

FD	Finite difference
PDE	Partial differential equation

References

Ern, A.; Guermond, J.-L. Theory and Practice of Finite Elements; Springer Science+Business Media: New York, NY, USA, 2004. [Google Scholar]
Thomas, J.W. Numerical Partial Differential Equations: Finite Difference Methods; Springer: New York, NY, USA, 1995. [Google Scholar]
Hesthaven, J.S.; Gottlieb, S.; Gottlieb, D. Spectral Methods for Time-Dependent Problems; Cambridge University Press: Cambridge, UK, 2007. [Google Scholar]
Fairweather, G.; Karageorghis, A. The method of fundamental solutions for elliptic boundary value problems. Adv. Comput. Math. 1998, 9, 69–95. [Google Scholar] [CrossRef]
Hoang, H.T.; Izsák, F.; Maros, G. Method of fundamental solutions: New approximation results and applications. J. Comput. Appl. Math. 2024, 448, 115934. [Google Scholar] [CrossRef]
Kálnay de Rivas, E. On the use of nonuniform grids in finite-difference equations. J. Comp. Phys. 1972, 10, 202–210. [Google Scholar] [CrossRef]
Gamet, L.; Ducros, F.; Nicoud, F.; Poinsot, T. Compact finite difference schemes on non-uniform meshes. Application to direct numerical simulations of compressible flows. Int. J. Numer. Meth. Fluids 2022, 29, 159–191. [Google Scholar] [CrossRef]
Rahman, T.; Valdman, J. Fast MATLAB assembly of FEM matrices in 2D and 3D: Nodal elements. Appl. Math. Comput. 2013, 219, 7151–7158. [Google Scholar] [CrossRef]
Cuvelier, F.; Japhet, C.; Scarella, G. An efficient way to assemble finite element matrices in vector languages. Bit Numer. Math. 2016, 56, 833–864. [Google Scholar] [CrossRef]
Fu, Z.; Lewis, T.J.; Kirby, M.R.; Whitaker, T.R. Architecting the finite element method pipeline for the GPU. J. Comput. Appl. Math. 2014, 257, 195–211. [Google Scholar] [CrossRef]
Krysl, P. Parallel assembly of finite element matrices on multicore computers. Comput. Meth. Appl. Mech. Eng. 2024, 428, 117076. [Google Scholar] [CrossRef]
Izsák, F.; Izsák, R. Neural-Network-Assisted Finite Difference Discretization for Numerical Solution of Partial Differential Equations. Algorithms 2023, 16, 410. [Google Scholar] [CrossRef]
Lunardini, V.J. Heat Transfer with Freezing and Thawing; Elsevier: Amsterdam, The Netherlands, 1991. [Google Scholar]
Visintin, A. Handbook of Differential Equations: Evolutionary Equations; Dafermos, M.C., Milan Pokorny, M., Eds.; North-Holland: Amsterdam, The Netherlands, 2008; Volume 4, Chapter 8; p. 377. [Google Scholar]
Vasil’ev, V.; Vasilyeva, M. An Accurate Approximation of the Two-Phase Stefan Problem with Coefficient Smoothing. Mathematics 2020, 8, 1924. [Google Scholar] [CrossRef]
Caldwell, J.; Kwan, Y.Y. Numerical methods for one-dimensional Stefan problems. Commun. Numer. Meth. Eng. 2004, 20, 535–545. [Google Scholar] [CrossRef]
Dalhuijsen, A.J.; Segal, A. Comparison of finite element techniques for solidification problems. Int. J. Numer. Meth. Eng. 1986, 23, 1807–1829. [Google Scholar] [CrossRef]
Segal, G.; Vuik, K.; Vermolen, F. A Conserving Discretization for the Free Boundary in a Two-Dimensional Stefan Problem. J. Comput. Phys. 1998, 141, 1–21. [Google Scholar] [CrossRef]
Bommes, D.; Levy, B.; Pietroni, N.; Puppo, E.; Silva, C.; Zorin, D. Quad-Mesh Generation and Processing: A Survey. Comput. Graph. Forum. 2013, 32, 12014. [Google Scholar] [CrossRef]
Patra, M.; Karttunen, M. Stencils with isotropic discretization error for differential operators. Numer. Methods Partial Differ. Eq. 2006, 22, 936–953. [Google Scholar] [CrossRef]
Brosa Planella, F.; Please, C.P.; van Gorder, R.A. Extended Stefan problem for the solidification of binary alloys in a sphere. Eur. J. Appl. Math. 2021, 32, 242–279. [Google Scholar] [CrossRef]
Escobedo, R.; Fernández, L.A. Classical one-phase Stefan problems for for describing polymer crystallization processes. SIAM J. Appl. Math. 2013, 73, 254–280. [Google Scholar] [CrossRef]
Tubini, N.; Gruber, S.; Rigon, R. A method for solving heat transfer with phase change in ice or soil that allows for large time steps while guaranteeing energy conservation. Cryosphere 2021, 15, 2541–2568. [Google Scholar] [CrossRef]
Zeneli, M.; Nikolopoulos, A.; Karellas, S.; Nikolopoulos, N. Ultra-High Temperature Thermal Energy Storage, Transfer and Conversion; Datas, A., Ed.; Woodhead Publishing: Cambridge, UK, 2021; Chapter 7; p. 165. [Google Scholar]
Fullana, T.; Le Chenadec, V.; Sayadi, T. Adjoint-based optimization of two-dimensional Stefan problems. J. Comput. Phys. 2023, 475, 111875. [Google Scholar] [CrossRef]
Brink, F.; Izsák, F.; van der Vegt, J.J.W. Hamiltonian Finite Element Discretization for Nonlinear Free Surface Water Waves. J. Sci. Comput. 2017, 73, 366–394. [Google Scholar] [CrossRef]
Arnold, D.N.; Boffi, D.; Falk, R.S. Approximation by quadrilateral finite elements. Math. Comput. 2002, 71, 909–922. [Google Scholar] [CrossRef]
Dias, F.; Dutykh, D. Dynamics of tsunami waves. In Extreme Man-Made and Natural Hazards in Dynamics of Structures; Ibrahimbegovic, A., Kozar, I., Eds.; NATO Security Through Science Series; Springer: Dordrecht, The Netherlands, 2007; p. 201. [Google Scholar]

Figure 1. A uniform (left) and a general (right) quadrilateral grid for finite difference approximations. A uniform rectangular grid (left) related to the approximation (4) and a non-uniform quadrilateral grid (right) related to (5).

Figure 2. Location of non-zero elements in the discretization matrix

D_{h, FD}

corresponding to (4) (left) and (5) (right), respectively. The total matrix size

36 \times 36

corresponds to the number of the interior grid points. The notation nz yields the number of non-zeros.

Figure 2. Location of non-zero elements in the discretization matrix

D_{h, FD}

corresponding to (4) (left) and (5) (right), respectively. The total matrix size

36 \times 36

corresponds to the number of the interior grid points. The notation nz yields the number of non-zeros.

Figure 3. Spatial discretization using structured grids in two cases. Number of uniformly distributed grid points: 20 (left) and 40 (right) in the horizontal direction. Number of grid points: 16 (left) and 33 (right) in the vertical direction.

Figure 4. Motion of the segments

x_{i - 1} x_{i}, x_{i} x_{i + 1},

and

x_{i + 1} x_{i + 2}

Γ^{n}

in step (ii) (a).

d_{i}

denotes the shift vectors obtained from the jump term in (26). Dashed lines show the shifting of the segments of

Γ^{n}

. Their average on the vertical lines gives the new position of the vertices of

Γ^{n + 2}

Figure 4. Motion of the segments

x_{i - 1} x_{i}, x_{i} x_{i + 1},

and

x_{i + 1} x_{i + 2}

Γ^{n}

in step (ii) (a).

d_{i}

denotes the shift vectors obtained from the jump term in (26). Dashed lines show the shifting of the segments of

Γ^{n}

. Their average on the vertical lines gives the new position of the vertices of

Γ^{n + 2}

Figure 5. The steps of shifting the interface

Γ^{n}

to obtain the new one

Γ^{n + 1}

Figure 5. The steps of shifting the interface

Γ^{n}

to obtain the new one

Γ^{n + 1}

Figure 6. Initial condition (left, with a heat map) and the corresponding initial grid (right) for the simulation of the problem in (11) and (12). Bottom region:

Ω_{1, t}

, depicted with red color (left) and red grid points (right); top region:

Ω_{0, t}

, depicted with blue color (left) and blue grid points (right). The common interface and the corresponding grid points are colored with green. The horizontal distribution of grid points is uniform, while the vertical distribution becomes more dense towards the common interface.

Figure 6. Initial condition (left, with a heat map) and the corresponding initial grid (right) for the simulation of the problem in (11) and (12). Bottom region:

Ω_{1, t}

, depicted with red color (left) and red grid points (right); top region:

Ω_{0, t}

Figure 7. The result of the simulation procedure for the Stefan problem after 7 time steps (left) and 17 time steps (right), respectively. The corresponding parameters are given in the list below. Again, blue color represents negative temperature, while red indicates positive temperature.

Table 1. Simulation time (ms) and discrete

L_{2}

-norm error for various numbers of horizontal (

N_{x}

) and vertical (

N_{y}

) grid points in the discretization of the model problem (10).

Table 1. Simulation time (ms) and discrete

L_{2}

-norm error for various numbers of horizontal (

N_{x}

) and vertical (

N_{y}

) grid points in the discretization of the model problem (10).

$N_{x}$	10	15	20	25	30	35	40	45	50	55
$N_{y}$	8	12	16	20	24	28	32	36	40	44
time (ms)	22	20	27	41	55	73	95	119	146	174
error $\cdot 10^{3}$	10.9	5.2	3.1	2.1	1.6	1.2	1	0.87	0.77	0.7

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Castillo, S.-J.; Izsák, F. Efficient Discretization of the Laplacian: Application to Moving Boundary Problems. Algorithms 2024, 17, 541. https://doi.org/10.3390/a17120541

AMA Style

Castillo S-J, Izsák F. Efficient Discretization of the Laplacian: Application to Moving Boundary Problems. Algorithms. 2024; 17(12):541. https://doi.org/10.3390/a17120541

Chicago/Turabian Style

Castillo, Sebastian-Josue, and Ferenc Izsák. 2024. "Efficient Discretization of the Laplacian: Application to Moving Boundary Problems" Algorithms 17, no. 12: 541. https://doi.org/10.3390/a17120541

APA Style

Castillo, S. -J., & Izsák, F. (2024). Efficient Discretization of the Laplacian: Application to Moving Boundary Problems. Algorithms, 17(12), 541. https://doi.org/10.3390/a17120541

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Efficient Discretization of the Laplacian: Application to Moving Boundary Problems

Abstract

1. Introduction

2. Materials and Methods

2.1. Problem Statement and Mathematical Background

2.2. The Pointwise Optimization

2.3. Numerical Experiments

2.3.1. The Model Problem

2.3.2. Spatial Discretization

2.3.3. Numerical Results

2.4. Application to the Numerical Solution of the Stefan Problem

2.4.1. General Setting

2.4.2. Specific Initial and Boundary Conditions

2.4.3. Construction of the Grid

2.5. Principles of the Algorithm for Solving (11)–(13)

2.5.1. STEP (i) (b)

2.5.2. STEP (i) (c)

2.5.3. STEP (ii) (a)

2.5.4. STEP (ii) (b)

2.5.5. Numerical Results

3. Discussion and Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

Abbreviations

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI