A Physics-Informed Meta-Learning Framework
for the Continuous Solution of Parametric PDEs
on Arbitrary Geometries

Abstract

-1cm-1cm In this work, we introduce implicit Finite Operator Learning (iFOL) for the continuous and parametric solution of partial differential equations (PDEs) on arbitrary geometries. We propose a physics-informed encoder-decoder network to establish the mapping between continuous parameter and solution spaces. The decoder constructs the parametric solution field by leveraging an implicit neural field network conditioned on a latent or feature code. Instance-specific codes are derived through a PDE encoding process based on the second-order meta-learning technique. In training and inference, a physics-informed loss function is minimized during the PDE encoding and decoding. iFOL expresses the loss function in an energy or weighted residual form and evaluates it using discrete residuals derived from standard numerical PDE methods. This approach results in the backpropagation of discrete residuals during both training and inference.

iFOL features several key properties: (1) its unique loss formulation eliminates the need for the conventional encode-process-decode pipeline previously used in operator learning with conditional neural fields for PDEs; (2) it not only provides accurate parametric and continuous fields but also delivers solution-to-parameter gradients without requiring additional loss terms or sensitivity analysis; (3) it can effectively capture sharp discontinuities in the solution; and (4) it removes constraints on the geometry and mesh, making it applicable to arbitrary geometries and spatial sampling (zero-shot super-resolution capability). We critically assess these features and analyze the network’s ability to generalize to unseen samples across both stationary and transient PDEs. The overall performance of the proposed method is promising, demonstrating its applicability to a range of challenging problems in computational mechanics.

Keywords Parameterized PDEs  $\cdot$ Physics-informed operator learning  $\cdot$ Conditional neural field  $\cdot$ PDE encoding  $\cdot$ Second-order meta-learning

1 Introduction

Numerical Solvers. Numerical methods for solving nonlinear PDEs have been central to scientific computing, enabling the simulation of complex physical phenomena. Classical approaches such as the Finite Element Method (FEM), Finite Difference Method (FDM), and Spectral Methods are widely used due to their robustness and versatility. FEM excels in handling complex geometries and boundary conditions, while FDM is known for its simplicity and efficiency on structured grids. Spectral methods, on the other hand, offer high accuracy for problems with smooth solutions by leveraging global basis functions. Despite their success, these methods often require significant computational resources, especially for high-dimensional problems, as they need to be repeated for any new set of input parameters.

Physics-Informed Neural Networks. Physics-Informed Neural Networks (PINN) have emerged as a promising alternative, leveraging neural networks to approximate solutions to PDEs by embedding the governing equations directly into the loss function Raissi et al. (2019). PINNs are particularly attractive for solving problems with sparse or incomplete data. In specific cases where all governing equations and boundary conditions are known, PINNs can be employed as forward solvers to compute the solutions. Here, the term solver refers to finding solutions that satisfy the governing equations, translating the problem into a constrained optimization task. However, this approach faces two significant challenges: (1) the training time for PINNs is still not competitive with traditional numerical solvers, except in very high-dimensional settings where classical methods struggle Rezaei et al. (2022); Grossmann et al. (2024), and (2) even at their best (considering recent advances), standard PINNs perform comparably to methods like FEM, and the solutions lack generalizability to other parametric inputs. Consequently, the core limitations of traditional numerical methods, such as computational inefficiency and lack of flexibility for parametric variations, persist when using PINNs solely for forward problems. The performance of PINNs in inverse problems, model discovery, and calibration differs from their application in forward problems, with several reported advantages in these contexts Faroughi et al. (2024).

It’s worth mentioning that, as a current trend, researchers are also combining ideas from FEM and neural networks to develop new deep learning-based solvers. See, for example, the following studies:Mitusch et al. (2021); Škardová et al. (2024); Xiong et al. (2025); Li et al. (2024). Although these approaches offer appealing features, particularly in terms of implementation simplicity, their computational time remains comparable to the FEM, and their applicability to highly nonlinear multiphysics problems has yet to be demonstrated.

Neural Operators. Neural Operators (NOs) extend the capabilities of traditional machine learning by learning mappings between infinite dimensional function spaces, making them powerful tools for solving parametric PDEs and modeling complex physical systems. NOs ideally should generalize across entire families of functions, enabling them to predict solutions for varying input conditions, such as material properties or boundary conditions. Rather well-known approaches for operator learning so far are Deep Neural Operator (DeepOnet) and its extensions Lu et al. (2021); Abueidda et al. (2025); He et al. (2024); Kumar et al. (2025); Yu et al. (2024), Physics-informed DeepOnet Wang et al. (2021); Mandl et al. (2025); Li et al. (2025), Fourier Neural Operator (FNO) and its extensions Li et al. (2021); Azizzadenesheli et al. (2024); Li et al. (2023a, b), Unet Ronneberger et al. (2015); Chen et al. (2019); Mendizabal et al. (2020); Mianroodi et al. (2022); Gupta et al. (2023); Najafi Koopas et al. (2025) and others.

The mentioned methods are also continuously improving their generalizability to complex geometries and their ability to handle complex solution fields. For example, Yin et al. (2024); Li et al. (2023a); Xiao et al. (2024); Zeng et al. (2025); Chen et al. (2024a) introduce various frameworks designed to learn geometry-dependent solution operators for different PDEs efficiently. Similar to the previous section, the combination of physics-informed NO and classical numerical solvers such as FEM, FDM, and FFT are also getting more attention, see for example Yamazaki et al. (2025a),Rezaei et al. (2025a), Xu et al. (2024), Eshaghi et al. (2025), Lee et al. (2025),Kaewnuratchadasorn et al. (2024), Franco et al. (2023) and Harandi et al. (2025). The latter helps avoid time-consuming data generation or reduces the required data, and delivers higher accuracy for unseen predictions.

Operator learning faces several key challenges that hinder its broader adoption in scientific computing.

• A major issue is generalization, as models often struggle to extrapolate beyond training data, especially in highly nonlinear systems with complex geometries.

• Capturing high-frequency features (discontinuous solution fields) remains challenging. In many industrial applications, one expects discontinuities in the solution fields which produce high gradients and therefore very relevant for design purposes (e.g. consider stress or flux values coming from derivation of the primary field such as displacement, or temperature).

• Data efficiency is a concern, as training requires large datasets from expensive simulations or experimental measurements. On that note, many operator-learning frameworks also lack physical constraints, leading to solutions that may violate fundamental laws.

• Scalability to large 3D problems is computationally expensive, and meta-learning techniques for adaptability are still underexplored. This also includes extending their applicability to irregular geometries or domains with complex topologies.

Implicit Neural Representations Implicit Neural Representations (INRs), also known as neural fields, represent a novel approach for encoding data as continuous functions parameterized by neural networks. Unlike traditional data storage methods, such as grids, meshes, or point clouds, INRs implicitly map input coordinates (e.g., spatial or temporal points) to corresponding output values (e.g., deformation, temperature, or color). This method enables efficient and flexible modeling of high-dimensional and complex data. Consequently, the application of INRs in scientific computing has been receiving increasing attention.

The literature and research on neural fields for scientific machine learning, particularly in computational mechanics, have gained momentum in recent years. Serrano et al. (2023) introduced a novel method that uses coordinate-based networks to solve PDEs on general geometries, removing input mesh constraints and excelling in diverse problem domains such as spatio-temporal forecasting and geometric design. Naour et al. (2024) proposed a continuous-time modeling approach for time series forecasting, leveraging implicit neural representations and meta-learning. Boudec et al. (2024)introduced a physics-informed iterative solver that learns to condition gradient descent for solving parametric PDEs, and improving optimization stability, accelerating convergence. See also contributions by Yeom et al. (2025) for achieving speed-up in training neural fields via weight scaling of sinusoidal neural fields. Hagnberger et al. (2024) proposed vectorized conditional neural fields to solve time-dependent PDEs, that have efficient inference, zero-shot super-resolution, and generalization to unseen PDE parameters through attention-based neural fields. Du et al. (2024a) discussed the conditional neural field latent diffusion model which is a generative deep learning framework that enables efficient and probabilistic simulation of complex spatiotemporal turbulence in 3D irregular domains, leveraging conditional neural fields and latent diffusion. Dupont et al. (2022a) proposed a neural compression framework that handles diverse data modalities by converting data into implicit neural representations. Their approach reduces encoding time by two orders of magnitude. Catalani et al. (2024) developed a methodology using INRs to learn surrogate models for steady-state fluid dynamics on unstructured domains, handling 3D geometric variations and generalizing to unseen shapes.

Despite these strengths, INRs face several challenges. Training INRs can be computationally expensive, requiring large amounts of data and careful regularization to prevent overfitting. Generalization remains a significant challenge, although meta-learning and transfer-learning approaches are emerging to address this limitation. Additionally, while INRs excel in representing continuous fields, their scalability to large-scale or high-dimensional data remains constrained by the capacity of the neural network and the computational cost of training. Representing disconnected or complex topologies without artifacts can also be challenging. Addressing these challenges will be critical for further advancing the utility and adoption of INRs in real-world applications.

Summary, open questions and our contributions. Traditional numerical solvers are highly optimized for solving challenging PDE problems. However, they are typically designed for one-time use, requiring costly recomputation whenever input parameters change. Deep learning methods, particularly neural operators, offer a promising alternative but are often highly data-demanding. Integrating physical constraints into neural operators has shown potential in mitigating these challenges, yet training such complex networks remains difficult, especially across a wide range of input parameters. Furthermore, the application of neural operators to real-world problems is still underdeveloped, with key challenges including scalability to complex geometries, accuracy in capturing sharp gradients, and efficient training, all of which remain active areas of research.

Conditional neural fields represent a promising yet straightforward network architecture that compactly and continuously maps function spaces across arbitrary domains. They offer significant flexibility for operator learning in parametrized PDEs, as their outputs can be sampled at arbitrary resolutions and modified by adjusting a set of latent variables. In contrast to discrete networks, where memory usage increases inefficiently with spatio-temporal resolution, neural fields require memory that scales primarily with the number of parameters in the network. Additionally, specialized architectures, such as sinusoidal activations or Fourier features, enhance their ability to capture high-frequency details, while integration with numerical methods extends their applicability to complex geometries and topologies. This work aims to address some of the above challenges by integrating concepts from the standard finite element method and conditional neural fields into operator learning. As a result, we introduce a physics-informed encoder-decoder network to establish the mapping between continuous parameter and solution spaces for parametrized PDEs. To the best of our knowledge, this is the first application of a single-network conditional neural field for operator learning, owing this capability to its unique loss formulation, which eliminates the need for the complex encode-process-decode pipeline (Serrano et al., 2023; Catalani et al., 2024). This approach is well-suited for complex geometries, offering data efficiency, requiring no labeled data, and demonstrating scalability to real-world industry applications, especially in the fields of materials engineering and science.

In Section 2, we summarize the main problem formulation and the types of PDEs and geometries chosen for this study. Next, in Section 3, we introduce implicit Finite Operator Learning as a promising approach for near-real-world problems. In Section 4, we present the study results, followed by conclusions and future directions in Section 5.

2 Parameterized PDEs

Let $u:\mathcal{T}\times\Omega\times\mathcal{C}\to\mathbb{R}^{d}$ be the solution field representing physical quantities such as temperature, displacement, velocity, or pressure. The solution field $u$ depends on the time $t\in\mathcal{T}\subset\mathbb{R}$, the spatial coordinates ${x}\in\Omega\subset\mathbb{R}^{d}$, and the control parameter $c\in\mathcal{C}\subset\mathbb{R}^{d^{\prime}}$. The solution field $u$ is parameterized by $c$ while satisfying the following PDE:

$$\mathcal{R}=\partial_{t}u(t,x;c)-f(x,u;c)=0, (t,x) \in\mathcal{T}\times\Omega,$$ (1)

$$u(0,x;c)-u_{0}(x;c)=0, (t,x) \in\{(0,x)\mid x\in\Omega\}.$$

Here, $u_{0}$ denotes the initial condition, and $f$ can be either linear or nonlinear with respect to both the solution and the control parameter, and it typically involves partial derivatives of $u$ with respect to $x$.

The control parameter $c$ can influence the solution field in various ways, such as defining the initial and boundary conditions, material properties, or other system characteristics. In more complex scenarios, $c$ may represent a spatially distributed field, accounting for variations in geometry or domain (e.g. material) heterogeneity.

Standard numerical methods, such as the finite element and finite volume methods, are widely used for solving Eq. 1. In practical applications and design processes, they are repeatedly executed for each variation in the control parameter $c$. These methods yield a numerical approximation $u(\cdot,\cdot;c):\mathcal{T}\times\Omega\to\mathbb{R}$ for varying $c\in\mathcal{C}$.

After applying a spatial discretization, such as the finite element method, we approximate the solution field as:

$$u(t,x;c)\approx u^{h}(t,x;c)=\sum_{i=1}^{M}u_{i}(t,c)\psi_{i}(x),$$ (2)

where $\{\psi_{i}(x)\}_{i=1}^{M}$ are spatial basis functions, $u_{i}(t,c)$ are the time-dependent, parameterized coefficients representing the discrete solution at spatial degrees of freedom, and $M$ is the number of grid points. Substituting this approximation into the governing PDE and applying the Galerkin projection, we obtain the semi-discrete system:

$$\mathbf{r}(\mathbf{u},\mathbf{c})=\mathbf{M}\frac{d\mathbf{u}}{dt}+\mathbf{K}(\mathbf{u},\mathbf{c})\mathbf{u}-\mathbf{g}(\mathbf{u},\mathbf{c})=\mathbf{0},$$ (3)

where:

• $\mathbf{u},\mathbf{r}\in\mathbb{R}^{M}$ are the vectors of unknowns (discrete solutions at spatial nodes) and nodal residuals, respectively,

• $\mathbf{M}\in\mathbb{R}^{M\times M}$ is the mass matrix, given by

  $$M_{ij}=\int_{\Omega}\psi_{i}(x)\psi_{j}(x)\,d\Omega,$$ (4)

• $\mathbf{K}(\mathbf{u},\mathbf{c})\in\mathbb{R}^{M\times M}$ represents the nonlinear stiffness, often resulting from the spatial derivatives of $u$. It may be dependent on both the solution as well as the problem parameters.

• $\mathbf{g}(\mathbf{u},\mathbf{c})$ is the discrete source term,

• $\mathbf{c}$ is the control vector that parametrizes the discrete system by governing any or a combination of initial and boundary conditions, material properties, domain geometry, and spatial heterogeneity.

This semi-discrete system is an ordinary differential equation (ODE) system in time, which must be further discretized using a time-stepping method. In this work, the implicit Euler method is used to approximate the time-dependent term in Eq. 3 with a chosen time step size $\Delta t$. Table 1 outlines the parametrized boundary value problems explored in this research for operator learning. These problems are selected from a wide range of applications, particularly in computational mechanics. The upper part of the table focuses on stationary problems, where the given PDE has no transient term or time evolution, while the lower part highlights two well-known nonlinear transient equations. For clarity, the FEM-based residual vector for each problem is specified in Table 3.

Table 1: Summary of the parametrized PDEs and their setup for operator learning.

2D,3D Steady-state problems	Mechanical problem hyper/linear elasticity	Nonlinear diffusion
PDE	$$\nabla\cdot\bm{P}+\bm{f}=0+\text{BCs},~{}\bm{P}=\frac{\partial W}{\partial\bm{F}}$$ $$W=\frac{\mu(\bm{X})}{2}(\bar{I_{1}}-3)-\frac{\kappa(\bm{X})}{4}(J_{F}^{2}-1-2\ln J_{F})$$ $$\nabla\cdot\bm{\sigma}+\bm{f}=0+\text{BCs},~{}\hat{\bm{\sigma}}=\bm{C}\hat{\bm{\varepsilon}}$$	$$-\nabla\cdot\bm{q}=0$$ $$\bm{q}=-k(\bm{x},T)~{}\nabla T$$ $$k(\bm{x},T)=k_{0}(\bm{x})(1+2T^{4})$$
Operator $\mathcal{O}$: In $\to$ Out	$\mathcal{O}:\bm{C}(\bm{x})\to\bm{U}(\bm{x})$,   $\mathcal{O}:\bm{U}_{b}\to\bm{U}(\bm{x})$	$\mathcal{O}:{k}(\bm{x})\to{T}(\bm{x})$
Geometry and BCs
Input Samples	$$\bm{C}(\bm{x})=\text{Fourier-based param.}$$ $$\bm{U}_{b}=\text{Random Gaussian}$$	$\displaystyle k_{0}=0.95\cdot\text{Sigmoid}\left(20\left(k_{f}-0.5\right)\right)$ $\displaystyle+0.05.$ $$k_{f}=\text{Fourier-based param.}$$
2D Transient problems	Non-linear thermal diffusion	Allen-Cahn equation
PDE	$$\rho c_{p}\frac{\partial T}{\partial t}=-\nabla\cdot\bm{q}+Q$$ $$\bm{q}=-K~{}\nabla T,~{}K=k_{0}(\bm{x})(1+\alpha T)$$	$$\frac{\partial\phi}{\partial t}=-M\left(\frac{\delta F}{\delta\phi}\right)$$ $$F=\int_{\Omega}\left(\frac{\epsilon}{2}|\nabla\phi|^{2}+f(\phi)\right)d\Omega$$
Operator $\mathcal{O}$, In $\to$ Out	$\mathcal{O}:T_{i}(\bm{x})\to T_{i+1}(\bm{x})$	$\mathcal{O}:\phi_{i}(\bm{x})\to\phi_{i+1}(\bm{x})$
Geometry and BCs
Input Samples	$$T_{0}\sim\mathcal{N}\left(0,\exp\left(-\tfrac{\|X_{i}-X_{j}\|^{2}}{2\varepsilon^{2}}\right)\right)$$	$$\phi_{0}\sim\mathcal{N}\left(0,\exp\left(-\tfrac{\|X_{i}-X_{j}\|^{2}}{2\varepsilon^{2}}\right)\right)$$

3 iFOL: implicit Finite Operator Learning

3.1 Implicit Neural Representations

Implicit Neural Representations (INRs) are multi-layer perceptron (MLP) networks that are coordinate-based and parameterized by $L$ layers of weights $\mathbf{W}_{i}$, biases $\mathbf{b}_{i}$, and nonlinear activation functions $\sigma_{i}$, with $\theta=(\mathbf{W}_{i},\mathbf{b}_{i})_{i=1}^{L}$. These networks model spatial fields as an implicit function that maps spatial coordinates to scalar or vector quantities $x\in\mathbb{R}^{d}\mapsto f_{\theta}(x)$. We adopt SIREN (Sitzmann et al., 2020) as the core INR architecture in our framework. This network utilizes sine activation functions, combined with a distinct initialization strategy.

$$\mathbf{SIREN}({x})=\mathbf{W}_{L}(\sigma_{L-1}\circ\sigma_{L-2}\circ\cdots\circ\sigma_{0}({x}))+\mathbf{b}_{L},\quad\text{with}\quad\sigma_{i}(\bm{\eta}_{i})=\sin\left(\omega_{0}(\mathbf{W}_{i}\bm{\eta}_{i}+\mathbf{b}_{i})\right)$$ (5)

Here, $\bm{\eta}_{0}=\bm{x}$, and $(\bm{\eta}_{i})_{i\geq 1}$ represent the hidden activations at each layer of the network. The parameter $\omega_{0}\in\mathbb{R}_{+}^{*}$ is a hyperparameter that governs the frequency bandwidth of the network. SIREN requires a specialized initialization of weights to ensure that outputs across layers adhere to a standard normal distribution.

3.2 Conditional Neural Fields

Our goal is to obtain a plausible solution to a parameterized partial differential equation using neural fields. This is achieved by conditioning the neural field on a set of latent variables $\bm{l}$, which can encode the solution field across arbitrary parameterizations and discretization of the underlying PDEs. By varying these latent variables, we can effectively modulate the neural solution field. $\bm{l}$ is typically a low-dimensional vector, and is also referred to as a feature code. For a comprehensive discussion and review of conditional neural fields, interested readers are referred to Xie et al. (2022).

iFOL utilizes Feature-wise Linear Modulation (FiLM Perez et al. (2018)), which conditions neural field in an auto-decoding manner. It employs a simple neural network without hidden layers (i.e., a linear transformation) to predict a shift vector from the latent variables $\bm{l}$ for each layer of the neural field network. This yields the shift-modulated SIREN:

$$u_{\theta,\gamma}({x},\bm{l}) =\text{Decode}(\bm{l})=\mathbf{W}_{L}\left(\sigma_{L-1}\circ\sigma_{L-2}\circ\cdots\circ\sigma_{0}({x})\right)+\mathbf{b}_{L},$$ (6)

$$\sigma_{i}(\bm{\eta}_{i},\bm{\phi}_{i}) =\sin\left(\omega_{0}(\mathbf{W}_{i}\bm{\eta}_{i}+\mathbf{b}_{i}+\bm{\phi}_{i})\right),$$

$$\bm{\phi}_{i}(\bm{l}) =\mathbf{V}_{i}\bm{l}+\mathbf{c}_{i}.$$

where $u_{\theta,\gamma}$ is the neural solution field, $\theta=(\mathbf{W}_{i},\mathbf{b}_{i})_{i=1}^{L}$ and $\gamma=(\mathbf{V}_{i},\mathbf{c}_{i})_{i=1}^{L-1}$ represent the trainable parameters of the SIREN network, referred to as the Synthesizer, and the FiLM hypernetworks, referred to as the Modulator, respectively.

3.3 PDE Loss Function

Building on physics-informed neural networks, we introduce a domain-integrated physical loss function whose variation with respect to the solution field yields the residuals of the governing partial differential equations. Among the functionals commonly used in computational mechanics, the total potential energy functional and the weighted residual functional naturally fulfill this property. Since the energy functional is physics-specific, not always explicitly known, and its derivation can be challenging, we adopt the well-established weighted residual functional as follows:

$$\mathcal{L}_{\mathbf{PDE}}=\int_{\Omega}u_{\theta,\gamma}\mathcal{R}\,d\Omega=0.$$ (7)

This formulation weights the residual using the neural field as a test function, enforcing its vanishing in an integral sense. Applying the chain rule and variational calculus, the gradient of the loss function with respect to the predicted solution is given by:

$$\delta_{u_{\theta,\gamma}}\mathcal{L}_{\mathbf{PDE}} =\int_{\Omega}\mathcal{R}\,d\Omega+\int_{\Omega}u_{\theta,\gamma}\,\delta_{u_{\theta,\gamma}}\mathcal{R}\,d\Omega$$ (8)

$$=\int_{\Omega}\mathcal{R}\,d\Omega$$

where $\delta_{u_{\theta,\gamma}}\mathcal{R}$ is zero due to the stationarity of the residual functional with respect to the solution field. To compute the loss function efficiently, we discretize the domain and employ the corresponding discrete residuals (Eq. 3) as follows:

$$\mathcal{L}_{\mathbf{PDE}}(\mathbf{u}_{\theta,\gamma}(t,\bm{x},\bm{l}),\mathbf{c})=\sum_{e=1}^{n_{el}}({\mathbf{u}_{\theta,\gamma}^{e}})^{T}\mathbf{r}^{e},~{}\text{where}~{}~{}\mathbf{r}^{e}=\mathbf{M}^{e}\frac{d\mathbf{u}_{\theta,\gamma}^{e}}{dt}+\mathbf{K}^{e}\mathbf{u}_{\theta,\gamma}^{e}-\mathbf{g}^{e}.$$ (9)

where $\mathbf{u}_{\theta,\gamma}$ represents the neural solution field vector evaluated at the mesh points, and the superscript $e$ indicates that the quantity is evaluated at the element level. Here note that, although the PDE loss presented here is based on the residual, one can also directly utilize the energy functional of a PDE as shown in Table 3 in Appendix A.4 as well. In fact, the energy functional-based loss in Table 3 was employed for the transient problems.

This strategy enables the seamless integration of well-established numerical methods into the learning process while minimizing computational overhead. Consequently, it allows for the construction of PINNs without relying on resource-intensive automatic differentiation to evaluate the components of the PDE.

3.4 Training

The training and inference of the conditional neural fields (CNFs) involve the computation of the latent variables $\bm{l}$, a step commonly referred to as encoding. In data-driven operator learning methodologies utilizing CNFs, encoding is performed on both the input and output fields, following the so-called Encode-Process-Decode framework(Dupont et al., 2022b; Yin et al., 2022; Serrano et al., 2023; Catalani et al., 2024). In contrast, this work uniquely encodes the PDE and the underlying physics, rather than the spatial fields. In each training step of iFOL, the latent codes for the sample batch $\mathcal{B}$ are first derived by minimizing the physical loss with respect to the latent codes in just a few steps of gradient descent:

$$\bm{l}^{*}(\mathbf{c}_{i}) =\text{Encode}(\mathbf{PDE})$$ (10)

$$=\arg\min_{\bm{l}}\mathcal{L}_{\mathbf{PDE}}=\sum_{e=1}^{n_{el}}({\mathbf{u}_{\theta,\gamma}^{e}})^{T}\mathbf{r}^{e}(\mathbf{u}^{e}_{\theta,\gamma},\mathbf{c}^{e}_{i}),\quad\forall i\in\mathcal{B}$$

Subsequently, the parameters of the Synthesizer and Modulator networks are optimized using the computed latent codes:

$$\theta^{*},\gamma^{*}=\arg\min_{\theta,\gamma}\mathcal{L}_{\mathbf{PDE}}=\sum_{e=1}^{n_{el}}({\mathbf{u}_{\theta,\gamma}^{e}})^{T}\mathbf{r}^{e}(\mathbf{u}^{e}_{\theta,\gamma},\mathbf{c}^{e}_{i}),\quad\forall i\in\mathcal{B}$$ (11)

Basically, we partition the model into context-specific parameters, which dynamically adapt to individual samples, and meta-trained parameters, which are globally optimized to enable knowledge transfer across diverse contexts. To the best of our knowledge, this work presents the first application of second-order meta-learning, inspired by the CAVIA algorithm (Zintgraf et al., 2019), for the parametric learning of PDEs in a physics-informed manner. The details of the approach are provided in Algorithms 1 and 2. Here, $\alpha$ denotes the encoding learning rate, while $\lambda$ is the training learning rate, which adjusts the weights of the modulator and synthesizer networks.

Algorithm 1 Encode PDE for sample batch $\mathcal{B}$

1:Initialize: Set codes to zero $\bm{l}_{i}\leftarrow 0,\quad\forall i\in\mathcal{B}$;

2:for each $i\in\mathcal{B}$ and step $\in\{1,\dots,K_{e}\}$ do

3:     $\bm{l}_{i}\leftarrow\bm{l}_{i}-\alpha\nabla_{\bm{l}_{i}}\mathcal{L}_{\mathbf{PDE}}(\mathbf{u}_{\theta,\gamma}(t,\bm{x},\bm{l}_{i}),\mathbf{c}_{i})$

4:end for

Algorithm 2 Training of iFOL

1:while no convergence do

2:     for each mini-batch $\mathcal{M}\subseteq\mathcal{B}$ do

3:         /* Encode PDE */

4:         for each $i\in\mathcal{M}$ and step $\in\{1,\dots,K_{e}\}$ do

5:              $\bm{l}_{i}\leftarrow\bm{l}_{i}-\alpha\nabla_{\bm{l}_{i}}\mathcal{L}_{\mathbf{PDE}}(\mathbf{u}_{\theta,\gamma}(t,\bm{x},\bm{l}_{i}),\mathbf{c}_{i})$

6:         end for

7:         /* Update Modulator & Synthesizer */

8:         $\theta,\gamma\leftarrow\theta,\gamma-\lambda\frac{1}{|\mathcal{M}|}\sum_{i\in\mathcal{M}}\nabla_{\theta,\gamma}\mathcal{L}_{\mathbf{PDE}}(\mathbf{u}_{\theta,\gamma}(t,\bm{x},\bm{l}^{*}_{i}),\mathbf{c}_{i})$;

9:     end for

10:end while

To provide further clarity on the architecture presented above, a comparative analysis between the proposed approach and two well-established operator learning algorithms is shown in Fig. 1. Furthermore, a more detailed description of each component of the iFOL framework is provided in Fig. 2. The upper section of the figure outlines the approach for stationary problems, while the lower section demonstrates the application of the same architecture to transient problems. Notably, once trained on transient problems, the network can be repeatedly invoked to generate solutions over time, without the need for additional retraining.

Figure 1: Comparison of different architectures for operator learning. The brownish arrows indicate the flow of gradients in each method, while the input parametric space and output continuous field are highlighted with reddish and bluish circles, respectively. Notably, in the proposed iFOL method, the input space first influences the loss term, after which the latent code is constructed.

Figure 2: Training and inference procedures in the iFOL framework. Top: Details of the iFOL architecture for quasi-static problems, where the goal is to predict the corresponding solution for a given input parameter space in a single step. Bottom: Details of the iFOL architecture for transient problems, where the only difference lies in the inference step—the trained network is repeatedly called to predict the solution field over time.

4 Results

4.1 Stationary problems: hyper-elastic mechanical equilibrium

We begin by evaluating the performance of iFOL in predicting the mechanical deformation of a hyperelastic heterogeneous solid when the property distribution varies, i.e., the operator maps $\mathcal{O}:\bm{C}(\bm{x})\to\bm{u}(\bm{x})$. The strong form and energy formulation of the solid in this context follows a nonlinear form, which is detailed in Tables 1 and 3 and Section A.1, along with the selected material properties.

Figure 3: Results for the 2D stationary mechanical equilibrium PDE considering a hyperelastic material model. The operator learns to map the elasticity distribution to the deformation fields (see also Table 1 and Section A.1).

The $8000$ random training samples are generated using a Fourier-based parametrization with $f_{x}=\{2,4,6\}$ and $f_{y}=\{2,4,6\}$ (See also Rezaei et al. (2025b) for details on this parameterization scheme). The chosen network parameters are listed in Table 4.

In Fig. 3, the results of iFOL are shown for four unseen test cases, all evaluated at more than twice the training resolution (i.e., approximately six times more grid points). The first row presents predictions for unseen higher-frequency components as well as low-frequency ones with unsymmetrical property distribution. In the second row, we further challenge the network by testing its performance on polycrystalline microstructures, which go well beyond the training samples used in this study. Notably, in the last case featuring a multiphase polycrystal, we even alter the material property values (e.g., Young’s modulus) to include ranges not encountered during training.

Despite significant changes in the input space and resolution, the network produces reasonable predictions. Across all test cases, the errors remain well below one order of magnitude compared to the applied displacement and maximum deformation in the solid.

At this point, a valid question is how we can ensure the quantity and quality of the initial training samples to perform a certain task. Although a solid and unique answer to this question requires much further intensive study, as it is influenced by many parameters, we attempt to address it here by altering the number of training samples while keeping the network hyperparameters the same. As expected and shown in Fig. 4, increasing the number of samples systematically reduces the errors up to a certain level. Interestingly, the entire framework appears to be very sample- or data-efficient, as extensive samples are not required to achieve reasonable results. The same pattern repeats in other case studies, but we omit them for the sake of brevity.

Figure 4: Studies on the influence of the number of training samples on the obtained errors for unseen test samples. Increasing the number of samples reduces the prediction errors while increasing the training time.

4.2 Stationary problems: elasticity and learning on boundary conditions

In this study, we examine the performance of iFOL in learning the solution for given BCs while keeping material properties and geometry fixed. More details on the chosen boundary values are provided in Section A.1, while details on the formulation of the loss term can be found in Tables 1 and 3. The selected network hyperparameters are listed in Table 4. Here, we define the operator as $\mathcal{O}:\bm{U}_{b}\to\bm{U}(\bm{x})$. We focus on gyroid surfaces, a specific class of metamaterials. The complex topology of the gyroid enables novel functionalities, which highlights its importance in the design of next-generation materials. Moreover, we consider applications in multiscale analysis, where the applied boundary conditions typically originate from the macroscale, and it is crucial to determine the microstructural response to the imposed macroscopic displacement field.

For this purpose, we generate 1000 random samples for the applied displacement vector on the front surface of the chosen metamaterial and train iFOL to learn the full deformation field for a given Dirichlet BC. In Fig. 5, the results of iFOL versus FEM are shown for three different test cases. For visualization purposes, the deformations are magnified by a factor of 10. The maximum pointwise errors for unseen test cases are one order of magnitude lower than the applied displacement field and are mainly concentrated in specific regions where Neumann BCs are applied. It is worth noting that, within the iFOL framework, Dirichlet BCs are enforced in a hard manner, which explains the zero error at the front and back surfaces.

Figure 5: Results for the 3D stationary mechanical equilibrium PDE considering a linear elastic material model. The operator learns to map the given applied Drichlet boundary conditions to the deformation fields (see also Table 1 and Section A.1).

4.3 Stationary problems: Nonlinear diffusion

In this study, we not only critically assess the applicability of iFOL in approximating nonlinear operators on complex geometries with irregular meshes but also evaluate its effectiveness for sensitivity analysis. We are specifically interested in the derivatives of domain-integrated functionals that explicitly depend on the solution field (the operator’s output) with respect to the parameter field (the operator’s input), e.g.,

$$\mathcal{J}=\int_{\Omega}T(k)\,d\Omega,~{}~{}~{}\delta_{k}\mathcal{J}=\int_{\Omega}\delta_{k}T\,d\Omega.$$ (12)

Here, $\delta_{k}T$ represents the Jacobian of the operator $\mathcal{O}:{k}({\bm{x}})\to{T}({\bm{x}})$ and can be efficiently computed by applying automatic differentiation (AD) to the trained iFOL network as a postprocessing step.

Figure 6: Results for stationary diffusion PDE parameterized by a temperature-dependent and spatially varying conductivity field (see also Table 1 and Section A.2). The training mesh is 7× coarser than the evaluation mesh.

For this example, 8000 samples were used for training, and 2000 samples were reserved for testing as unseen data. The conductivity field is parametrized using the Fourier parametrization as introduced in (Rezaei et al., 2024), with $f_{x}=\{1,2,3\}$, $f_{y}=\{1,2,3\}$, and $f_{z}=\{0\}$. The so-called sigmoidal projection was applied to generate dual-phase conductivity fields within the range of 0.5 to 1.0. For further details, refer to the corresponding column in Table 1. A discussion and details on the network architecture are provided in Section B and Table 4. For the summary of the mathematical formulation see also Section A.2.

Figure 6 showcases how iFOL performs in predicting temperature fields for unseen conductivity distributions compared to the reference finite element method. The complicated 3D geometry of FOL is discretized using a fully unstructured tetrahedral mesh. The input parameter exhibits a complex spatial distribution, while the temperature field is governed by a thermal diffusion PDE with high nonlinearity arising from the temperature-dependent conductivity field. iFOL and reference FE solutions are presented for two samples that were not included in the training set. We observe that the learned operator exhibits generalization both to the parameter field and to the evaluation grid, which is 7× finer than the training mesh.

Figure 7 displays the sensitivity maps of $\delta_{k}\mathcal{J}$ for two unseen samples, where one map is obtained via automatic differentiation applied to iFOL’s inference function, and the other is computed using adjoint-based sensitivity analysis. Qualitative comparisons show a high degree of agreement between the two maps, with both displaying similar spatial sensitivity patterns. To the authors’ knowledge, this is the first time that the Jacobian of neural network-based operators for spatially distributed input and output fields has been evaluated and compared against classical sensitivity analysis techniques, such as adjoint methods. Based on the results presented here, we can conclude that iFOL is not only able to accurately approximate PDE-based operators in a mesh- and geometry-agnostic manner, but also to reliably estimate the operator’s Jacobian.

Figure 7: Sensitivity maps of $\delta_{k}\mathcal{J}$ for two unseen samples, comparing AD-iFOL and adjoint-based analysis.

4.4 Transient problems: nonlinear thermal diffusion in heterogeneous domain

As mentioned before, for transient problems, we adopt a slightly different strategy, where we train the network to learn the dynamics of the transient problem from step $i$ to $i+1$. Upon successful training, the trained network can then be called multiple times in a loop to predict the evolution of the desired solution field. In this section, we present the results for nonlinear thermal diffusion in a heterogeneous domain where we have a map between $T_{i}(x,y)$ and $T_{i+1}(x,y)$. The heterogeneity is intentionally introduced to create sharp solution jumps, which typically correspond to higher frequencies in the solution and are usually not easily captured by a naive network setup. From a physical perspective, this can be related to different phases in a material microstructural system. The additional nonlinearity arises from the fact that the material conductivity is not only a function of space but also varies with temperature (i.e., the solution field). See also Section A.2 for the summary of the mathematical formulation and further details on the material properties. The results of this study are summarized in Fig. 8.

Figure 8: Results for the transient nonlinear diffusion PDE considering a temperature-dependent and heterogeneous conductivity. The operator learns to map the initial conditions to the next evolved temperature field (see also Table 1 and Section A.2).

Here, the training fields are randomly generated using the Gaussian process with varying length scales of the radial basis function kernel, and the results of the temperature profile evolution for different unseen initial conditions are provided. At this point, the resolution for both training and testing is the same (a grid of $51\times 51$). Note that any obtained temperature profile as an output serves as a new input for the next step. Therefore, for transient problems, errors can accumulate over time, which has also been reported in previous studies Du et al. (2024b); Dingreville et al. (2024); Yamazaki et al. (2025b). The network’s capability to handle unseen temperature fields is also worth mentioning. For example, from the fifth time step onward, the temperature profile is primarily dominated by BCs and does not resemble the randomly generated fields process used to train the network. Yet, the predictions remain reasonably accurate. We emphasize that the predictions of up to 50 steps shown in this figure are purely based on the one-time-trained network, with no labeled data or hybrid solver used to further support the network.

As a final remark, we should point out our fundamentally different approach to handling transient problems, which involves breaking the time axis as explained. In the classical PINN approach, the time variable is yet another input in addition to space. Although this approach initially makes sense and appears natural, training with respect to both time and space can become very challenging for high-dimensional problems, as the network can simply violate the so-called causality, requiring additional enhancements that are not needed here Wang et al. (2024).

4.5 Transient problems: Allen-Cahn equation on the complex domain

In this section, we present the results for the nonlinear Allen-Cahn equation in a homogeneous irregular domain, where we establish a mapping between $\phi_{i}(x,y)$ and $\phi_{i+1}(x,y)$. See also Section A.3 for the summary of the mathematical formulation. The irregular domain is intentionally chosen to demonstrate the approach’s capability to handle complex geometries. From a practical perspective, such shapes are commonly observed in active material particles on the cathode side of battery systems, where phase transitions in the material govern ion diffusion Chen et al. (2024b). As before, training samples are generated from the Gaussian process to train the network. The results in Fig. 9 illustrate two different initializations, leading to distinct final phase profiles. The same strategy for the nonlinear thermal problem applies here, so we do not repeat every detail. Predictions remain accurate up to 10 time steps, after which errors gradually accumulate, a point that will be analyzed further. It is important to note that the resolution for both training and testing is kept the same, as it was found to be sufficient based on the length scale parameter in the phase field equation. Further discussion on super-resolution capabilities is provided in Section 4.6. We did not leverage any precomputed solution fields or hybrid solver coupling. However, such techniques could be incorporated in future developments to further reduce errors and maintain accuracy within expected margins. See investigations by Li et al. (2023c); Dingreville et al. (2024).

4.6 Studies on zero-shot super resolution

Zero-shot super-resolution (ZSSR) in operator learning refers to the ability of a trained model to enhance the resolution of function mappings without requiring additional high-resolution training data. Unlike traditional super-resolution techniques that rely on supervised learning with paired low- and high-resolution samples, zero-shot approaches leverage the inherent structure of the learned operator to generalize across different resolutions Li et al. (2021); Koopas et al. (2024).

ZSSR is particularly useful in scientific computing, where performing calculations on very fine discretizations is computationally expensive. It is worth noting the similarities to model order reduction techniques, where the dimensionality of the problem is reduced by capturing the dominant solution modes Quarteroni and Rozza (2011).

The ZSSR capability of iFOL is also demonstrated in Figures 3 and 6. Here, we further validate this claim by systematically testing the proposed iFOL approach and quantitatively assessing the error accumulation across different mesh resolutions. Our main focus is on stationary nonlinear diffusion in 3D, where we evaluate the solution on both finer and coarser grids, as well as on the 2D nonlinear Allen-Cahn equation, where we analyze error accumulation over time and for varying resolutions. Similar behavior is observed for other investigated problems, and therefore, for brevity, we do not report all details.

In Fig.10, we evaluate the network on two additional mesh resolutions. Since we use unstructured meshes and nodes in the 3D domain, we aim to assess whether the network can generate results for both coarser and finer meshes. For error measurement in this particular example, we adopt a strict criterion based on the maximum pointwise error. To ensure the results are not case-dependent, we systematically test multiple samples, with five representative cases shown for each resolution in Fig.10. The averaged values are highlighted in red.

Figure 9: Results for the transient Allen-Cahn PDE considering in a complex-shape domain. The operator learns to map the initial conditions to the next evolved phase-field variable (see also Table 1 and Section A.3).

As expected, the prediction error increases when moving to resolutions different from the training resolution. However, as demonstrated in Fig. 6, the overall results remain acceptable. A potential way to mitigate this issue is to train the network on multiple resolutions across different samples, reducing bias toward a specific mesh resolution.

Figure 10: Performance of iFOL in ZSSR for nonlinear thermal diffusion in a heterogeneous 3D problem with unstructured meshes. The maximum pointwise error is lowest for the trained resolution and increases for other mesh resolutions, yet remains within an acceptable range.

Figure 11: Performance of iFOL in ZSSR for the nonlinear transient Allen-Cahn equation in an arbitrary 2D domain with unstructured meshes. The relative L2 error is lowest for the trained resolution but increases for other mesh resolutions and accumulates over time during further network evaluations.

In Fig. 11, we analyze the transient Allen-Cahn problem and compare results using relative L2 error measurements across three different resolutions. On the far left, we show the trained resolution, followed by two additional cases with two and four times more nodes, respectively. As expected, and similar to other physical problems, the error increases slightly when moving to different resolutions, particularly finer ones. Moreover, this diagram illustrates how errors accumulate over time for a given resolution (observe the spacing between circles, squares, and triangles along the y-axis). A potential remedy is to train the model using multiple resolutions.

Figure 12: Prediction of iFOL for different time steps and resolutions for the nonlinear Allen-Cahn equation.

4.7 Computational cost analysis

Accuracy and computational cost are two key criteria for comparing classical techniques, such as the FEM and adjoint-based sensitivity analysis, with machine learning-based approaches, including PINNs and operator learning methods. In this section, we analyze the computational cost of the iFOL technique using the 3D nonlinear diffusion problem described in Section 4.3. Leveraging Google JAX Bradbury et al. (2018), we developed a unified framework that seamlessly integrates classical solver operations—such as element-wise computations, global Jacobian assembly, and right-hand side (RHS) construction—with machine learning techniques, including loss function evaluation and gradient-based optimization. This framework is implemented within the FOL paradigm Rezaei et al. (2024), utilizing key JAX features such as jax.vmap for efficient vectorized operations and jax.jit for just-in-time (JIT) compilation, enabling the translation of computationally demanding tasks into optimized machine code. The loss function is implemented exactly as formulated in Eq. 9. Furthermore, the physical loss function not only returns the element-wise weighted residual or its energy counterpart as a scalar loss function to be minimized during the learning process but also provides the corresponding elemental Jacobian and RHS. This unified implementation enables the same framework to be used for both FEM solvers and parameterized learning approaches such as iFOL.

Table 2 summarizes the computational times for the primal and sensitivity analyses performed using iFOL and FEM. or the sake of brevity, we focus on the example reported in Section 4.3, where we deal with a nonlinear PDE in a 3D complex geometry. A similar trend of results was also observed for other cases. We note that training iFOL on 8000 samples with 3800 iterations required 14 hours on NVIDIA A100 GPU with 40 GB of RAM (see Table 4). However, during inference, iFOL achieves speedups of 200×, 1000×, and 8300× compared to FEM on coarse (1069 nodes), normal (5485 nodes), and fine (37362 nodes) meshes, respectively. Beyond the speedup in the primal inference, iFOL also achieved 1.3× and 13× speedups for sensitivity analysis on normal and fine meshes, respectively. It should be noted that the reported times for AD-iFOL and Adjoint-FEM correspond exclusively to sensitivity analysis, assuming the primal solution is already available. Furthermore, the computational cost of iFOL-based sensitivity analysis remains independent of the number of response functions. In contrast, adjoint-based sensitivity analysis scales directly with the number of responses, as a separate adjoint system of equations must be solved for each response.

Table 2: Computational costs of the studies in the 3D nonlinear diffusion problem Sec. 4.3.

	Training (h)	Batch size	Inference (ms/sample)			# Iterations/sample
Mesh res. (# nodes)			1069	5485	37362
iFOL	14	64	1.26	1.60	3.59	3
FEM	-	64	256	1599	29961	10
AD-iFOL	-	32	169	173	228	0
Adjoint-FEM	-	32	175	233	3152	1

5 Conclusion

We introduced the implicit Finite Operator Learning (iFOL) technique, a novel framework for solving parametric PDEs on arbitrary geometries. iFOL is a physics-informed automatic encoder-decoder that combines three state-of-the-art techniques: conditional neural fields for spatial encoding, second-order meta-learning for PDE encoding, and a physics-informed loss function inspired by the original FOL concept. This operator learning technique is capable of providing both a precise continuous mapping between parameter and solution spaces and a reliable estimation of the operator’s Jacobian. The computed Jacobian can be directly utilized for sensitivity analysis and gradient-based optimization without incurring additional computational costs.

The applicability of the methodology is tested on several classes of PDEs, with special attention to problems in computational mechanics and materials engineering, namely hyperelasticity, linear elasticity, nonlinear thermal diffusion, and the Allen-Cahn equation. We address stationary and transient problems in 2D and 3D, including both regular and complex geometries. The variety of examples is intended to cover most potential applications. In all cases, iFOL is able to parametrically learn the solution and predict unseen test samples, even beyond the training regime, demonstrating its generalization capability to new scenarios without any retraining.

It is worth highlighting that no labeled data is used throughout this work. Instead, the PDEs in the weighted residual form or in their energy counterparts are directly defined as loss functions. This approach improves efficiency, as no automatic differentiation operator is required to construct the loss terms. Moreover, it allows for the direct integration of residuals from other numerical solvers, such as FEM or spectral methods like FFT. Finally, since the network learns the equation parametrically for a class of problems, the inference time of iFOL can be highly competitive with direct numerical solvers like FEM. The speed-up factor is not only problem-dependent but also resolution-dependent and influenced by implementation details, making it difficult to provide a universal comparison, especially considering differences in hardware setups for DL-based solvers and classical numerical methods. In specific cases and under comparable conditions, we observe at least a 100$\times$ speed-up for nonlinear problems, which can further increase for higher resolutions.

Outlook and future research direction Physics-informed neural operators, and more specifically the iFOL framework, offer significant potential for the next generation of numerical solvers by leveraging their ability to learn solutions parametrically. A natural next step is to apply the introduced methodologies to more complex multiphysics problems, where multiple physical processes influence each other through one- or two-way coupling. In such scenarios, classical solvers often struggle and require staggered approaches. Moreover, in this work, we primarily rely on a single input parameter space at a time. A natural extension would be to incorporate multiple input parameter spaces, a direction that has been explored to some extent by other authors, primarily using data-driven operator learning techniques. Finally, since we have access to meaningful and reliable sensitivity information independently of the number of objectives, gradient-based optimization becomes a compelling approach. This capability has significant potential for efficiently solving complex optimization and inverse problems, even in time-dependent settings, while maintaining reasonable computational and implementation costs.

Code Availability

The implementation of iFOL and the examples presented in this study will be made publicly available under the FOL framework repository on GitHub at https://github.com/RezaNajian/FOL, following the publication of this work.

Authors contributions

R.N.A.: Conceptualization, Methodology, Software, Writing - Review & Editing. Y.Y.: Software, Writing - Review & Editing. K.T.: Software. M.M.: Funding, Review & Editing. M.A.: Funding, Review & Editing. S.R.: Supervision, Methodology, Software, Writing - Review & Editing.

Appendix A Review on physical formulations and their discretization using FEM

A.1 Mechanical problem

Here, we summarize the mechanical problem where the position of material points is denoted by $\bm{X}^{T}=[x,~{}y,~{}z]$. We denote the displacement components by $U$, $V$, and $W$ in the $x$, $y$, and $z$ directions, respectively. The kinematic relation defines the strain tensor $\bm{\varepsilon}$ in terms of the deformation vector $\bm{U}^{T}=[U,~{}V,~{}W]$ and reads:

$$\bm{\varepsilon}\, =\,\text{sym}(\text{grad}(\bm{U}))=\nabla^{s}\bm{U}=\dfrac{1}{2}\left(\nabla\bm{U}+\nabla\bm{U}^{T}\right),$$ (13)

$$\bm{F} =\dfrac{\partial\bm{x}}{\partial\bm{X}},~{}\bm{C}=\bm{F}^{T}\bm{F},~{}J_{F}=\det(\bm{F}).$$ (14)

Here, $\bm{F}$ is the deformation gradient, and $J_{F}$ is the Jacobian determinant, representing the local volume change. The elastic energy of the solid for linear and nonlinear (hyperelastic Neo-Hookean material) case reads as

$$\psi_{lin} =\frac{\lambda}{2}(\text{tr}\bm{\varepsilon})^{2}+\mu\,\text{tr}(\bm{\varepsilon}^{2}),$$ (15)

$$\psi_{nonlin} =\frac{\mu(\bm{X})}{2}(\bar{I_{1}}-3)-\frac{\kappa(\bm{X})}{4}(J_{F}^{2}-1-2\ln J_{F}).$$ (16)

The Lamé constants are denoted by $\lambda$ and $\mu$, and $\text{tr}(\bm{\varepsilon})$ is the trace of the strain tensor, representing the volumetric strain. $\bar{I}_{1}=\text{tr}(\mathbf{\bar{C}})$ is the first invariant of the isochoric right Cauchy-Green deformation tensor $\mathbf{\bar{C}}=\mathbf{J}^{-2/3}\mathbf{C}$ where $\mathbf{C}$ is the right Cauchy-Green tensor. We then define the Cauchy stress tensor $\bm{\sigma}=\frac{\partial\psi_{lin}}{\partial\mathbf{\epsilon}}={\mathbb{C}}~{}{\bm{\varepsilon}}$ and the Piola-Kirchhoff stress $\bm{P}=\frac{\partial\psi_{nonlin}}{\partial\mathbf{F}}$, which are obtained by differentiating the strain energy function. Here $\mathbb{C}(\lambda,\mu)=\lambda~{}\mathbf{I}\otimes\mathbf{I}+2\mu~{}\mathbb{I}^{s}$ or $\mathbb{C}(E,\nu)$ is the fourth-order elasticity tensor by defining $\mathbf{I}$ as the second-order identity tensor and $\mathbb{I}^{s}$ as the symmetric fourth-order identity tensor. Finally, the mechanical equilibrium in the absence of body force, as well as the Dirichlet and Neumann boundary conditions, are written as:

$$\text{div}({\bm{\sigma}})=\text{div}(\mathbb{C}\nabla^{s}\bm{U})\ =\bm{0}~{}~{}~{}~{}\text{in}~{}~{}~{}\Omega$$ (17)

$$\bm{U} =\bar{\bm{U}}~{}~{}~{}\text{on}~{}~{}\Gamma_{D}$$ (18)

$$\bm{\sigma}\cdot\bm{n}=\bm{t} =\bar{\bm{t}}~{}~{}~{}~{}\text{on}~{}~{}\Gamma_{N}$$ (19)

In the above relations, $\Omega$ and $\Gamma$ denote the material points in the body and on the boundary area, respectively. Moreover, the Dirichlet and Neumann boundary conditions are introduced in Eq. 18 and Eq. 19, respectively. In the case of hyperelasticity, we then have $\text{Div}({\bm{P}})=\bm{0}$ instead of Eq. 17.

A.2 Thermal problem

The heat equation describes how the temperature $T(\bm{x},t)$, with $\bm{x}$ representing position and $t$ the time, evolves within the domain $\bm{x}\in\Omega$ over time. Let the heat source be $Q:\Omega\times(0,\tau)\to\mathbb{R}$, the boundary temperature $T_{D}(\bm{x}):\Gamma_{D}\times(0,\tau)\to\mathbb{R}$, and the boundary heat flux be $q_{N}:\Gamma_{N}\times(0,\tau)\to\mathbb{R}$, where $\Gamma_{D}$ and $\Gamma_{N}$ are the regions where the Dirichlet and Neumann boundary conditions are applied, respectively. Here, $t\in(0,\tau)$ represents the temporal domain, with $\tau$ denoting the end time. The strong form of the heat equation is given as:

$$c\rho\frac{\partial T(\bm{x},t)}{\partial t} =-\text{div}(\bm{q})+Q\quad\text{in}\ \Omega\times(0,\tau),$$ (20)

$$T(\bm{x},t) =T_{D}(\bm{x})\quad\text{on}\ \Gamma_{D}\times(0,\tau),$$ (21)

$$\nabla T(\bm{x},t)\cdot\bm{n} =q_{N}(\bm{x})\quad\text{on}\ \Gamma_{N}\times(0,\tau),$$ (22)

$$T(\bm{x},0) =T_{0}(\bm{x})\quad\text{in}\ \Omega.$$ (23)

In the above equations, $\bm{n}$ is the outward normal vector, $c$ is the specific heat capacity, and $\rho$ is the density. The heat flux is given by $\bm{q}=-K\nabla T$, where $K=k_{0}(\bm{x})(1+\alpha T)$ represents the thermal conductivity, which depends on both position and temperature.

The above set of equations forms the basis for the results and studies presented in Section 4.4. For the stationary thermal problem (see Section 4.3), the time derivative in Eq. 20 naturally vanishes.

A.3 Allen-Cahn equation

The Allen-Cahn equation describes the evolution of an order parameter $\phi(\bm{x},t)$, where $\bm{x}$ represents position and $t$ is time, within the domain $\bm{x}\in\Omega$. Let the source term be $S:\Omega\times(0,\tau)\to\mathbb{R}$, the boundary value be $\phi_{D}(\bm{x}):\Gamma_{D}\times(0,\tau)\to\mathbb{R}$, and the boundary flux be $q_{N}:\Gamma_{N}\times(0,\tau)\to\mathbb{R}$, where $\Gamma_{D}$ and $\Gamma_{N}$ are the regions where the Dirichlet and Neumann boundary conditions are applied, respectively. Here, $t\in(0,\tau)$ represents the temporal domain, with $\tau$ denoting the end time. The strong form of the Allen-Cahn equation is given as:

$$\frac{\partial\phi(\bm{x},t)}{\partial t} =-M\left(\frac{\delta F}{\delta\phi}\right)\quad\quad\quad\quad\quad\quad\text{in}\ \Omega\times(0,\tau),$$ (24)

$$F =\int_{\Omega}\left(\frac{\epsilon}{2}|\nabla\phi|^{2}+f(\phi)\right)d\Omega~{}~{}~{}~{}\text{in}~{}~{}~{}\Omega\times(0,\tau),$$ (25)

$$\phi(\bm{x},t) =\phi_{D}(\bm{x})~{}~{}\quad\quad\quad\quad\quad\quad\quad\quad\text{on}\ \Gamma_{D}\times(0,\tau),$$ (26)

$$\nabla\phi(\bm{x},t)\cdot\bm{n} =\phi_{N}(\bm{x})~{}~{}\quad\quad\quad\quad\quad\quad\quad\quad\text{on}\ \Gamma_{N}\times(0,\tau),$$ (27)

$$\phi(\bm{x},0) =\phi_{0}(\bm{x})~{}~{}~{}\quad\quad\quad\quad\quad\quad\quad\quad\text{in}\ \Omega.$$ (28)

In the above equations, $\bm{n}$ is the outward normal vector, $M$ is the mobility coefficient, and $\frac{\delta F}{\delta\phi}$ represents the functional derivative of the free energy, typically given by:

$$\frac{\delta F}{\delta\phi}=-\epsilon^{2}\nabla^{2}\phi+f^{\prime}(\phi),$$ (29)

where $\epsilon$ is a small parameter controlling the interface thickness and $f^{\prime}(\phi)=(\phi^{2}-1)\phi$ denotes the derivative of the double-well potential enforcing phase separation. The above set of equations forms the basis for the results and studies presented in Section 4.5. The Allen-Cahn equation is a fundamental phase-field model and plays a crucial role in modeling microstructure evolution in alloys Allen and Cahn (1979); Karma and Rappel (1998). The equation captures the dynamics of interfaces and transitions between phases (including applications in phase-field fracture Bourdin et al. (2000); Rezaei et al. (2021)). Its applications extend to image processing, tumor growth modeling, and topology optimization.

A.4 A short review on FEM formulation

Next, we briefly summarize the FEM discretization techniques, particularly for a 2D quadrilateral element, to clarify the loss terms in the iFOL formulation, especially for implementation purposes. The details provided here are fairly standard, and the extension to other types of complex element formulations with different polynomial orders should be straightforward. Readers are also encouraged to see the standard procedure in any finite element subroutine Hughes (2000); Bathe (1996). The corresponding linear shape functions $\bm{N}$ and the deformation matrix $\bm{B}$ used to discretize the mechanical weak form in the current work.

$$\bm{N}=\begin{bmatrix}N_{1}&0&\dots&N_{4}&0\\ 0&N_{1}&\dots&0&N_{4}\end{bmatrix},~{}\bm{B}=\begin{bmatrix}N_{1,x}&0&\dots&N_{4,x}&0\\ 0&N_{1,y}&\dots&0&N_{4,y}\\ N_{1,y}&N_{1,x}&\dots&N_{4,y}&N_{4,x}\end{bmatrix}.$$ (30)

The notation $N_{i,x}$ and $N_{i,y}$ represent the derivatives of the shape function $N_{i}$ with respect to the coordinates $x$ and $y$, respectively. To compute these derivatives, we utilize the Jacobian matrix

$$\bm{J}=\partial\bm{X}/\partial\bm{\xi}=\begin{bmatrix}\frac{\partial x}{\partial\xi}&\frac{\partial y}{\partial\xi}\\ \frac{\partial x}{\partial\eta}&\frac{\partial y}{\partial\eta}\end{bmatrix}.$$ (31)

Here, $\bm{X}=[x,y]$ and $\bm{\xi}=[\xi,\eta]$ represent the physical and parent coordinate systems, respectively. It is worth mentioning that this determinant remains constant for parallelogram-shaped elements, eliminating the need to evaluate this term at each integration point. Finally, for the B matrix, we have

$$\bm{B}=\bm{J}^{-1}\begin{bmatrix}\frac{\partial N_{1}}{\partial\xi}&0&\frac{\partial N_{2}}{\partial\xi}&0&\frac{\partial N_{3}}{\partial\xi}&0&\frac{\partial N_{4}}{\partial\xi}&0\\ 0&\frac{\partial N_{1}}{\partial\eta}&0&\frac{\partial N_{2}}{\partial\eta}&0&\frac{\partial N_{3}}{\partial\eta}&0&\frac{\partial N_{4}}{\partial\eta}\\ \frac{\partial N_{1}}{\partial\eta}&\frac{\partial N_{1}}{\partial\xi}&\frac{\partial N_{2}}{\partial\eta}&\frac{\partial N_{2}}{\partial\xi}&\frac{\partial N_{3}}{\partial\eta}&\frac{\partial N_{3}}{\partial\xi}&\frac{\partial N_{4}}{\partial\eta}&\frac{\partial N_{4}}{\partial\xi}\end{bmatrix}.$$ (32)

For each element the deformation field and stress tensor $\hat{\bm{\sigma}}$ are approximated as

$$U=\sum{N}_{i}u_{i}=\bm{N}\bm{U}_{e},~{}V=\sum{N}_{i}V_{i}=\bm{N}\bm{V}_{e},~{}\hat{\bm{\sigma}} =\bm{C}\bm{B}\bm{U}_{e},$$ (33)

Here, $\bm{U}^{T}_{e}=\left[\begin{matrix}u_{1}&\cdots&u_{4}\\ \end{matrix}\right]$ and $\bm{V}_{e}$ are the nodal values of the deformation field of element $e$ in the $x$ and $y$ directions, respectively. The same procedure applies to both the thermal problem and the Allen-Cahn equation. However, one must account for additional time derivatives. We summarize the implemented loss terms in Table 3.

Table 3: Summary of the FEM-based residuals of the investigated problems.

PDE type	FEM residual / iFOL loss function
Stationary nonlinear mechanical $\bm{C}(x,y)\to\bm{U}(x,y)$ $\mathbf{U}_{\theta,\gamma}\equiv\bm{U}_{e}$, (Sec. 4.1)	$r_{e}=-\int_{\Omega_{e}}\bm{S}\colon\delta\bm{E}dV-\int_{\Omega_{e}}(\Delta\bm{S}\colon\delta\bm{E}+\bm{S}\colon\Delta\delta\bm{E})dV+\int_{\Omega_{t}}\delta\bm{U}.~{}\bm{T}dA$
	$\mathcal{L}_{iFOL}=\sum_{e=1}^{n_{el}}\sum_{k=1}^{n_{int}}~{}\left(\frac{\mu(\bm{\xi}_{k})}{2}~{}\left(\text{tr}(\bar{\bm{C}}(\bm{U_{e}}))-3\right)~{}+~{}\frac{\kappa(\bm{\xi}_{k})}{4}~{}(J_{F}^{2}(\bm{U_{e}})-2\ln J_{F}(\bm{U_{e}})-1)\right)~{}\text{det}(\bm{J})~{}w_{k}$
Stationary linear mechanical $\bm{U}_{b}\to\bm{U}(\bm{x})$ $\mathbf{U}_{\theta,\gamma}\equiv\bm{U}_{e}$, (Sec. 4.2)	$r_{e}=-\int_{\Omega_{e}}[\bm{B}]^{T}\bm{C}~{}[\bm{B}]\bm{U}_{e}~{}dV+\int_{\Omega_{t}}[\bm{N}]^{T}\bm{C}~{}[\bm{B}\bm{U}_{e}]^{T}~{}\bm{n}~{}dS$
	$\mathcal{L}_{iFOL}=\sum_{e=1}^{n_{el}}\bm{U}^{T}_{e}\left[\bm{\left missing}(\sum_{k=1}^{n_{int}}\dfrac{w_{k}}{2}~{}\text{det}(\bm{J})~{}[\bm{B}(\bm{\xi}_{k})]^{T}\bm{C}_{e}(\bm{\xi}_{k})\bm{B}(\bm{\xi}_{k})\right)~{}\bm{U}_{e}\right].$
Stationary nonlinear thermal $K(x,y)\to T(x,y)$ $\mathbf{U}_{\theta,\gamma}\equiv T_{e}$, (Sec. 4.3)	$r_{e}=-\int_{\Omega_{e}}[\bm{B}]^{T}\bm{K}(\bm{T}_{e})~{}[\bm{B}]\bm{T}_{e}~{}dV+\int_{\Omega_{t}}[\bm{N}]^{T}\bm{t}~{}dS+\int_{\Omega_{e}}[\bm{N}]^{T}Q~{}dV$
	$\mathcal{L}_{iFOL}=\sum_{e=1}^{n_{el}}\left[\bm{T}_{e}^{n+1}\right]^{T}\left[\bm{\left missing}(\sum_{k=1}^{n_{int}}\dfrac{w_{k}}{2}~{}\text{det}(\bm{J})~{}[\bm{B}(\bm{\xi}_{k})]^{T}\bm{K}(\bm{T}_{e}^{n+1},\bm{\xi}_{k})\bm{B}(\bm{\xi}_{k})\right)~{}\bm{T}_{e}^{n+1}\right]$
Transient nonlinear thermal $T^{n}(x,y)\to T^{n+1}(x,y)$ $\mathbf{U}_{\theta,\gamma}\equiv T_{e}^{n+1}$, (Sec. 4.4)	$r_{e}=-\int_{\Omega_{e}}[\bm{B}]^{T}\bm{K}(\bm{T}_{e}^{n+1})~{}[\bm{B}]\bm{T}_{e}^{n+1}~{}dV+\int_{\Omega_{t}}[\bm{N}]^{T}\bm{t}~{}dS+\int_{\Omega_{e}}[\bm{N}]^{T}Q~{}dV$            $-\int_{\Omega_{e}}[\bm{N}]^{T}\rho c_{p}~{}[\bm{N}]\frac{\bm{T}_{e}^{n+1}-\bm{T}_{e}^{n}}{\Delta t}~{}dV$
	$\mathcal{L}_{iFOL}=\sum_{e=1}^{n_{el}}\left[\bm{T}_{e}^{n+1}\right]^{T}\left[\bm{\left missing}(\sum_{k=1}^{n_{int}}\dfrac{w_{k}}{2}~{}\text{det}(\bm{J})~{}[\bm{B}(\bm{\xi}_{k})]^{T}\bm{K}(\bm{T}_{e}^{n+1},\bm{\xi}_{k})\bm{B}(\bm{\xi}_{k})\right)~{}\bm{T}_{e}^{n+1}\right]$                 + $\sum_{e=1}^{n_{el}}\left[\bm{\sum}_{k=1}^{n_{int}}\dfrac{w_{k}}{2}~{}\text{det}(\bm{J})~{}\rho c_{p}\left[\bm{N}(\bm{\xi}_{k})\right]^{T}\frac{\left(\bm{N}(\bm{\xi}_{k})\bm{T}_{e}^{n+1}-\bm{N}(\bm{\xi}_{k})\bm{T}_{e}^{n}\right)^{2}}{\Delta t}\right].$
Transient Allen-Cahn $\phi^{n}(x,y)\to\phi^{n+1}(x,y)$ $\mathbf{U}_{\theta,\gamma}\equiv\phi_{e}^{n+1}$, (Sec. 4.5)	$r_{e}=-\int_{\Omega_{e}}[\bm{B}]^{T}\bm{[}\bm{B}]\bm{\phi}_{e}^{n+1}~{}dV-\int_{\Omega_{e}}[\bm{N}]^{T}~{}[\bm{N}]\frac{\bm{\phi}_{e}^{n+1}-\bm{\phi}_{e}^{n}}{\Delta t}~{}dV$            $-\int_{\Omega_{e}}[\bm{N}]^{T}\frac{1}{\varepsilon^{2}}\left[\left(\bm{N}\bm{\phi}_{e}^{n+1}\right)^{2}-1\right]\bm{N}\bm{\phi}_{e}^{n+1}~{}dV+\int_{\Omega_{t}}[\bm{N}]^{T}\left(\nabla\bm{\phi}_{e}^{n+1}\cdot\bm{n}\right)~{}dS.$
	$\mathcal{L}_{iFOL}=\sum_{e=1}^{n_{el}}\left[\bm{\phi}_{e}^{n+1}\right]^{T}\left[\bm{\left missing}(\sum_{k=1}^{n_{int}}\dfrac{w_{k}}{2}~{}\text{det}(\bm{J})~{}[\bm{B}(\bm{\xi}_{k})]^{T}\bm{B}(\bm{\xi}_{k})\right)~{}\bm{\phi}_{e}^{n+1}\right]$                $+\sum_{e=1}^{n_{el}}\left[\bm{\sum}_{k=1}^{n_{int}}\dfrac{w_{k}}{2}~{}\text{det}(\bm{J})~{}[\bm{N}(\bm{\xi}_{k})]^{T}\frac{\left(\bm{N}(\bm{\xi}_{k})\bm{\phi}_{e}^{n+1}-\bm{N}(\bm{\xi}_{k})\bm{\phi}_{e}^{n}\right)^{2}}{\Delta t}\right]$                $+\sum_{e=1}^{n_{el}}\left[\bm{\sum}_{k=1}^{n_{int}}w_{k}~{}\text{det}(\bm{J})~{}\frac{1}{\varepsilon^{2}}\dfrac{\left[(\bm{N}(\bm{\xi}_{k})\bm{\phi}_{e}^{n+1})^{2}-1\right]^{2}}{4}\right].$

In all the reported studies, we utilized simple first-order shape functions. The element type used in the studies varies based on the application, ranging from structured to unstructured meshes, including quadrilateral and tetrahedral elements, demonstrating the flexibility of the approach in handling problems with complex geometries.

Moreover, in the above set of equations, $n_{\text{int}}$ and $n_{\text{el}}$ represent the number of Gaussian integration points and number of elements, respectively. $\bm{\xi}_{k}$ and $w_{k}$ denote the coordinates and weight of the $k$-th integration point. The determinant of the Jacobian matrix is denoted by $\text{det}(\bm{J})$. Readers are encouraged to refer to Bathe (1996); Hughes (2000) as well as Rezaei et al. (2025b); Yamazaki et al. (2025b) and the references therein for more details on the FE formulation.

Appendix B Choice of Hyperparameters

A comprehensive summary of the network’s hyperparameters and configurations for all the models introduced above, as well as the various analyses conducted, is provided in Tables 4 and 5. NVIDIA Quadro RTX 6000s with 24 GB of RAM and NVIDIA A100 with 40 GB of RAM were utilized for the training of the stationary problems. The studies on the transient problems were performed on NVIDIA GeForce RTX 4090 with 24 GB of RAM.

We conducted extensive studies on multiple parameters, including latent size, learning rate, neural network structure, and the frequency of sinusoidal functions, among others. The values reported in Tables 4 and 5 correspond to the best-performing configurations. In general, we did not observe significant improvements with frequency parameter $\omega_{0}$ other than $30$, which interestingly performed well across various problem classes. Increasing the depth and latent size proved consistently beneficial with the sinusoidal activation function used in this study, though it comes with a higher training cost. Depending on problem complexity, reducing the number of epochs and latent iterations can further accelerate both training and inference in the proposed iFOL method. As discussed in the results section, having a larger number of samples also leads to better performance. On the other hand, the training cost may increase depending on the chosen batch size, and one must adapt the network architecture and latent size accordingly.

Table 4: Network hyperparameters for stationary problems

Training parameter	Nonlin. elas. (Sec. 4.1)	Lin. elas. (Sec. 4.2)	Staion. Nonlin. therm. (Sec. 4.3)
Number of samples	8000	1000	8000
Grid in training	$41\times 41$	$5605$	$5485$
Grid in evaluation	$101\times 101$	$5605$	37362
Synthesizers	[64]*4	[256]*6	[128]*6
$\omega_{0}$	30	30	30
Modulators	Linear (FiLM)	Linear (FiLM)	Linear (FiLM)
Latent size	512	1024	512
Number of latent iterations	3	3	3
Latent/encoding learning rate	$10^{-2}$	$10^{-2}$	$10^{-2}$
Training learning rate	$10^{-5}$	$10^{-4}$ to $10^{-7}$	$10^{-5}$
Batch size	320	120	100
Gradient normalization	Yes	Yes	Yes
Number of epochs	10000	10000	3800
Optimizer	Adam	Adam	Adam
Total trainable parameters	143938	330755	476417

Table 5: Network hyperparameters for the transient problems

Training parameter	Transient nonlin. thermal (Sec. 4.4)	Transient nonlin. Allen-Cahn (Sec. 4.5)
Number of samples	6000	8000
Grid in training	$51\times 51$	2624
Grid in evaluation	$51\times 51$, $101\times 101$	2624, 5678, 9873
Synthesizers	[256]*6	[256]*6
$\omega_{0}$	30	30
Modulators	Linear (FiLM)	Linear (FiLM)
Latent size	256	256
Number of latent iterations	3	3
Latent/encoding learning rate	$1.0\times 10^{-2}$	$1.0\times 10^{-2}$
Training learning rate	$1.0\times 10^{-4}$ to $1.0\times 10^{-7}$	$1.0\times 10^{-4}$ to $1.0\times 10^{-7}$
Batch size	120	100
Gradient normalization	Yes	Yes
Number of epochs	10000	10000
Optimizer	Adam	Adam
Total trainable parameters	723457	723201