Deep neural operator for free boundary problems

Problem set-up

FBP formulation

Consider a family of evolving domains \({\{{\varOmega }_{t}\}}_{t\in [0,T]}\), where each \({\varOmega }_{t}\subset {{\mathbb{R}}}^{s}\) is a measurable compact set representing the spatial region at time t. We consider the effective spatio-temporal domain of the FBP to consist of concatenated sets \(\varOmega \,:={\cup }_{t\in [0,T]}{\varOmega }_{t}\times \{t\}\subset {{\mathbb{R}}}^{s}\times [0,T]\), where Ω is a (s + 1)-dimensional bounded closed set. The free boundary Γ := ⋃_{t∈[0, T]}∂Ω_t × {t} is assumed to be a compact Lipschitz hypersurface in \({{\mathbb{R}}}^{s-1}\times [0,T]\). Unlike a classical PDE where the domain Ω and Γ are given, the main feature of the FBP is the free motion property of the boundary ∂Ω_t. This motion is intrinsically coupled with the variation of internal physical quantities that collectively determine the solution. Consequently, both Ω and Γ exhibit a priori unknown properties, with solely the initial configuration Ω₀ and requisite input conditions being given, as illustrated in Fig. 1a. This intrinsic geometry–physics coupling is epitomized by the classical Stefan problem, as detailed in ‘Illustration of the Stefan problem’.

Description of the solution operator

Within the realm of FBPs, the solution of the FBP family we attempt to obtain is

which is a family of Sobolev spaces, where \(\mathrm{Comp}({{\mathbb{R}}}^{s+1})\) denotes compact subsets of \({{\mathbb{R}}}^{s+1}\). Here the superscript s + 1 in U^s+1( ⋅ ) indicates that the domain of the function u ∈ U^s+1( ⋅ ) is embedded in \({{\mathbb{R}}}^{s+1}\). Mathematically, we can derive an accurate solution operator as depicted in Fig. 1a:

where U₀(Ω₀) denotes the Banach space of initial condition functions (for example, temperature distributions), and ψ ∈ Ψ refers to an auxiliary function defined on the domain \({\mathcal{A}}\) other than Ω₀ (for example, heat sources). U₀ and Ψ can be realized under different frameworks, such as Sobolev spaces W^k,2(Ω₀) and \({W}^{k,2}({\mathcal{A}})\), respectively.

The foundational universal approximation theorem (UAT) for neural operators, such as DeepONet¹⁵, is established in ref. ³⁵. A key assumption in this work is that the codomain of the operator must be C(K), meaning the output functions are continuous over a prescribed compact set K. This requirement ensures that the operator satisfies the conditions of the Tietze extension theorem³⁵, uniform continuity and boundedness, properties essential in the proof of the universal approximation. However, for the FBP solution operator \({\mathcal{F}}\), the codomain U^s+1( ⋅ ) consists of functions defined on domains that are not fixed a priori. This lack of a prescribed compact domain introduces a fundamental incompatibility with the assumptions of the UAT, making the existing theory inapplicable to \({\mathcal{F}}\).

Methodology

Overview

The solution operator pertaining to the FBP family can be envisaged as an infinite-dimensional dynamical system. The phase space is denoted by U^s( ⋅ ) and can be analogously defined as in equation (1). The flow map \({{\mathcal{F}}}^{t}={\mathcal{F}}({u}_{0}| \psi )(\cdot ,t)\) represents the solution operator at time t. It is natural to construct a conjugate dynamical system over a fixed reference domain Ω_ref that is diffeomorphic to Ω_t. The conjugate dynamics are governed by a flow map \({{\mathcal{G}}}^{t}\) acting on a phase space V(Ω_ref) that consists of functions v_t( ⋅ ) defined on Ω_ref. The key is to establish a topological conjugacy between \({{\mathcal{F}}}^{t}\) and \({{\mathcal{G}}}^{t}\) via a conjugate homeomorphism \({\mathcal{H}}\) such that \({{\mathcal{F}}}^{t}={\mathcal{H}}\circ {{\mathcal{G}}}^{t}\circ {{\mathcal{H}}}^{-1}\) for t∈[0,T]. This facilitates the indirect computation of the dynamics of the original system \({{\mathcal{F}}}^{t}\) through \({{\mathcal{G}}}^{t}\) subsequent to the transformation by \({\mathcal{H}}\). Figure 1c is a schematic representation.

Operator construction

To ensure that the target operator satisfies the assumptions of the UAT, it is imperative that v_t exhibits continuity, that the input function is established a priori and that operator \({\mathcal{H}}\) exhibits reversibility. Under these stipulated conditions, we employed these subsequent constructions:

where \({\mathcal{M}}:{u}_{0}| {v}_{t},\psi \mapsto {\chi }_{t}^{-1}\) generates the inverse diffeomorphism \({\chi }_{t}^{-1}:\,{\varOmega }_{t}\to {\varOmega }_{\mathrm{ref}}\), which maps the evolving domain Ω_t to the reference domain Ω_ref. Also, \({{\mathcal{M}}}^{{\prime} }:{u}_{t}| {u}_{0},\psi \mapsto {\chi }_{t}\). It can be verified that \({\mathcal{H}}\) is invertible: \({\mathcal{H}}\circ {{\mathcal{H}}}^{-1}({u}_{t}| {u}_{0},\psi )\)\(={\mathcal{H}}[{u}_{t}\circ {{\mathcal{M}}}^{{\prime} }({u}_{t}| {u}_{0},\psi )({\xi })| {u}_{0},\psi ]\)\(={u}_{t}\circ {\chi }_{t}\circ {\chi }_{t}^{-1}={u}_{t}\), which shows that the construction is correct. Nevertheless, according to the UAT, it is not feasible to approximate \({\mathcal{H}}\), \({{\mathcal{G}}}^{t}\) or \({{\mathcal{H}}}^{-1}\) using neural networks. Consequently, we propose the application of indirect representation operators:

where χ(ξ, t) = χ_t(ξ). Ultimately, it is presumed that the operators exhibit continuity with respect to time, as indicated by \({\mathscr{G}}({u}_{0},\psi )({\xi },\tau )={{\mathscr{G}}}^{\tau }({u}_{0},\psi )({\xi })\), where \({\mathscr{G}}:{u}_{0},\psi \mapsto v({\xi },t)\). Consequently, the solution operator can be approximated by using neural operators \(\widehat{{\mathscr{G}}}\) and \(\widehat{{\mathscr{H}}}\):

To ensure that the output χ(ξ, t) of \(\widehat{{\mathscr{H}}}\) constitutes a diffeomorphism, we impose the diffeomorphic constraint detailed in ‘Diffeomorphic constraint’. The detailed conceptual framework of the construction is provided in Supplementary Section A. Given that the UAT cannot be directly applied to \({\mathcal{F}}\) when the solution functions are defined on predefined domains, we present the following approximation theorem tailored to our methodology:

Theorem 1

(UAT for FBPs) Let X₁(Σ) and X₂(Ω₀) be Banach function spaces defined on Σ and Ω₀, respectively. Suppose K_i ⊂ X_i(i = 1, 2) are compact sets, and X_i(i = 1, 2) have a Schauder basis. U^s+1( ⋅ ) is defined in equation (1), and K_u ⊂ U^s+1( ⋅ ) is a compact set. Assuming that for all the domains of the solution function Ω:

For an operator \({\mathcal{F}}:{K}_{1}\times {K}_{2}\to {K}_{u}\), then for any ϵ > 0, there exist two neural networks \(\widehat{{\mathscr{G}}}\) and \(\widehat{{\mathscr{H}}}\) with the forms given in equation (24), equation (25), respectively, such that their composition \(\widehat{{\mathcal{F}}}\) as defined in equation (7), satisfies

for all ψ ∈ K₁, u₀ ∈ K₂ and (x, t) ∈ Ω.

A complete proof and discussion are provided in Supplementary Section C.

Model architecture

Consider a collection of n observations \({\{{u}_{0}^{i},{\psi }^{i}\}}_{i=0}^{n}\), each of which is randomly sampled within its function spaces \({U}_{0}({\varOmega }_{0})\) and \(\varPsi (A)\), respectively, where uⁱ ∈ U^s+1( ⋅ ) represents the target solutions. It is noteworthy that the Ω₀ may vary among different observations, which will be elaborated upon in ‘Adaptive handling of the initial domain’. The training process involves simultaneously learning the two models \(\widehat{{\mathscr{G}}}\) and \(\widehat{{\mathscr{H}}}\) in equation (7). The model \(\widehat{{\mathscr{H}}}\) operates on parameterized representations of the input data \({\{{u}_{0}^{i},{\psi }^{i}\}}_{i=0}^{n}\) augmented with a spatial coordinate ξ sampled from the reference domain Ω_ref and a temporal coordinate t ∈ [0, T] to generate a diffeomorphism χ(ξ, t). Concurrently, the model \(\widehat{{\mathscr{G}}}\) processes the same parameterized inputs, along with ξ and t, to compute the conjugate phase field vⁱ(ξ, t). Upon completion of the training processes for both models, they should be amalgamated as delineated in equation (7) to derive the solution operator \(\widehat{{\mathcal{F}}}\). The application of this operator for domain evaluation and validation is discussed in ‘Evaluation of the solution operator’.

The theoretical framework for FBNO is decoupled from specific architectural implementations of neural operators, subject to the fulfilment of the key requirements above. To maintain theoretical generality, we employ MIONet as the foundational model, as detailed in ‘Neural operator’. Additionally, we developed an architecture that enables the embedding of χ to be embedded in \(\widehat{{\mathscr{G}}}\), although this is not strictly necessary, as shown in Fig. 1d.

Numerical results

Stefan problem

Problem statement

The Stefan problem is a classical FBP that models phase transitions driven by heat diffusion and latent heat effects, such as melting or solidification. It originated in the late 19th century with the study of ice formation^36,37. A mathematical framework was developed to describe the growth of sea ice in polar regions. We consider a Stefan problem of the form:

where Γ_f = {(x, t) ∈ Ω∣x = s(t)} and Γ_d = Γ\Γ_f. The temperature field is denoted by u. The free boundary separating the two phases is represented by s(t), and Ω₀ is defined by a given s(0). The thermal diffusivity is given by the parameter α, and the spatially dependent heat source term is modelled as Q(x). The parameter k represents the thermal conductivity of the material. Additionally, C(t) accounts for the effect of heat absorption, and v_n denotes the normal velocity of the free boundary s(t). The solution u and the free boundary Γ_f are determined under appropriate initial conditions u₀ and Dirichlet boundary conditions u_c and u_b. The source term Q(x) and the absorption term C(t) are modelled as random fields, as elaborated in ‘Configuration of the numerical experiment’.

Method

In this example, the mapping χ( ⋅ , t) constitutes a diffeomorphism provided the spatial Jacobian determinant remains positive. Thus, specialized bijection network architectures are unnecessary in this setting. The training process relies solely on a physics-informed method. To ensure numerical stability and prevent local over-expansion or compression, the training process jointly optimizes both interior points and free boundary points. The free boundary condition χ(s(0), τ) = s(τ) is enforced. Interior points at time τ are constrained such that: \(\delta < | {J}_{{\mathcal{S}}}| < s(\tau )\), where δ is a small positive constant. Figure 2d highlights how the determinant bound \(| {J}_{{\mathcal{S}}}| < s(\tau )\) mitigates excessive local expansion.

Network performance

Figure 2a displays the solutions predicted by the neural network u_pred for three test cases, along with the ground truth solutions u_ref obtained via a finite element solver. The figure further presents the absolute error ∣u_ref − u_pred∣ and the pointwise relative absolute error for the free boundary, defined as \(\frac{|{{\varGamma }_{{\rm{f}}}}_{{\rm{ref}}}-{{\varGamma }_{{\rm{f}}}}_{{\rm{pred}}}|}{|{{\varGamma }_{{\rm{f}}}}_{{\rm{ref}}}|}\). The results demonstrate the capability of FBNO to accurately infer both the domain Ω and the associated temperature field, even in the absence of empirical training data. A further quantitative assessment is provided in Fig. 2c, which illustrates the relative L² errors \(\frac{\parallel {u}_{\mathrm{ref}}-{u}_{\mathrm{pred}}\parallel }{\parallel {u}_{\mathrm{ref}}\parallel }\) across different forcing functions. Given that cases where \(| {J}_{{\mathcal{S}}}| < 0\) constitute fewer than 1 × 10⁻¹⁰ in total, the neural network can be considered to satisfy the bijection constraint with near-perfect accuracy. Consequently, the relative L² errors of Γ_f serve as a reliable metric for assessing the precision of χ. The results indicate that most of the errors in the temperature field u are within 1.5%, whereas those in Γ_f are predominantly below 1%. The marginally higher error in u was anticipated, given that u represents the temperature field acting on Ω, whereas Γ_f is a parameter that characterizes the domain Ω itself. To further assess the performance of FBNO, the domain distance d(Ω, Ω_ref) between Ω and Ω_ref was quantified using a normalized Riemannian metric, as detailed in ‘Normalized Riemannian metric’. As depicted in Fig. 2b, the resulting losses for both u and Γ_f exhibit negligible sensitivity to the domain geometry, with deviations not exceeding 0.03%. These findings collectively confirm that FBNO effectively captures the functional mapping between various domain geometries.

Thermal–structural coupling problem

Problem statement

To further demonstrate the capability of FBNO in solving complex FBPs, we tested a more complicated problem with several physical fields: a thermal–structural coupling problem. This system describes the thermal expansion, heat transfer and heat convection processes. We consider equations of the following form:

The system describes the evolution of density ρ and temperature T, where ρ is determined by a given function F(T) of temperature. The velocity v is governed by a continuity equation and a heat equation with advection. Here c_p is the specific heat capacity, λ is the thermal conductivity and Q(x, t) is a randomly sampled heat source, as detailed in ‘Configuration of the numerical experiment’. The boundary conditions include a free boundary Γ_f = {(x, t) ∈ Ω∣x = s(t)} with normal velocity \({v}_{{\rm{n}}}{| }_{{\varGamma }_{{\rm{f}}}}=v\), a fixed temperature T_b, a prescribed velocity v_b on Γ_d = Γ⧹Γ_f, and an initial condition T₀ at t = 0 with Ω₀ = [0, s(0)]. The variables to be solved for are the density ρ(x, t), temperature T(x, t), velocity v(x, t) and the free boundary Γ_f. All requisite data were generated following the methodology presented in ‘Configuration of the numerical experiment’.

Method

To illustrate the high flexibility and inclusiveness of FBNO, under the premise of mass conservation ∂_tρ + ∂_xρv = 0, we introduce particle tracking via the diffeomorphism χ. The essence of this method lies in the arbitrary Lagrangian–Eulerian framework, which maps the particles from the initial moment to their evolved positions at subsequent times, thereby enabling the extraction of the corresponding physical fields. To ensure that the neural network \(\widehat{{\mathscr{H}}}\) approximates this mapping effectively, we impose the constraint:

where ρ₀ = ρ(χ⁻¹(x, 0), 0) = ρ(ξ, 0) represents the initial density distribution. We can reasonably assume that \(\rho \in \left\{g:\varOmega \to {\mathbb{R}}\right.\)\(\left.| \exists M\in {\mathbb{R}},\,\forall {x}\in \varOmega ,\,0\le g({x})\le M\,\right\}\). Then χ is naturally a diffeomorphism. Consequently, the neural network need not enforce bijectivity explicitly. This approach can also be interpreted as augmenting the original system with another PDE defined on a fixed reference domain Ω_ref. In this example, the models were trained using a mere ten low-resolution simulations, in conjunction with physics-informed constraints. Performance was assessed on a test set comprising 200 high-resolution simulations.

Network performance

Given that the diffeomorphism χ is uniquely determined, we treated χ itself as a physical field. Figure 3a presents a comparative visualization of the ground truth solution u_ref, the neural network prediction u_pred and the corresponding absolute error ∣u_ref − u_pred∣. The absolute error field demonstrates that the neural network achieves high fidelity in approximating the physical fields. As depicted in Fig. 3e, a quantitative analysis of the relative L² error distributions for FBNO across multi-physics fields reveals a consistent positive skewness. All distributions exhibit positive skewness, indicative of a systematic bias towards lower error values, with mean relative error confirming that most errors reside within a satisfactory tolerance range. Furthermore, as illustrated in Fig. 3b,c, the incorporation of physics-informed constraints substantially enhances model performance, even when most training errors alone are far less favourable compared with purely data-driven approaches. Notably, this method achieves a reduction of more than 95% in testing errors across all physical quantities relative to exclusively data-based approaches, demonstrating its efficacy in mitigating data scarcity. Additionally, a temporal stability analysis (Fig. 3d) demonstrated the robustness of the model over extended simulation periods. The results confirm that prediction errors remain bounded within consistent margins throughout the temporal domain, with no evidence of substantial error accumulation. Furthermore, experiments conducted in two-dimensional settings, as detailed in Supplementary Section G, collectively demonstrate that FBNO achieves high efficiency and accurate predictive performance.

Tumour growth problem

Problem statement

Here we demonstrate the capability of FBNO to handle non-convex geometries, with potential applications in tumour growth modelling. Tumour growth kinetics presents a complex challenge due to the intricate interplay between cellular proliferation and microenvironmental nutrient dynamics, and it is governed by highly nonlinear interactions^38,39. To generate a comprehensive dataset encompassing arbitrary growth durations and repeated trials, we simulated tumour expansion within surrounding tissue (in vivo) or passive polymer gels (in vitro) using the established mathematical framework from ref. ⁴⁰. The governing dynamics incorporate a forcing function proposed by ref. ⁴¹, which describes the net cell proliferation rate:

where ϵ denotes the tumour volume fraction, α, β, λ and γ are kinetic parameters regulating birth and death rates, and u represents the nutrient concentration that needs to be solved. The first term models mitotically driven growth, and the second captures nutrient-dependent apoptosis. The various death rates engender distinct dynamic processes, as illustrated in Supplementary Section E.2.1. Our aim is to ascertain the nutrient distribution u by employing various governing equations alongside diverse initial domain configurations Ω₀. Thus, we randomly sample f and examine three representative initial domains Ω₀, as depicted in ‘Configuration of the numerical experiment’. Next, we simulate the system dynamics, as detailed in ‘Numerical simulations of tumour growth’.

Method

To comprehensively demonstrate the capability of FBNO in modelling dynamical systems with incomplete physical information, we adopt a purely data-driven approach. Furthermore, to ensure the physical plausibility of biological tissue deformation and prevent unphysiological folding or tearing artefacts, we incorporate another constraint via the reconciliation mapping, expressed as ∥Δχ∥, while maintaining \(| {J}_{{\mathcal{S}}}| > \delta\). Given the sufficiency of the available data, we perform pretraining for the neural operator \(\widehat{{\mathscr{H}}}\). To guarantee that the output function of \(\widehat{{\mathscr{H}}}\) is bijective, we impose an autoencoder-based regularization \(\parallel \widehat{{\mathscr{H}}}(\cdot )\circ \widehat{{\mathscr{H}}}{(\cdot )}^{-1}-I\parallel\) on the neural network, following the methodology introduced by ref. ⁴². Moreover, as the initial domain geometries vary across examples in this study, we employ a parametric representation of the initial domains, which are subsequently integrated into the neural network using the approach described in ‘Adaptive handling of the initial domain’.

Network performance

The evaluation of the pretrained model \(\widehat{{\mathscr{H}}}\) is shown in Fig. 4a. The figure illustrates the reference free boundary Γ, model predictions and absolute errors across three representative test cases. The maximum absolute error observed in all instances confirms that our framework accurately constructs the domains. We further analyse the performance of the model sequence in Supplementary Section E.1.1. Figure 4b presents the distribution of pointwise absolute errors for various initial domain configurations, as detailed in Supplementary Section E.2.2. Notably, the error distributions exhibit a consistent shape and dispersion range, regardless of the geometry of the initial domain (including non-convex structures). This demonstrates that the approximating capability of FBNO is invariant to the topological complexity of the input domain. As depicted in Fig. 4c, the mean error curve and 95% confidence interval remain stable over time, with fluctuations confined to a narrow range. Concurrently, the green bar plot highlights the monotonic increase in tumour surface area during growth. These findings indicate that FBNO maintains consistent predictive accuracy across both spatial and temporal scales, underscoring its robustness to dynamic morphological changes. For each quantity, the distributions of both the relative L² error and the mean squared error are illustrated. The right-skewed error distributions, coupled with satisfactory mean error values, indicate that FBNO is more likely to produce predictions better than the average error threshold. This further validates the strong generalizability of the model in capturing tumour growth dynamics. Figure 4d highlights the computational advantages of FBNO. After training, FBNO achieves a speed-up factor of 10⁴ compared with traditional numerical methods, with memory and energy efficiency ratios of the order of 7 and 5 × 10³, respectively. Specifically, FBNO predicts several tumour growth trajectories in seconds on a single graphics processing unit, whereas equivalent simulations require over 2 days on 16 central processing unit nodes. This notable efficiency gain enables real-time tumour growth prediction, with meaningful implications for clinical applications.