Multilayer Level-Set Method (MLSM)

• MLSM is a computational framework that represents multiple interfaces in multiphase systems using a single scalar level-set function with prescribed thresholds.
• It leverages advanced numerical schemes and deep learning techniques to accurately compute curvature and normals, enhancing simulation fidelity in complex geometries.
• MLSM is applied in geophysical imaging, multiphase flow, and topology optimization, providing efficient interface reconstruction and robust error correction.

The Multilayer Level-Set Method (MLSM) extends classical level-set frameworks to represent and compute evolving interfaces in systems containing more than two regions or phases. Traditional level-set methods identify interfaces by the zero level-set of a scalar field; in contrast, MLSM introduces multiple prescribed level values (e.g., $\{i_n \mid n = 0, 1, ..., N\}$) within a single level-set function, allowing simultaneous evolution and recovery of an arbitrary number of interfaces and subregions. MLSM is applicable to free-boundary problems, multiphase systems, topology optimization, and inverse imaging problems where rapid topological changes, discontinuities, or piecewise-constant fields occur.

1. Conceptual Foundations of MLSM

The key innovation of MLSM is encoding the phase structure via a single scalar field $\phi(x)$, in which each threshold $i_n$ defines an interface $\phi(x) = i_n$. Adjacent intervals between threshold values correspond to distinct physical regions or phases. This generalizes two-phase representations without introducing multiple independent level-set functions, simplifying both the mathematical formalism and the computational cost in problems with $N > 2$ regions.

In inverse imaging (e.g., eikonal-based traveltime tomography), physical model parameters (such as slowness $S(x)$) are represented by compositions of Heaviside functions applied to the multilayer level set. A typical formula in MLSM is

$$S(x) = p_0(x)[1 - H(\phi(x) - i_0)] + p_1(x)[H(\phi(x) - i_0) - H(\phi(x) - i_1)] + \cdots + p_N(x) H(\phi(x) - i_{N-1})$$

where $H$ is a smoothed Heaviside function and the $p_n(x)$ parameterize phase-specific properties (Li et al., 18 Oct 2025).

2. Interface Geometry and Numerical Schemes

Accurate estimation of interface curvature and normals is critical for MLSM, especially in simulations of multiphase flow, interface-driven physics, and topology optimization. Standard discretizations derive curvature as $\kappa = \nabla \cdot (\nabla\phi / |\nabla\phi|)$, but can fail near "kinks" where the signed-distance property of $\phi$ breaks down—leading to large spurious errors near closely spaced interfaces (Lervåg, 2014).

Improved schemes detect such regions using quality measures (e.g., $Q(x) = |1 - |\phi(x)||$) to identify non-smooth zones and apply:

• Directionally one-sided finite differences for normal computations;
• Local spline-based curve fitting (e.g., monotone cubic Hermite splines) to reconstruct a smooth local interface and enable robust curvature evaluation.

These techniques minimize topological artifacts in simulations with frequent interface merging and breakup, essential for robust MLSM applications.

3. Deep Learning and Hybrid Inference for Curvature Estimation

Recent developments have incorporated supervised feedforward neural networks (multilayer perceptrons) to estimate curvature directly from local stencils of level-set values (2002.02804, Larios-Cárdenas et al., 2021, Larios-Cárdenas et al., 2022, Larios-Cárdenas et al., 2022). The general framework is:

• The numerical scheme produces an initial curvature $h\kappa$ on a stencil.
• A neural network, trained on synthetic data from exact solutions (circles, sinusoidal surfaces, hyperbolic paraboloids), predicts the systematic error $\bar{\varepsilon}$ in this estimate.
• The corrected curvature is $h\kappa^* = h\kappa + \bar{\varepsilon}$, suppressing numerical deficiencies.

For MLSM on adaptive or coarse grids, this hybrid method is activated selectively in regions where the numerical error is likely high (e.g., curvature magnitude above a resolution-dependent threshold). Data preprocessing methods—rotational canonicalization, sign normalization, reflection-based augmentation—are essential for training efficiency and model accuracy, with regularization and dimensionality reduction (e.g., via PCA) minimizing outlier effects (Larios-Cárdenas et al., 2022, Larios-Cárdenas et al., 2022).

4. Multiphase Representations and Topology Optimization

Extensions of MLSM to multiphase topology optimization problems can be found in the “extended level set” (X-LS) method (Noda et al., 2021). Here, interfaces between every pair of materials $(i,j)$ are represented by zero-level surfaces of antisymmetric level-set functions $\phi_{ij} = -\phi_{ji}$, with phase membership encoded by

$\psi_m(x) = \prod_{i \neq m} H(\phi_{im}(x))$

The evolution of interfaces is governed by reaction-diffusion equations incorporating extended topological derivatives and Laplacian regularization, solved in a finite element setting. The antisymmetric property ensures perfect interface equivalence and unbiased phase treatment, reducing initial condition dependence and simplifying sensitivity analysis in optimization routines.

5. Regularization in MLSM Inverse Imaging Applications

In MLSM-based inverse problems (e.g., eikonal tomography (Li et al., 18 Oct 2025)), the stability and accuracy of interface recovery depend on several regularization strategies:

• Multilayer Reinitialization: Level-set function $\phi$ is reinitialized near each prescribed threshold $i_n$ to maintain the signed-distance property locally.
• Arc-Length Penalization: Interface complexity is discouraged with terms of the form $$E_r(\phi) = \int_{\Omega} \sum_n \delta(\phi - i_n) |\nabla \phi| \, dx$$, where $\delta$ approximates the Dirac measure.
• Sobolev Smoothing: Phase parameter updates $p_n(x)$ are regularized via PDEs, e.g., $(I - \gamma \Delta)P = p$, to encourage smooth field reconstructions.

6. Error Assessment and Reconstruction Quality

MLSM recognizes that partial-ray coverage ("illumination") may hinder reliable inversion in certain domains. An illumination-based error measure weights local errors according to a labeling $F(x)$, quantifying the density/energy of traversing rays. The error metric

$$e_F(x) = F(x) \cdot |S_{\text{recovered}}(x) - S_{\text{true}}(x)|$$

enables physically meaningful assessment, more stringent where data is abundant and more tolerant in shadowed regions.

7. Practical Impact and Applications

MLSM finds strong application in geophysical imaging (seismic tomography), multiphase flow simulation, solidification modeling, tumor growth tracking, and multi-material topology optimization. By encoding arbitrary phase numbers with a single function and leveraging advanced geometric, numerical, and machine learning techniques for curvature and interface regularization, MLSM enables robust and computationally efficient treatment of problems featuring complex geometry and rapid topological change.

Recent numerical experiments demonstrate that MLSM approaches can efficiently recover discontinuous fields with multiple interfaces, outperforming baseline methods in accuracy and computational cost for under-resolved regions. The synergy of error-correcting neural networks, adaptive regularization, and geometric processing marks MLSM as an extensible framework for next-generation computational physics and inverse problems (Li et al., 18 Oct 2025, Larios-Cárdenas et al., 2022, Larios-Cárdenas et al., 2022, Noda et al., 2021, Lervåg, 2014).