I-FENN with DeepONets: accelerating simulations in coupled multiphysics problems

TL;DR

This paper introduces a hybrid FEM-neural network framework using DeepONets within I-FENN to efficiently simulate coupled multiphysics problems, achieving high accuracy and generalization with reduced computational costs.

Contribution

It presents a novel I-FENN architecture integrating DeepONets for multiphysics simulations, enhancing generalizability and efficiency over existing methods.

Findings

Achieves over 95% accuracy in unseen scenarios.

Reduces computational cost by simplifying FEM calculations.

Demonstrates scalability with model complexity.

Abstract

Coupled multiphysics simulations for high-dimensional, large-scale problems can be prohibitively expensive due to their computational demands. This article presents a novel framework integrating a deep operator network (DeepONet) with the Finite Element Method (FEM) to address coupled thermoelasticity and poroelasticity problems. This integration occurs within the context of I-FENN, a framework where neural networks are directly employed as PDE solvers within FEM, resulting in a hybrid staggered solver. In this setup, the mechanical field is computed using FEM, while the other coupled field is predicted using a neural network (NN). By decoupling multiphysics interactions, the hybrid framework reduces computational cost by simplifying calculations and reducing the FEM unknowns, while maintaining flexibility across unseen scenarios. The proposed work introduces a new I-FENN architecture with extended generalizability due to the DeepONets ability to efficiently address several combinations of natural boundary conditions and body loads. A modified DeepONet architecture is introduced to accommodate multiple inputs, along with a streamlined strategy for enforcing boundary conditions on distinct boundaries. We showcase the applicability and merits of the proposed work through numerical examples covering thermoelasticity and poroelasticity problems, demonstrating computational efficiency, accuracy, and generalization capabilities. In all examples, the test cases involve unseen loading conditions. The computational savings scale with the model complexity while preserving an accuracy of more than 95\% in the non-trivial regions of the domain.

Tables

Table 1. Table 1: Hyperparameters of the DeepONet model trained for the 3D cube example

	Input Size	$N_{GRU}$	$N_H$	$N_{ch}$	$N_{FC}$	$D_{out}$
Branch # 1	$N_l=$512	2	200	-	1	200
Branch # 2	$N_b=$512	2	50	25	1	50
Trunk	$N_d=$3	-	200	-	4	250

Table 2. Table 2: Comparison of FEM and I-FENN computation times in minutes across various domain sizes.

Domain size $[$Cells$]$	FEM time $[$Minutes$]$	I-FENN time $[$Minutes$]$	Savings
$30\times 30\times 30$	5.1	3.3	35%
$45\times 45\times 45$	20.0	12.5	38%
$60\times 60\times 60$	56.9	34.2	40%
$75\times 75\times 75$	134.0	75.8	43%

Table 3. Table 3: Hyperparameters of the DeepONet model trained for the 3D tube example

	Input Size	$N_{GRU}$	$N_H$	$N_{ch}$	$N_{FC}$	$D_{out}$
Branch # 1	$N_l=$128	2	64	-	1	64
Branch # 2	$N_b=$512	2	64	32	1	64
Trunk	$N_d=$3	-	256	-	4	128

Table 4. Table 4: Hyperparameters of the surrogate DeepONet model trained for simulation substitution for the 3D tube example

	Input Size	$N_{GRU}$	$N_H$	$N_{ch}$	$N_{FC}$	$D_{out}$
Branch # 1	$N_l=$128	2	84	-	1	84
Trunk	$N_d=$3	-	332	-	4	84

Table 5. Table 5: Hyperparameters of the DeepONet model trained for the 2D excavation example

	Input Size	$N_{GRU}$	$N_H$	$N_{ch}$	$N_{FC}$	$D_{out}$
Branch # 1	$N_l=$1	2	64	-	1	64
Branch # 2	$N_b=$584	2	64	32	1	64
Trunk	$N_d=$2	-	256	-	4	128

Equations

$${\int }_{\Omega }{\nabla }^s\hat{\mathbf{w}}:(\hat{\mathbf{C}}:{\mathbf{\epsilon }}_{n+1})d\Omega -{\int }_{\Omega }{\nabla }^s\hat{\mathbf{w}}:(\hat{\mathbf{C}}:\alpha (T_{n+1}-T_0)I)d\Omega$$

$$={\int }_{{\Gamma }_t}\hat{\mathbf{w}}\cdot {\overline{\mathbf{t}}}_{n+1}d\Gamma +{\int }_{\Omega }\hat{\mathbf{w}}\cdot {\mathbf{b}}_{n+1}d\Omega \forall \hat{\mathbf{w}}\in {\mathcal{W}}_u$$

$${\int }_{\Omega }\hat{T}\rho C_{\epsilon }\frac{T_{n+1}-T_n}{\Delta t}d\Omega +{\int }_{\Omega }\nabla \hat{T}\cdot \left(k\nabla T_{n+1}\right)d\Omega +{\int }_{\Omega }\hat{T}\alpha (n_{dim}\lambda +2\mu )\text{tr}\left(\frac{{\mathbf{\epsilon }}_{n+1}-{\mathbf{\epsilon }}_n}{\Delta t}\right)T_{0}d\Omega$$

$$=-{\int }_{{\Gamma }_q}\hat{T}{\overline{Q}}_{n+1}d\Gamma +{\int }_{\Omega }\hat{T}{r}_{n+1}d\Omega \forall \hat{T}\in {\mathcal{W}}_T$$

$${\int }_{\Omega }{\nabla }^s\hat{\mathbf{w}}:\hat{\mathbf{C}}:{\mathbf{\epsilon }}_{n+1}d\Omega -{\int }_{\Omega }{\nabla }^s\hat{\mathbf{w}}:({\alpha }^{\prime }{p}_{n+1}I)d\Omega$$

$$={\int }_{{\Gamma }_t}\hat{\mathbf{w}}\cdot {\overline{\mathbf{t}}}_{n+1}d\Gamma +{\int }_{\Omega }\hat{\mathbf{w}}\cdot {\mathbf{b}}_{n+1}d\Omega \forall \hat{\mathbf{w}}\in {\mathcal{W}}_u$$

$${\int }_{\Omega }\hat{p}\frac{p_{n+1}-p_n}{\Delta tM}d\Omega +{\int }_{\Omega }\hat{p}{\alpha }^{\prime }\text{tr}\left(\frac{{\mathbf{\epsilon }}_{n+1}-{\mathbf{\epsilon }}_n}{\Delta t}\right)d\Omega +{\int }_{\Omega }\nabla \hat{p}\cdot \left(\frac{K_{I}I}{{\mu }_f}\cdot \nabla p_{n+1}\right)d\Omega$$

$$=-{\int }_{\Omega }\nabla \hat{p}\cdot \left(\frac{K_{I}I}{{\mu }_f}\cdot {\gamma }_f{\mathbf{i}}_g\right)d\Omega -{\int }_{{\Gamma }_q}\hat{p}{{\overline{q}}_f}_{n+1}d\Gamma +{\int }_{\Omega }\hat{p}{Q_f}_{n+1}d\Omega \forall \hat{p}\in {\mathcal{W}}_p$$

$$\mathcal{G}(F_I)(x_O)=\sum _{k=1}^{D_{out}}{b}_k(F_I)t_k(x_O)+b_0$$

$$F_I:D_I\ni x_I$$

$$F_O:D_O\ni x_O$$

$$\mathcal{G}:{\mathcal{F}}_{\mathcal{I}}\ni F_I$$

$${\mathbf{r}}_t$$

$${\mathbf{z}}_t$$

$${\mathbf{h}}_t^{\prime }$$

$${\mathbf{h}}_t$$

$$\mathcal{G}(F_I)(x_O)$$

$$\ell (x_O)$$

$$\mathcal{G}(F_I)(x_O)$$

$${\ell }_i(x_O)$$

$${\varphi }^{\ast }=\arg \underset{\varphi }{\min }\mathcal{L}(\varphi )$$

$$e$$

$${\mathcal{L}}_{L_2}$$

$${\mathcal{L}}_{L_{2}Norm}$$

$${\mathcal{L}}_{\text{SSE}}$$

$${\mathcal{L}}_{\text{MSE}}$$

$$L_2^t$$

$$L_2^{LC}$$

$$L_2^{AllLCs}$$

$${ϵ}_{rel}$$

$$\nabla \cdot \mathbf{\sigma }+\mathbf{b}$$

$$\rho T_o\dot{\eta }+\nabla \cdot \mathbf{q}-r$$

$$\mathbf{u}$$

$$\mathbf{t}$$

$$T$$

$$Q$$

$$\mathbf{q}=-k\nabla T,$$

$$\mathbf{\sigma }$$

$$\rho T_o\dot{\eta }$$

$$\rho C_{\epsilon }\dot{T}+\alpha \mathbf{I}:\hat{\mathbf{C}}:\dot{\mathbf{\epsilon }}+\nabla \cdot \mathbf{q}-r$$

$${\int }_{\Omega }{\nabla }^s\hat{\mathbf{w}}:(\hat{\mathbf{C}}:\mathbf{\epsilon })d\Omega -{\int }_{\Omega }{\nabla }^s\hat{\mathbf{w}}:(\hat{\mathbf{C}}:\alpha (T-T_0)I)d\Omega$$

Full text

*** Fouad Amin, \linkablefma9357@nyu.edu ** ${\dagger}$ Diab W. Abueidda, \linkableda3205@nyu.edu

${\ddagger}$ Panos Pantidis, \linkablepp2624@nyu.edu **§ Mostafa E. Mobasher, \linkablemostafa.mobasher@nyu.edu **

1 Introduction

1.1 Literature review

A wide range of natural phenomena incorporates multiphysics interactions, often characterized by the coupling between the mechanical behavior and other physical fields. Thermomechanical coupling, for example, plays a critical role in diverse applications such as additive manufacturing [1], thermal protection systems [2], and metal solidification [3]. Hydromechanical coupling is essential in seepage and consolidation problems [4], fluid-structure interaction [5], and slope stability analysis [6]. More complex multiphysics processes exist in several domains, including infrastructure, energy, and manufacturing [7, 8, 9]. The large-scale and intricate character of these phenomena naturally limits the scope of experimental approaches for their investigation, rendering computational methods the primary pathway for studying the evolution of the phenomena. Modeling of these multiphysics systems is often based on conventional numerical approaches such as Finite Element Method (FEM), Finite Difference/Volume Method (FDM/FVM) [10], among others [11, 12, 13, 14]. These techniques are often associated with high computational costs, which can be prohibitively expensive for complex, high-dimensional, large-scale problems [15, 16]. To address these challenges, researchers have developed a range of numerical remedies, including domain decomposition, parallel computing, staggered solvers, and adaptive remeshing techniques [17, 18, 19, 20]. Advancing existing computational capabilities remains a vibrant and rapidly evolving area of research, with a profound impact on our understanding of multiphysics real-world problems.

Numerous machine learning (ML) approaches have been proposed to address the computational challenges associated with coupled mechanics problems. These approaches can be categorized into two main groups based on training strategies: (1) data-driven methods, and (2) physics-informed methods. Data-driven methods rely solely on supervised learning using labeled datasets to train the ML models. On the other hand, physics-informed methods incorporate physics into loss functions or manipulate the ML model output to obey physical laws [21]. Another categorization approach for ML adoption in computational mechanics can be based on the methodology: (1) simulation substitution, and (2) simulation enhancement [22]. Simulation substitution methods aim to replace the entire numerical simulation with an ML (surrogate) model to predict all response fields, whereas simulation enhancement focuses on improving a specific component of the simulation, including preprocessing, modeling, and postprocessing.

Neural networks (NNs) play a central role in the successful application of ML-based methods to computational mechanics problems. NNs are nonlinear approximators of functions, mapping input to output variables through a set of learnable parameters[23]. Several types of networks have been developed over the years, including multi-layer perceptrons (MLPs), sequence models (RNNs, GRUs, and transformers), and others [24, 25]. Recently, to further enhance the expressivity and generalization capabilities of ML techniques, operator-based architectures have emerged and are being actively developed. Lu et al. [26] proposed the deep operator network (DeepONet), showcasing its ability to learn both explicit operators (such as integrals) and implicit operators (such as differential equations). Alongside DeepONet, alternative deep operator networks have also been introduced [27, 28, 29]. A typical DeepONet architecture consists of two main components: (1) a branch network, and (2) a trunk network. The branch encodes the discretized input function at fixed scattered locations called ”sensors”, while the trunk encodes the domain of the output function described by a set of target locations that can be selected arbitrarily. The branch and trunk outputs are then aggregated through a dot product operation to produce the final discretized output at the target locations. Given its advanced generalization capabilities, DeepONet has been used in a wide range of applications, including large-scale carbon storage operations [30] and domain decomposition [31], using an innovative hybrid NN-FEM framework. In addition, several variations of DeepONet have been introduced, aiming to enhance its generalization capabilities and overcome training challenges. For instance, the multiple-input deep neural operators (MIONet) [32] extends DeepONet to handle multiple input functions. Another example is Fusion DeepONet[33], which implements data transfer between the branch and trunk hidden layers.

Despite the efficient computational performance of ML models, they suffer from significant limitations like limited interpretability [34, 35], generalizability [36], reliability [37, 38], sensitivity to hyperparameters [37, 39], network convergence [40, 41], and training stability [42, 43]. Therefore, the simulation substitution approach struggles in large-scale applications because it demands massive datasets and expensive training. In addition, without enforcing governing equations, time-history predictions can suffer from uncontrolled error propagation. In this context, it is evident that there is a need for more robust and reliable hybrid approaches that can combine the strengths of both numerical methods and ML models, while addressing their limitations. To this end, several hybrid formulations have been developed, such as the Integrated Finite Element Neural Network (I-FENN) proposed by the authors [44], as well as other approaches which integrate neural networks and operator networks with FEM for multiple applications such as FEMIN [45], hybrid FEM-NN [46, 47], NNFE [48], and hybrid PI-DeepONet with FEM [31]. I-FENN is a framework designed to accelerate and improve the accuracy of multiphysics simulations, and it is conceptually distinct from the frameworks mentioned above. The core novelty of I-FENN is to employ NNs as PDE approximators in conjunction with a conventional FEM solver, establishing a hybrid scheme that is conceptually similar to staggered approaches for coupled problems. This approach is different from the other hybrid approaches [45, 46, 47, 48, 31], where a NN is used to approximate solutions for a specific part of the mesh or a specific term in a balance law, or as a surrogate model trained using a differentiable FEM representation. Over the past years, the authors have demonstrated the applicability of I-FENN across several linear and non-linear problems, including phase-field fracture[49], non-local gradient damage[50], and thermoelasticity[51, 52], while investigating the error convergence analysis and model performance[37]. However, each of these efforts had its unique limitations and challenges, including the need for expensive training, data-driven dependencies, and others. One of the main issues was the poor generalization capabilities, with methods often restricted to a single snapshot in time or a specific load history.

1.2 Scope and Outline

In this study, we expand the scope and modeling capabilities of I-FENN, targeting two- and three-dimensional thermoelasticity and poroelasticity problems under varying (body and surface), unseen loading conditions. For the first time, we demonstrate the feasibility and efficiency of integrating a neural operator (DeepONet) as the network of choice within I-FENN, a development that significantly expands the framework’s generalization capability. We propose a modified MIONet architecture to accommodate multiple input functions. In addition, we introduce a simplified approach for enforcing boundary conditions across different boundary segments. The proposed architecture comprises a fully connected network (MLP) in the trunk and a gated recurrent unit (GRU) in the branch, with the latter capturing temporal dependencies without requiring a predetermined number of prior input time steps. [24]. The framework is integrating a finite element solver built upon the deal.II library [53]. Deal.II is a widely used open-source C++ library that supports creating finite element solvers with the capability of parallel processing. The framework implementation enables high-performance computing (HPC) deployment while maintaining efficient data transfer and memory management.

The proposed framework is applied to three key case studies: (1) a thermoelasticity problem with a spatially and temporally varying thermal body load, (2) a thermoelasticity problem with a spatially and temporally varying thermal surface load, and (3) a poroelasticity problem featuring a time-dependent fluid flux. The examples are provided for various geometric shapes spanning both two-dimensional (2D) and three-dimensional (3D) domains, demonstrating the proposed framework’s ability to predict over different geometries and loading conditions. Importantly, all testing examples are provided for load conditions that were never used during training, demonstrating the proposed framework’s ability to generalize to unseen scenarios. Additionally, examples of mesh refinement are provided to illustrate the framework’s scalability and independence from specific mesh resolutions.

The rest of the document is structured as follows: Section 2 presents a concise overview of the mathematical formulation for thermoelasticity and poroelasticity problems. Section 3 details the theoretical foundation of the I-FENN framework, including its mathematical formulation for thermoelasticity and poroelasticity problems. In addition, this section outlines the implementation aspects, including data management and transfer between the DeepONet model and the FEM solver. Section 4 delves into the DeepONet architecture developed in this work. Section 5 presents the numerical examples used to test the proposed framework. Finally, Section 6 summarizes this work and offers a future outlook.

2 Problem Statement

In this section, we introduce the mathematical formulation for thermoelasticity and poroelasticity problems. For the sake of brevity, here we introduce only the time-discretized weak forms, and a more detailed description of the strong and weak formulations is provided in Appendix A.

2.1 Thermoelasticity

For the thermoelasticity problem under consideration, the unknown variables are the displacement field $\boldsymbol{u}$ and the temperature profile $T$ across the domain, which are coupled through the governing system of equations shown in Eq. (A.1)-(A.2). The finite element analysis for the thermoelasticity problem is implemented using the discretized weak form, where the simulation is incremented from time step $n$ to time step $n+1$. Using an implicit Euler scheme, the incremental weak form of the balance of linear momentum equation reads as:

[TABLE]

where $\boldsymbol{\nabla}$ is the gradient operator and $\boldsymbol{\nabla}\cdot$ is the divergence operator, $T_{o}$ is the reference temperature  and $\boldsymbol{\varepsilon}$ is the strain defined as a function of displacement ($\boldsymbol{u}$) as $\boldsymbol{\varepsilon}=\boldsymbol{\nabla}^{s}\boldsymbol{u}=\frac{1}{2}(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{T})$ for small deformations, $\boldsymbol{\overline{t}}$ is the boundary traction, $\boldsymbol{b}$ is the mechanical body force vector, $\boldsymbol{\hat{C}}$ is the elasticity tensor, $\alpha$ is the coefficient of thermal expansion, $\Omega$ is the physical domain, $\Gamma$ is the domain boundary, $\widehat{\boldsymbol{w}}$ is the displacement test function, and $\mathcal{W}_{u}$ denotes the displacement function space.

The time-discretized weak form of the conservation of energy equation can be written as:

[TABLE]

where $\overline{Q}$ is the boundary heat flux, $r$ is the heat body source or sink, $\rho$ is the material mass density, $\mu$ and $\lambda$ are Lamé constants, $k$ is the thermal conductivity of the material, $C_{\varepsilon}$ is the specific heat per unit mass at constant strain, $n_{dim}$ is the number of spatial dimensions, $\widehat{T}$ is the temperature test function, and $\mathcal{W}_{T}$ is the temperature function space.

For the given thermoelasticity problem (described in Eq. (2.1) and Eq. (2.1)), an FEM solver is developed using the deal.II[53]. We implement the monolithic approach to solve the coupled system of equations, where the solver concurrently computes the displacement ($\boldsymbol{u}_{n+1}$) and temperature ($T_{n+1}$) fields given the response from the previous time step.

2.2 Poroelasticity

Poroelasticity describes the interaction between the mechanical and fluid flow behaviors in a porous medium, where deformations and fluid pressure mutually influence one another. The poromechanical properties of the porous medium can depend on both space and time. For the considered poroelasticity formulation, the unknown variables are the displacements ($\boldsymbol{u}$) and fluid pressure ($p$), coupled through the governing system in Eq. (A.20)-(A.21). Adopting the implicit Euler scheme for time discretization and considering incremental analysis, the weak form of the linear momentum equilibrium equation can be written as:

[TABLE]

where $\alpha^{\prime}$ is the Biot’s coefficient. All other quantities are the same as in the thermoelasticity formulation.

Assuming isotropic properties, the weak form of the fluid flow continuity equation can be written as:

[TABLE]

where $\overline{q}_{f}$ is the boundary fluid flux, $Q_{f}$ represents any fluid source or sink, $M$ is the Biot’s modulus, $K_{I}$ is the isotropic permeability, $\gamma_{f}$ is the weight density, $\mu_{f}$ is the fluid viscosity, $\boldsymbol{i}_{g}$ is the unit vector parallel to gravity, but in the opposite direction, $\widehat{p}$ is the fluid pressure test function, and $\mathcal{W}_{p}$ is the fluid pressure function space. An FEM solver is built using the deal.II library to solve the poroelasticity problem in a monolithic approach. For each time step, the solver, monolithically, satisfies Eq.(2.2) and Eq.(2.2) to compute the displacement ($\boldsymbol{u}_{n+1}$) and pressure ($p_{n+1}$) fields, given the response at the previous time step.

3 I-FENN Formulation and Implementation

3.1 I-FENN Formulation

The Integrated Finite Element Neural Network (I-FENN) framework is a computational methodology introduced to accelerate the numerical simulation of multiphysics problems [44]. An illustration of the I-FENN setup is shown in Fig. 1. For a problem with two governing equations and two coupled variables (displacement field $u$ and another variable $z$), the matrix form of the governing equations can be expressed as shown in Fig. 1. One standard option to solve the system is the FEM monolithic approach (refer to Fig. 1), where the two equations are solved simultaneously. This approach requires assembling the system matrix and the right-hand side vector, including all degrees of freedom, which can be computationally expensive and/or memory-intensive, especially for large-scale problems. I-FENN, on the other hand, is conceptually similar to FEM-staggered approaches, where the two equations are decoupled and solved iteratively (see Fig. 1). I-FENN utilizes a conventional FEM solver for the balance of linear momentum to compute the displacement variable $u$, while the other variable is predicted using a Neural Network (NN). The second coupled field ($z$) is explicitly introduced to the first equation at each time step, therefore constantly informing the displacement field of its updated profile. This hybrid formulation yields several advantages. First, for the FEM-based part, it results in a symmetric system matrix and right-hand side vector that only include the degrees of freedom related to the displacement variable $u$. The symmetry of the system matrix $K_{uu}$ facilitates the convergence of the linear solver. In addition, given the reduced matrix size and degrees of freedom, I-FENN provides faster assembly and lower memory requirements for the FEM solver, while still ensuring the minimization of the displacement residuals. Second, by focusing on predicting only one unknown variable ($z$), which constitutes a single learning task, the NN achieves greater training efficiency than approaches that attempt to predict all variables concurrently [54] (i.e., simulation substitution). Third, since the NN-based prediction of the second variable is much faster than an FEM-based computation, I-FENN yields considerable savings overall compared to either an FEM-monolithic or an FEM-staggered method, which is demonstrated through previous studies [50].

3.2 I-FENN for Thermoelasticity and Poroelasticity

In this section, the I-FENN formulation is applied to two multiphysics problems, where the second coupled variable, denoted $z$, represents the temperature field ($T$) in thermoelasticity and the fluid pressure field ($p$) in poroelasticity. The governing equations of the thermoelasticity problems are listed in Section 2.1. Utilizing FEM, I-FENN solves the equilibrium equation (Eq. (2.1)) and computes the displacement field $\boldsymbol{u}$, while the temperature field $T$ in Eq. (2.1) is predicted using the NN. Equation (2.1) shows that the temperature $T$ is coupled with the displacement $u$ through the trace of the strain tensor $\text{tr}(\boldsymbol{\varepsilon})$. Utilizing this remark enables the use of simplified inputs to the NN. Rather than relying on the full displacement field $\boldsymbol{u}$ (a vector quantity), the model can be trained to predict $T$ based on the trace of the strain tensor $\text{tr}(\boldsymbol{\varepsilon})$, which is a scalar quantity. This simplification reduces the complexity of the input parameters for the NN, making it more efficient and easier to train. In addition to strain trace $\text{tr}(\boldsymbol{\varepsilon})$, thermal loads depicted in Eq. (2.1) are also included as input training parameters. The thermal loads can be domain loads (body loads $r$) or surface loads (Neumann boundary conditions $\overline{Q}$) . Depending on the specific problem setup, either or both types of thermal loads can serve as input features for the NN.

I-FENN for poroelasticity employs a methodology similar to that used in thermoelasticity, illustrating the framework’s adaptability for various multiphysics problems. First, the balance of linear momentum equation (Eq. (2.2)) is solved using an FEM solver, and the continuity equation (Eq. (2.2)) is solved by the NN to predict the fluid pressure field ($p$). By inspecting the coupled equations (Eq. (2.2) and Eq. (2.2)), it can be seen that the pressure field is coupled with the displacement field through strain trace. This implies that the NN needs only this scalar quantity to sufficiently represent the mechanical effects (deformations). In addition to the strain trace, the NN inputs should also represent the fluid flux $(\overline{q}_{f})$ over the boundaries and the fluid sinks or sources $(Q_{f})$ throughout the domain. The design of the NN model (DeepONet) will be presented in detail in Section 4. In this paper, it is important to note the difference in numerical strategies employed in I-FENN versus the fully coupled FEM solver used for data generation and performance comparisons. The I-FENN setup inherently follows a staggered scheme, whereas the fully coupled FEM solver employs a monolithic approach. The choice of a monolithic solver for FEM is driven by the need to achieve optimal performance for the strongly coupled transient multiphysics problems.

3.3 I-FENN Implementation

In this section, we illustrate the integration of the NN of choice (in this case, a DeepONet) within the I-FENN framework. Fig. 2 depicts the framework architecture, which consists of three main layers: (1) Python DeepONet layer, (2) Operating system communication layer, and (3) C++ FEM layer. The DeepONet layer, written in Python, works as a server running on the machine’s graphics processing unit (GPU). It contains the trained DeepONet model, which is loaded once and receives inputs from the communication layer at every time step. The second layer serves as an intermediate communication layer that leverages the operating system’s shared memory. Both the first and third layers utilize libraries for interprocess communication via shared memory and semaphore-based synchronization. In this context, a semaphore is a coordination mechanism that controls access to shared resources by blocking some processes until another process signals that it has completed its operation.

The third layer primarily consists of an FEM solver implementation based on the deal.II library [53], with PETSc [55] integration for large-scale systems parallel processing. The third layer’s FEM solver is used to solve mechanical systems governed by the coupled balance of linear momentum equations, as discussed in Section 3. Specifically, an FEM solver is created to solve Eq. (2.1) and Eq. (2.2) for thermoelasticity and poroelasticity problems, respectively. Depending on the problem, temperature or pressure fields are provided to a specific FEM solver as inputs. These fields are then internally incorporated when solving for the displacement field.

The execution of the framework initiates with the activation of the first computational layer, which monitors a shared memory space and waits for input from the FEM solver. At each discrete time step $n+1$, the FEM solver outputs strain trace data from the previous time step $\text{tr}(\boldsymbol{\varepsilon}_{n})$ and writes it into the shared memory. In a parallel computing environment, the computational mesh is distributed among multiple processors, each contributing its local strain response to the shared memory. Once all processors have successfully completed this data exchange, a synchronization mechanism, implemented via a semaphore, is triggered to signal the DeepONet module that it may proceed.

Upon activation, the DeepONet module performs inference using the strain data $\text{tr}(\boldsymbol{\varepsilon}_{n})$ retrieved from the shared memory in conjunction with the prescribed input loads for step $n+1$. The resulting outputs ($T_{n+1}$ or $p_{n+1}$) are then written back into the shared memory. Another semaphore is subsequently activated to notify the FEM solver that it can resume computation. The FEM solver then reads the updated input fields ($T_{n+1}$ or $p_{n+1}$), computes the corresponding displacement field ($\boldsymbol{u}_{n+1}$), and performs post-processing operations to obtain the strain trace field $\text{tr}(\boldsymbol{\varepsilon}_{n+1})$, thus completing one iteration of the coupled workflow.

This framework is designed with high versatility and flexibility, effectively leveraging both GPU acceleration for the DeepONet model (using PyTorch library) and parallel processing capabilities for the FEM solver (using deal.II library). Notably, the implementation achieves efficient resource utilization by loading the DeepONet inference model into memory only once, even when FEM is executed across multiple processors. This design choice optimizes memory utilization, avoiding memory overhead, and enabling scalability and computational efficiency in distributed systems typically required for computationally demanding multiphysics problems.

4 DeepONet for I-FENN

In this section, we present a detailed overview of the proposed DeepONet, including the operator architecture, GRU design and integration, boundary conditions handling, training setup, loss functions, and testing metrics.

4.1 DeepONets General Architecture

The DeepONet architecture was introduced as a network approximating an operator $(\mathcal{G})$ that can map an input function $(F_{I})$ to an output function $(F_{O})$ [56, 26]. A DeepONet model consists mainly of two sub-networks, denoted as ”branch” and ”trunk”, with each network assigned a specific functionality. As shown in Fig. 3, the branch network is used to encode the discretized input function $F_{I}(x_{I})$, while the trunk encodes the target locations in the output domain $x_{O}$. The final output of the DeepONet network is computed through an element-wise multiplication process defined as:

[TABLE]

where $b_{0}\in\mathbb{R}$ is a bias, $D_{out}$ is the number of outputs from both branch and trunk networks, $b(F_{I})$ and $t(x_{O})$ are the output vectors, of size $D_{out}$, from the branch and trunk, respectively.

For an input function $F_{I}$ and an output function $F_{O}$ that are defined over domains $D_{I}\subset\mathbb{R}^{d_{I}}$ and $D_{O}\subset\mathbb{R}^{d_{O}}$, respectively, the operator mapping can be defined as:

[TABLE]

where $\mathcal{F_{I}}$ and $\mathcal{F_{O}}$ are the spaces for $F_{I}$ and $F_{O}$ defined over domains $D_{I}$ and $D_{O}$, respectively. It should be noted that in this approach, the input function $F_{I}$ is discretized by substitution at a set of locations $\{{x_{I}}_{1},{x_{I}}_{2},...,{x_{I}}_{m}\}$ to get a set of evaluations $\{F_{I}({x_{I}}_{1}),F_{I}({x_{I}}_{2}),...,F_{I}({x_{I}}_{m})\}$. MIONet[32] extends the DeepONet architecture by introducing multiple branches, each processing a separate input function, with their outputs subsequently combined with the trunk output through a tensor product. A further development, Fourier-MIONet [57], introduces element-wise summation, instead of multiplication, to merge branch outputs. More recently, the Image Generator Enhanced Deep Operator Network (IGE-DeepONet) [58] explored hybrid strategies, combining data concatenation with data fusion (transfer) between the hidden layers of the branch and trunk. Notably, in IGE-DeepONet, concatenation-based fusion was shown to achieve higher accuracy than multiplication-based approaches.

Building on these ideas, our method employs multiple branches to accommodate multiple inputs and concatenates their outputs into a unified representation. For clarification, three different architectures are depicted in Fig. 4. First, the stacked DeepONet approach (see Fig. 4a), introduced by Lu et al. [26], employs several branches to encode the same input function, with each branch producing a scalar quantity. Second, MIONet (Fig. 4b) aggregates branch outputs through a tensor product with the trunk output. Finally, our proposed architecture (Fig. 4c), a modified version of MIONet, concatenates the outputs of the individual branches. This design provides flexibility in tailoring the branch network size for each input. Moreover, it preserves the full output of each branch prior to reduction, thereby offering a richer representation for integration with the trunk output. For clarity, we will hereafter refer to this modified version of MIONet as the ”DeepONet model”.

4.2 Gated Recurrent Unit (GRU) Networks

The gated recurrent unit (GRU) architecture was introduced as an extension to the recurrent neural network (RNN). GRU uses a relatively sophisticated approach, named the gating mechanism, to process the temporal dependencies along sequences of data [24, 59]. For an input $\mathbf{x}_{t}\in\mathbb{R}^{d}$, at time step $t$, the GRU output hidden state $\mathbf{h}_{t}\in\mathbb{R}^{h}$, is computed through the following computational steps [60, 61]:

[TABLE]

where symbol $\odot$ is the Hadamard (element-wise) product operator, $\mathbf{r}_{t}$ is the reset gate, $\mathbf{z}_{t}$ is the update gate, $\mathbf{h}^{\prime}_{t}$ is the candidate hidden state, $\sigma$ is the sigmoid activation function, $\mathbf{W}_{xr},\mathbf{W}_{xz},\mathbf{W}_{xh}\in\mathbb{R}^{d\times h}$ and $\mathbf{W}_{hr},\mathbf{W}_{hz},\mathbf{W}_{hh}\in\mathbb{R}^{h\times h}$ are weight parameters, $\mathbf{b}_{r},\mathbf{b}_{z},\mathbf{b}_{h}\in\mathbb{R}^{1\times h}$ are the bias.

4.3 DeepONet with GRUs: detailed implementation

Based on the analysis of the coupled problems and the I-FENN approach introduced in Section 3, a generic architecture of the DeepONet is proposed as depicted in Fig. 5. Load and strain inputs are processed through separate branches, allowing for greater flexibility in network size and enabling independent control over the network’s learning dynamics for each input. A comparable strategy was adopted in the literature [30], involving the use of two distinct branches to differentiate between homogeneous and heterogeneous inputs.

In the proposed model (refer to Fig. 6), each branch incorporates a stacked GRU (multiple consecutive GRUs) to capture temporal dependencies across time steps. The GRU outputs are subsequently passed through a fully connected feedforward multilayer perceptron (MLP) to effectively encode spatial relationships among the discretized inputs. The trunk network only consists of a fully connected feedforward MLP [56].

For the first branch, a load input of size $N_{l}$ is encoded through the GRU, which outputs a vector of size $N_{H\_1}$, which is then processed through the MLP with a final layer size $D_{out\_1}$. The second branch encodes the strain input of size $N_{s}$ through the GRU to get the hidden output of size $N_{H\_2}$. The hidden output in the second branch is split into $N_{ch}$ channels to be normalized through a group normalization layer, then reshaped before proceeding to the MLP, as depicted in Fig. 6. The additional group normalization layer is crucial for the I-FENN framework stability, as discussed in Section 5.1.5.

If the network is trained to predict multiple $N_{c}$ components, not a single component, the output of the branches is set to be $D_{out\_1}\times N_{c}$ and $D_{out\_2}\times N_{c}$, respectively. For example, if the network is trained to predict the temperature field only, $N_{c}$ is set to unity, while if it is trained to predict temperature and all displacement components, $N_{c}$ will be four. Noting that, for I-FENN implementation in this paper, $N_{c}$ will always be unity; however, for a case study on simulation substitution, $N_{c}$ will be set to four. Finally, for each branch, $N_{GRU\_Bi}$ and $N_{FC\_Bi}$ represent the number of stacked GRUs and fully connected layers, respectively, for the $ith$ branch.

Regarding the trunk, an input of size $N_{d}$ (number of coordinates per node) is processed through an MLP with $N_{FC\_T}$ layers, ending with an output of size $D_{out}$. The output of the branches is reshaped and concatenated to have the same dimension as the trunk output $D_{out\_1}+D_{out\_2}=D_{out}$, as shown in Fig. 5. The final output is computed by applying the element-wise multiplication (Eq. (5)) along the output vectors of length $D_{out}$ for $N_{c}$ times to compute the output for all components. Further details on the input and output dimensions, as well as the training procedure, are provided in Section 4.5.

4.4 Enforcing boundary conditions

Lu et al. [26] introduced the following approach to enforce Dirichlet boundary conditions (BCs). A Dirichlet BC $g(x_{O})$ defined over the boundary part $\Gamma_{D}$ can be automatically satisfied by the DeepONet as follows:

[TABLE]

where $\mathcal{N}(F_{I})(x_{O})$ is the resultant of the tensor product of the trunk and branch outputs.

In case of complex BCs where different functions are defined over different boundary parts, additional extensions are required. Assume we have $n$ boundary parts $\{\Gamma_{D_{1}},\Gamma_{D_{n}},...,\Gamma_{D_{n}}\}$ with well defined Dirichlet BC functions $\{g_{1}(x_{O}),g_{2}(x_{O}),...,g_{n}(x_{O})\}$. Sukumar and Srivastava [62] proposed enforcing these functions as follows:

[TABLE]

where $\phi_{i}$ is the approximate distance function for the boundary segment $\Gamma_{D_{i}}$, and $w_{i}$ is the corresponding transfinite interpolation weight. Further details on the construction of these distance functions and the interpolation scheme can be found in literature [62, 63].

Inspired by approaches in Eq. (4.4) and Eq. (14) we propose another extension as follows:

[TABLE]

Fig. 7 depicts the main difference between the original approach introduced by Lu et al. [26] and our approach. In the original approach (Eq. (4.4)), $\ell_{i}(x_{O})$ is zero on the boundary part $\Gamma_{D_{i}}$ and can be any positive value on other boundaries ($\Gamma\setminus\Gamma_{D_{i}}$). In our approach, the function $\ell_{i}(x_{O})$ is constructed to vanish on the boundary part $\Gamma_{D_{i}}$, thereby canceling the contribution of $\mathcal{N}(F_{I})(x_{O})$ on $\Gamma_{D_{i}}$. Simultaneously, it is set to be unity on the remaining boundary $\Gamma\setminus\Gamma_{D_{i}}$, effectively nullifying the influence of $g_{i}(x_{O})$ outside of $\Gamma_{D_{i}}$. However, it should be mentioned that enforcing boundary conditions through the proper selection of $\ell_{i}(x_{O})$ and $g_{i}(x_{O})$ functions can be challenging in the presence of complex geometries, highlighting the importance of the alternative methods[64, 62, 65], and the need for more versatile approaches in future research.

4.5 DeepONet training and testing

The DeepONet model is trained in a data-driven approach where labeled pairs of inputs and outputs are generated before model training. Using the fully coupled FEM solvers, explained in Section 2.1 and Section 2.2, datasets are generated for several load histories (load cases). Afterwards, the DeepONet model training is conducted using the PyTorch library written in Python. For a feed-forward path in the DeepONet, three arrays are used as inputs for the load branch, the strain branch, and the trunk, as depicted in Fig. 5. The shapes of the three input tensors for the first branch, the second branch, and the trunk are $[N_{b},N_{t},N_{l}]$, $[N_{b},N_{t},N_{s}]$, and $[N_{n},N_{d}]$, respectively, where $N_{b}$ is the batch size (number of load cases per batch), $N_{t}$ is the number of time steps, $N_{l}$ is the number of load points (load sensors), $N_{s}$ is the number of strain points (strain sensors), $N_{n}$ is the number of target nodes (locations), $N_{d}$ is the number of spatial dimensions. The trunk encodes the same set of node coordinates for every load case and time step. To eliminate redundancy and enhance the model performance[56], the batch size ($N_{b}$) and number of time steps ($N_{t}$) are excluded from the input dimensions of the trunk. The overall network output for a given batch is a tensor with dimensions $[N_{b},N_{t},N_{n},N_{c}]$. Each feed-forward operation is followed by backpropagation and weight updates to minimize the errors.

With backpropagation, the training process can be described as an optimization task, defined as follows:

[TABLE]

where $\phi$ represents the model trainable parameters, and $\phi^{*}$ is the optimized set of weights to reduce the loss function $\mathcal{L}$. Several definitions for the loss function were tested. Here, we list the main ones that were used for the following examples in Section 5:

[TABLE]

where $\left\|.\right\|_{2}$ is the second norm of a vector, $y_{pred}$ is the predicted response, $y_{true}$ is the fully coupled FEM response considered as the true value.

A clear distinction should be made between the training loss function $\mathcal{L}$ and the testing loss metric. To make consistent and fair comparisons, a unified system for testing error reporting is adopted as follows:

[TABLE]

where $L_{2}^{t}$ is the normalized response error for a specific component at a specific time step $t$, $L_{2}^{LC}$ is the error for a specific load case (LC) over its whole time history, $L_{2}^{All\ LCs}$ is the overall error for all testing load cases. Finally, $\epsilon_{rel}$ is the relative error, with the tolerance $\epsilon_{tol}$ added for numerical stability. Due to having different response components that vary widely in magnitude, ranging from $10^{-5}$ to $10^{5}$, a consistent $\epsilon_{tol}$ is set as $10^{-5}$ of the maximum absolute $y_{true}$ value per component.

5 Numerical Examples

In this section, we assess the proposed framework through different examples demonstrating its flexibility, generalizability, and scalability. Three examples are provided: (1) thermoelasticity model of a 3D cube with domain thermal loads, (2) thermoelasticity model of a 3D thick-walled tube with thermal surface loads, and (3) poroelasticity problem of a 2D excavation setup under fluctuating dewatering flux.

For all examples, a systematic approach is adopted to reach the I-FENN implementation example. First, a sensitivity analysis is conducted to define the initial problem setup, including total time, time step increment, and mesh size. Afterwards, labeled datasets are generated for different load cases. Datasets are split into training, validation, and testing subsets. Extensive optimization trials are then undergone to reach an optimized DeepONet network with reasonable accuracy. However, it should be noted that despite persistent efforts to optimize the DeepONet models, there remains room for improvement and further optimization, often constrained by limitations in resources and time.

After training, certain testing load cases exhibiting unique physics attributes are selected for I-FENN implementation and further model assessment. In the following discussions, these load cases will be referred to as ”featured” load cases. In addition, the DeepONet testing error is evaluated for all the load cases within the testing dataset. The 10th, 50th (median), and 90th error-percentile cases are also chosen for I-FENN implementation. For performance comparisons, I-FENN and fully coupled FEM simulations are conducted on a workstation running Ubuntu 24.04, equipped with an Intel Xeon W3-2435 processor (16 cores), 128 GB of RAM, and an NVIDIA RTX A2000 GPU with 12 GB of dedicated memory.

5.1 Thermoelasticity example: 3D Cube with thermal body load

5.1.1 Problem setup

The first example is a 3D cube with a unit length along each direction. The material properties are set as $\lambda=40$GPa, $\mu=27$ GPa, $\alpha=2.31\times 10^{-5}$ m/(m.K), $C_{\varepsilon}=910$ J/(kg.K), $\rho=2700$ kg/m3, and $k=237$ W/(m.K). As shown in Fig. 8, the cube has zero essential boundary conditions ($\boldsymbol{u}_{0}{}=\boldsymbol{0}$ and $\widetilde{T}=0$) on the left face, where $\widetilde{T}=T-T_{0}$, and the reference temperature $T_{0}=293$ K. On the top face, a temperature boundary condition, function of space and time, is imposed as $\widetilde{T}=10x\gamma$, where $\gamma=\text{time in seconds}/1800\leq 1$. The total duration is $18,000$ seconds, which means that the temperature on the top face stabilizes after 10% of the total duration. A Gaussian random thermal body load is applied throughout the whole domain. The simulation is conducted over uniformly spaced 100 time increments.

5.1.2 Training and testing results

A total of 1000 load cases were generated: 900 of these cases used a coarse mesh of $11\times 11\times 11$ nodes ($N_{n}=1331$) and were split into training and validation sets with a ratio of 4:1. The remaining 100 load cases were reserved for testing using a fine mesh $(31\times 31\times 31)$ nodes. Table 1 lists the hyperparameters describing the DeepONet model architecture. An input of size 512 is used for both load and strain branches. The input points (sensors) are spaced equally in a space grid, with 8 points per direction. A group normalization layer is added to the second branch to stabilize the framework’s performance, as discussed in Subsection 5.1.5.

Training is done over $24,000$ epochs with the $L_{2}$ norm function $\mathcal{L}_{L_{2}}$ selected for loss evaluation. Training vs. validation loss values are depicted in Fig. 9a, where both curves are decreasing fairly smoothly, showing that the model is learning effectively. The training loss is relatively lower than the validation curve, which is typically expected. Training required a total time of 9 hours and 4 minutes. Although increasing the batch size could further reduce training time, we kept the current configuration (batch size = 16) because earlier experiments showed that larger batches negatively affected accuracy. That said, there are still several avenues for improving training efficiency. Potential directions include recalibrating the training hyperparameters, applying mixed precision, experimenting with alternative optimizer variants, or adjusting data loading configurations [66, 67].

After training on a coarse mesh, the model is tested for 100 load cases (LCs) using a fine mesh of equally spaced $N_{n}=29791$ nodes (31 per direction). A histogram of the normalized $L_{2}$ norm of error per load case ($L_{2}^{LC}$) is plotted in Fig. 9b. The median, 10th, and 90th error-percentile load cases are selected for I-FENN integration. The results for the featured and the median load cases are presented in the following subsection, while the 10th and the 90th error-percentile load cases are presented in the appendix (B.1.1).

5.1.3 I-FENN results

The trained model is incorporated into the I-FENN framework as described in Section 3. The I-FENN results for the median load case are presented in Fig. 10. In addition, a featured load case is chosen based on the highest thermal load integral, representing the maximum amount of heat change introduced to the body. Throughout all the testing load cases, the selected featured load case exhibited the maximum absolute temperature change ($\widetilde{T}$) and maximum absolute strain trace. The I-FENN results for the featured load case are presented in Fig. 11, with maximum response values detected at the 86th time step.

In both figures (Fig. 10 and Fig. 11), the coupled FEM results are considered as the true values ($y_{true}$), which are plotted in the top row, followed by the I-FENN response ($y_{pred}$) in the second row. The third row presents the relative error ($\epsilon_{rel}$) capped at 5%. The first column depicts the temperature changes ($\widetilde{T}$), followed by the displacement components ($u_{x},u_{y},\text{and }u_{z}$). The I-FENN results demonstrate strong agreement with the FEM results. Overall, the results indicate excellent consistency, particularly for zones with extreme response values, which exhibit relative errors well below 5%. The relative error zones exceeding 5% are attributed to very small response values close to zero, which magnifies the relative error values. A consistent behavior of I-FENN is found in the results for the 10th and 90th error-percentile load cases, depicted in Fig. B.1 and Fig. B.2, respectively.

5.1.4 Computational savings and scalability

The network was trained using a coarse mesh of $(11\times 11\times 11)$ nodes, and then accurately predicted the response for testing on a fine mesh of $(31\times 31\times 31)$ nodes. A complementary analysis is provided in Appendix B.1.2, where Fig. B.3 confirms that the framework maintains its accuracy on further mesh refinement. To assess performance and scalability, the framework’s efficiency in computing responses on finer meshes is examined. The results are tabulated in Table 2, showing the computational time of the I-FENN framework and the fully coupled FEM solver for different domain sizes (cells per direction).

For all tested domain sizes, the I-FENN framework exceeded the performance of the fully coupled FEM solver, showing a consistent computational efficiency. To better visualize the results, computational times and savings across different domain sizes are plotted in Fig. 10. In addition to the demonstrated computational savings, this plot highlights the scalability of the I-FENN framework as the domain is refined. This scalability is a critical feature for a framework designed to improve performance in multiphysics simulations, ensuring its applicability to more complex and computationally demanding problems in the future. It is important to highlight that the reported computational times of I-FENN exclude the training phase of the DeepONet. For just a few simulations, on a coarse mesh, the overall time (including DeepONet training and I-FENN execution) can exceed that of running a fully coupled FEM solver. However, thanks to I-FENN’s demonstrated scalability, when performing many simulations for different load cases, especially on fine meshes, the combined cost becomes significantly lower than that of the traditional FEM approach.

I-FENN outperforming the fully coupled FEM is a particularly noteworthy observation, given that the deal.II library, written in C++, is highly optimized for performance and designed with high-performance computing (HPC) compatibility in mind. In contrast, the I-FENN framework is still under active development, with numerous performance optimizations yet to be implemented.

5.1.5 I-FENN stability

During the initial stages of developing the network, it is evident that there was an issue with the output stability. There was an error accumulation pattern through time steps, resulting in higher errors towards the end of the simulation. To address this issue, a modified framework design was established to investigate the root cause of the error.

In this approach, true values (computed using the fully coupled FEM solver) were incorporated into the framework. The simulation starts with using true strain trace and temperature values as inputs to the DeepONet and the mechanical FEM solver, respectively. After step number 33, the DeepONet output replaces the true temperature and is sent to the FEM solver. At step number 66, the strain trace output from the mechanical solver replaces the true strain and is used as input to the DeepONet.

This testing approach was applied to different networks with different architectures, with results summarized in Fig. 13. The figure depicts the normalized error norm per time step $L_{2}^{t}$ for strain trace and temperature values across time steps as the simulation progresses incrementally from one step to another. Four design options are tested for the strain branch: a) GRU hidden output of size 200 ($N_{H\_2}=200$) without group normalization, b) $N_{H\_2}=200$ with group normalization channels of $N_{ch}=50$, c) $N_{H\_2}=200$ with $N_{ch}=100$, and d) $N_{H\_2}=50$ with $N_{ch}=25$. The reader is referred to Section 4.3 for more details on the adopted group normalization approach.

For the first design option with no normalization, the error values of strain trace and temperature show different patterns, which can be described along three phases. Phase one starts from step 1 to step 33. The strain trace relative error is very low, given that this is a response of the mechanical FEM solver using true temperature inputs. In the first stage, also, the DeepONet output (temperature), shows a relatively reasonable error which is dependent on the network training.

Moving to the second phase between steps 33 and 66, only one parameter is changed. The temperature predictions from DeepONet are used as inputs to the mechanical FEM solver, resulting in a corresponding stable increase in the strain trace error. The third stage starts at the 66th step by using the FEM mechanical strain trace as input to the DeepONet. Therefore, it can be concluded that the DeepONet model is susceptible to minor changes in the strain trace inputs.

To overcome this sensitivity to changes in the strain input, the group normalization approach is introduced in the other three design options. Design options (c) and (d) show enhanced stability in the network performance. Option (d) provided an optimized network size and, hence, is selected for further implementation (refer to Table 1). Finally, an additional study on the 3D cube example is conducted to investigate the impact of not enforcing boundary conditions. For brevity, the detailed analysis and corresponding results are presented in Appendix B.1.3.

5.2 Thermoelasticity example: 3D Thick-walled tube with thermal surface load

5.2.1 Problem setup

The second example is a 3D thick-walled tube, as depicted in Fig. 14, with inner and outer radii of 1.0 m and 2.0 m, respectively. The tube has a length of 1.0 m along its longitudinal axis (Z-axis). Zero essential boundary conditions are imposed on the base plane ($z=0$). Material properties are the same as those selected for example 1 in Section 5.1.1. However, for this example, no body loads were incorporated. Instead, boundary flux values are applied on the inner and outer walls of the tube. Both inner flux $(q_{in})$ and outer flux $(q_{out})$ values are defined as functions of time and angular location $(\theta=\arctan(y/x))$. Flux functions are defined as $q_{in}(\theta,t)=q_{in-0}+q_{in-1}\sin{(\omega_{t}t)}\sin{(\omega_{in-r}\theta)}$, and $q_{out}(\theta,t)=q_{out-0}+q_{out-1}\sin{(\omega_{t}t)}\sin{(\omega_{out-r}\theta)}$.

5.2.2 Training and testing results

The parameters $q_{in-0}$, $q_{in-1}$, $q_{out-0}$, $q_{out-1}$, $\omega_{in-r}$, $\omega_{out-r}$, and $\omega_{t}$ were randomly sampled to generate $1000$ unique loading cases. Each case is simulated for $60000$ seconds, divided into 120 equal increments. A subset of $900$ load cases was split into training and validation sets with a ratio of $4:1$, while the remaining $100$ load cases were used for testing. Comprehensive parametric studies, detailed in Appendix  B.2, are conducted to optimize the model’s hyperparameters. Based on these studies, the hyperparameters listed in Table 3 are selected for training the DeepONet model employed in I-FENN.

A total of 128 equally spaced points are selected for surface load input, with 64 points per surface. Within the whole domain, 512 structurally spaced points are selected to represent the strain trace field. A total number of $5184$ nodes is utilized for training the network ($N_{n}=5184$). A sum of square error loss function $\mathcal{L}_{\text{SSE}}$ is utilized for calculating training loss for a total number of $24,000$ epochs, requiring 12 hours and 35 minutes for training. Training vs validation loss curves are provided in Fig. 15a, showing a stable descending learning curve with no signs of overfitting.

A histogram of testing error is depicted in Fig. 15b, showing a roughly normal error distribution with slight skew towards the left. The median, 10th, and 90th error-percentile load cases are selected for I-FENN implementation. The median load case is covered in the following subsection, while the 10th and 90th error-percentile load cases are provided in Appendix B.2.7.

5.2.3 I-FENN results

I-FENN implementation is done as described in Section 3, incorporating the trained DeepONet model. I-FENN results and the fully coupled FEM response for the median load case are depicted in Fig. 16 along with the relative error values. The coupled FEM results are plotted in the top row, as the true values ($y_{true}$), followed by the I-FENN response ($y_{pred}$) in the second row, and the relative error ($\epsilon_{rel}$) capped at 5% at the third row. All I-FENN response fields show a good agreement with the fully coupled FEM response, with relative error values much lower than 5% at the zones with high response values. Only zones with response values close to zero exhibited relative error exceeding 5%.

In addition to the median load case, a featured load case is chosen based on the highest discrepancy in the spatial frequency between inner and outer boundary fluxes. Figure. 17 presents the results of the featured load case at the 90th time step, showing the ability of I-FENN to capture extreme response spacial discrepancies with high accuracy. The I-FENN framework continues to demonstrate consistent and robust performance across this example and the previous one (Section 5.1). The current example further highlights the versatility of the I-FENN approach, showcasing its ability to handle flux values that vary both spatially and temporally within a more complex geometry than that of the previous example.

5.2.4 I-FENN vs. simulation substitution

As discussed in the literature review (Section 1.1), a possible approach for machine learning in computational mechanics is simulation substitution, where the whole response is predicted using a neural network as a surrogate model. To test this approach, a DeepONet model is trained for predicting all components as a surrogate model. The surrogate model has only one branch for load, as the strain trace will not exist as an output from the mechanical solver. For the surrogate model, $N_{c}$ is set as four to represent the four response components ($\widetilde{T},u_{x},u_{y},\text{and }u_{z}$). The surrogate model hyperparameters (listed in Table 4) are calibrated to provide a roughly close number of trainable parameters as those of I-FENN integrated DeepONet.

The DeepONet for I-FENN and the surrogate model have a number of trainable parameters of 372,032 and 375,740, respectively. Labeled pairs of input/output data were used for training the model, with four output components per node. Then, the surrogate model is trained using the same parameters used for training the DeepONet for I-FENN, such as training/testing load cases, the loss function, the learning rate, and the number of epochs.

The trained surrogate model results for the median load case at the last time step are depicted in Fig. 18. Noting that the median load case here refers to the same load case detected from the frequency distribution reported in Fig. 15b. Compared to I-FENN results (Fig. 16), the surrogate model exhibits higher error values across all components, with larger zones exceeding 5% relative error. For a more comprehensive assessment, the testing error $L_{2}^{t}$ values for I-FENN and the surrogate model for all time steps and response components are plotted in Fig. 19. Error plots show that I-FENN errors are almost one order of magnitude less than those of the surrogate model, especially for the displacement components. It is worth noting that the higher accuracy in displacement components is related to the mechanical solver integrated within the I-FENN framework, highlighting the accuracy edge provided by I-FENN compared to surrogate models.

Compared to I-FENN, the surrogate model requires four times larger output training datasets, making it considerably more demanding, particularly for more complex or high-dimensional problems. In addition to higher memory requirements, training the surrogate model took twice as long as training the I-FENN integrated DeepONet. Additionally, the surrogate model functions as a black box, offering neither interpretability nor flexibility, and lacks the capability for deeper analysis, such as the stability study presented in Section 5.1.5. While the surrogate model may be acceptable in specific contexts, the superior accuracy, interpretability, and adaptability of the I-FENN framework make it a significantly more appropriate choice for multiphysics simulations.

Finally, it is essential to comment on the relatively early prediction errors, shown in Fig.19. These errors are primarily concentrated within the initial steps and can be attributed to the transient phase introduced by the GRU network. This transient phase arises from the sensitivity of the GRU to its initial hidden states, which are typically initialized to random values or zeros [68]. To mitigate this issue, some researchers employ a technique named as the washout technique, where the first few prediction time steps are excluded from training, allowing the network to update the hidden states and washout the influence of the arbitrary initial conditions [69, 70, 71] However, the washout approach introduces its own drawbacks, including unreliable early predictions and potential training instability [68]. In our approach, the entire sequence is used during training, and although moderate errors appear in the early steps, they remain within acceptable bounds, as discussed in Appendix B.2.1. Importantly, as shown in Fig. 19, these initial errors did not propagate or degrade performance in subsequent predictions.

5.3 Poroelasticity example: 2D Excavation problem with hydraulic boundary flux (Dewatering)

5.3.1 Problem setup

Previous examples have demonstrated I-FENN’s performance and accuracy in thermoelasticity problems. In this example, we extend our validation set to poroelasticity problems. We propose a poroelasticity problem setup inspired by geotechnical engineering problems in literature [72, 73]. Figure 20 depicts an excavation setup where a dewatering system is applied to reduce the water level within the excavated part.

The porous medium properties are set as $\lambda=8.375$ MPa, $\mu=5.58$ MPa, $K_{s}=55556$ MPa, $K_{f}=2200$ MPa, $\phi=0.4$, $\mathbf{K}_{H}=I\times 0.0001$ m/s. A fluctuating water flow $q_{W}(t)$ is assumed at the excavated part water level. For this purpose, Gaussian random histories are generated to represent different scenarios of dewatering flux. The simulation is done over a time span of $3\times 10^{6}$ seconds, discretized over 60 equally spaced increments.

5.3.2 Training and testing results

A total of 700 dewatering flux histories are utilized as follows: 600 load cases for training and validation with ratio of 4:1, and 100 load cases for testing. A model was trained with the hyperparameters listed in Table 5. The input to the load branch varying in time, while spatially it is a scalar field representing a uniform flux value along the excavated surface. The input to the strain branch is a vector 584 strain trace points uniformly distributed along the 2D mesh. The model is trained to predict output for $N_{n}=2309$ Nodes, with trunk input size $N_{d}$ representing the nodes’ coordinates. Notably, no enforcement of boundary conditions was applied for this model. Further details on the training stabilization are provided in Appendix B.3.1.

A total of 8000 epochs are employed, requiring 18 minutes, for model training. Training and validation loss curves ($\mathcal{L}_{\text{SSE}}$) are depicted in Fig. 21a. Both curves show enhanced learning, with the curves descending through epochs with no significant signs of overfitting. The testing loss distribution, shown in Fig. 21b, indicates relatively low error values overall, with an apparent left skew. This suggests that the majority of cases exhibit minimal error, while a smaller number of instances have higher loss values. The load cases with the median, 10th, and 90th error-percentile loss values are selected for I-FENN implementation. The following section will provide the results for the median load case, while the 10th and 90th error-percentile ones are presented in Appendix B.3.2.

5.3.3 I-FENN results

Figure 22 depicts the I-FENN and the fully coupled FEM solver outputs for the median load case at the last time step. The top row displays the coupled FEM results as the ground truth values ($y_{true}$), the second row shows the I-FENN predictions ($y_{pred}$), and the third row presents the relative error ($\epsilon_{rel}$), which is limited to a maximum of 5%. In addition to the median load case, a featured load case is chosen for I-FENN implementation. The selected featured load case exhibited the highest and lowest recorded pressure values among all the testing load cases. Results at time steps with maximum and minimum pressure values are depicted in Fig. 23 and Fig. 24, respectively. Overall, all I-FENN results show a good match with the fully coupled FEM results, with relative errors much lower than 5%. Only zones with response values close to zero show high relative error exceeding 5%.

Figure 25 presents the normalized dewatering flux and the corresponding error for all time steps of the featured load case. The dewatering flux (shown in Fig. 25a) is almost monotonic with low frequency, enabling the system to exhibit the extreme response values. Error curves in Fig. 25b show small error values across all time steps, with displacement errors one to two orders of magnitude less than the pressure errors. This pattern is consistent with I-FENN errors depicted in Fig. 19a. The displacement field, computed using the mechanical FEM solver, shows lower errors than the coupled field predicted by DeepONet. This highlights the effectiveness of the I-FENN approach in maintaining high accuracy for the mechanical field (predicted using FEM), even when the coupled field predictions (using DeepONet) are comparatively less accurate. These observations offer basis for future research on understanding how errors scale in I-FENN simulations for multiphysics problems.

6 Summary and conclusions

In this paper, we present the I-FENN framework for multiphysics problems using DeepONets for thermoelasticity and poroelasticity simulations. A GRU-enabled DeepONet architecture is introduced for integration within the I-FENN framework, handling sequence inputs of time history excitations. The integrated I-FENN framework is introduced in a versatile approach, suitable for deployment on high-performance computing (HPC) systems with support for parallel processing. Three distinct examples covering thermoelasticity and poroelasticity problems are introduced with different domain and surface loading conditions. Based on the analysis and discussions, the key findings can be summarized as follows:

• •

Integration of the GRU-enabled DeepONets within the I-FENN framework provided high flexibility in terms of modeling different geometries and loading conditions in different multiphysics problems.

• •

The I-FENN accuracy for different problems was less than 5% of relative error for all non-trivial (non-zero) points in the domain.

• •

I-FENN predictions for mechanical field displacements demonstrated superior accuracy, even in the presence of higher errors in the coupled field (temperature or pressure predicted by DeepONet). This underscores the I-FENN’s robustness, as mechanical equilibrium is rigorously enforced through FEM, ensuring minimal error in the mechanical field despite inaccuracies in the coupled field. These observations open up new avenues for research to investigate the error and scalability of I-FENN simulations for large-scale multiphysics problems.

• •

I-FENN simulations required significantly less training time and training data and yielded significantly lower errors than simulation substitution (surrogate model) approaches.

• •

The I-FENN framework shows the ability to accurately predict responses for finer meshes despite being trained on a coarse mesh.

• •

I-FENN demonstrated superior performance compared to the fully coupled finite element solver, achieving up to 40% savings in computational cost when applied to finer meshes, highlighting its scalability and efficiency.

Acknowledgments

This work was partially supported by the Sand Hazards and Opportunities for Resilience, Energy, and Sustainability (SHORES) Center, funded by Tamkeen under the NYUAD Research Institute Award CG013. The authors also acknowledge the support of the NYUAD Center for Research Computing, which provided resources, services, and staff expertise that contributed to this work.

Data Availability

The data and source code supporting the findings of this study will be made publicly available upon publication of the article.

Appendix A Mathematical formulation for multiphysics problems

A.1 Thermoelasticity

A.1.1 Strong Form

The governing system of differential equations can be defined as follows [74, 75, 76]:

[TABLE]

where $\boldsymbol{\sigma}$ denotes the Cauchy stress tensor, $\boldsymbol{b}$ is the mechanical body force vector, $\rho$ is the material mass density, $T_{o}$ is the reference temperature, $\eta$ is the entropy per unit mass, $\boldsymbol{q}$ is the heat flux vector, and $r$ is the heat body source or sink. $\Omega$ denotes the physical domain (depicted in Fig. A.1), while $\boldsymbol{\nabla}$ is the gradient operator and $\boldsymbol{\nabla}\cdot$ is the divergence operator.

The boundary conditions are defined as follows:

[TABLE]

where the domain boundary ($\Gamma$) is defined as $\Gamma=\Gamma_{t}\cup\Gamma_{u}=\Gamma_{T}\cup\Gamma_{q}$. Symbols $\boldsymbol{u}$ and $T$ represent the displacement and temperature fields, with their boundary conditions imposed on boundary parts $\Gamma_{u}$ and $\Gamma_{T}$, respectively. Traction ($\boldsymbol{t}$) is defined as $\boldsymbol{t}=\boldsymbol{\sigma}\cdot\boldsymbol{n}$ on the traction boundary part $\Gamma_{t}$, where $\boldsymbol{n}$ represents the normal unit vector at the boundary $\Gamma$. Finally, the flux boundary condition is defined as $Q=\boldsymbol{q}\cdot\boldsymbol{n}$ on the flux boundary part $\Gamma_{q}$.

The heat flux vector $\boldsymbol{q}$ can be defined using Fourier’s law as [74, 75, 76]:

[TABLE]

where $k$ is the thermal conductivity of the material. Thermoelasticity constitutive laws for stress ($\boldsymbol{\sigma}$) and entropy ($\eta$) are expressed as [74, 75, 76]:

[TABLE]

where $C_{\varepsilon}$ is the specific heat per unit mass at constant strain, and $\boldsymbol{\hat{C}}$ is the elasticity tensor defined as $\hat{C}_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)$ for isotropic linear materials. While $\boldsymbol{\varepsilon}$ is the strain defined as $\boldsymbol{\varepsilon}=\boldsymbol{\nabla}^{s}\boldsymbol{u}=\frac{1}{2}(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{T})$ for small deformations. Parameters $\mu$ and $\lambda$ are Lamé constants, and $\alpha$ is the coefficient of thermal expansion.

Substituting Eq. (A.6) in Eq. (A.2) yields:

[TABLE]

A.1.2 Weak Form

Using the balance of linear momentum equation in (A.1) and the stress constitutive relation (A.5), the mechanical weak form can be written as:

[TABLE]

where $\widehat{\boldsymbol{w}}$ is the displacement test function, and $\mathcal{W}_{u}$ denotes the displacement function space. Substituting the Fourier definition of flux (A.4) into the energy equation (A.7), the thermal weak form can be written as:

[TABLE]

where $\widehat{T}$ is the temperature test function, and $\mathcal{W}_{T}$ is the temperature function space. Using the definition of the elasticity tensor for isotropic linear materials $(\hat{C}_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}))$, the weak form of the energy equation (Eq. (A.1.2)) can be written as

[TABLE]

where $n_{dim}$ is the number of spatial dimensions.

A.2 Poroelasticity

A.2.1 Strong Form

For a representative elementary volume composed of solid grains and pores filled with a fluid, the governing equations can be written as [77, 78, 79, 80, 81]:

[TABLE]

where $\boldsymbol{\sigma}_{t}$ is the total Cauchy stress due to deformation and fluid pressure ($p$), $\zeta$ is the increment of fluid content, $\boldsymbol{v}_{f}$ is the fluid velocity vector, and $Q_{f}$ represents any fluid source or sink. Domain boundary is defined as $\Gamma=\Gamma_{u}\cup\Gamma_{t}=\Gamma_{p}\cup\Gamma_{q}$ (see Fig. A.2), with boundary conditions described as follows:

[TABLE]

where $\boldsymbol{\overline{u}}$ and $\boldsymbol{\overline{t}}$ are the displacement and traction boundary conditions defined over boundaries $\Gamma_{u}$ and $\Gamma_{t}$, respectively. $\overline{p}$ is the fluid pressure over boundary $\Gamma_{p}$, while $\overline{q}_{f}$ is the fluid flux over boundary $\Gamma_{q}$, defined as $\overline{q}_{f}=\boldsymbol{v}_{f}\cdot\boldsymbol{n}$.

Assuming fluid flow is controlled by Darcy’s law, fluid velocity and increment can be defined as [77, 79, 80, 81]:

[TABLE]

where $\boldsymbol{K}_{H}$ is the Hydraulic conductivity tensor defined as $\boldsymbol{K}_{H}=\gamma_{f}\boldsymbol{K}_{I}/\mu_{f}$. Given that, $\boldsymbol{K}_{I}$ is the permeability tensor, $\gamma_{f}$ is the weight density, and $\mu_{f}$ is the fluid viscosity. $H_{f}$ is fluid total head defined as $H_{f}={p}/{\gamma_{f}}+\boldsymbol{x}\cdot\boldsymbol{i}_{g}$, where $\boldsymbol{x}$ is the coordinate vector, and $\boldsymbol{i}_{g}$ is the unit vector parallel to gravity, but in the opposite direction.

Finally, $\alpha^{\prime}$ and $M$ are the Biot’s coefficient and modulus, respectively, defined as [77, 79, 80, 81]:

[TABLE]

where $\phi$ is porosity of the medium, $K_{s}$ is the solid grain bulk modulus, $K_{f}$ is the fluid bulk modulus, and $K_{b}$ is the bulk moduli of the solid skeleton defined as [79]:

[TABLE]

Finally, the total stress $\boldsymbol{\sigma}_{t}$ is computed as a function of the solid effective stress $\boldsymbol{\sigma}^{\prime}$ and fluid pressure $\boldsymbol{\sigma}_{f}$, as follows:

[TABLE]

Substituting Eq. (A.2.1) into Eq. (A.11), the linear momentum equilibrium equation for the two-phase mixture can be written as:

[TABLE]

Incorporating definitions in Eq. (A.14) and Eq. (A.15) into Eq. (A.12), the fluid flow continuity equation can be written as:

[TABLE]

A.2.2 Weak Form and FEM

Given the final definitions, the linear momentum equilibrium equation Eq. (A.20) and the fluid flow continuity equation Eq. (A.21), their corresponding weak forms can be respectively listed as:

[TABLE]

[TABLE]

where $\widehat{\boldsymbol{w}}$ is the displacement test function, $\mathcal{W}_{u}$ is the displacement function space, $\widehat{p}$ is the fluid pressure test function, and $\mathcal{W}_{p}$ is the fluid pressure function space.

Appendix B Additional results and parametric studies

This section presents additional results and parametric studies conducted for the three examples in Section 5, offering further insights that may be of interest to the reader.

B.1 Thermoelasticity example: 3D Cube with thermal body load

B.1.1 I-FENN results for the 10th and 90th error-percentile load cases

I-FENN results for the 10th and 90th error-percentile load cases are depicted in Fig. B.1 and Fig. B.2, respectively. The results show consistent accuracy compared to the median load case results presented in Fig. 10.

B.1.2 I-FENN results for the median load case considering a finer mesh

The I-FENN results for the median load case applied on a fine mesh of $75\times 75\times 75$ cells are presented in Fig. B.3. The error limits and distributions closely match that of the coarser mesh of $30\times 30\times 30$ cells (see Fig. 10). This observation further confirms the framework’s ability to preserve accuracy under mesh refinement.

B.1.3 I-FENN results for the median load case without enforcing the boundary conditions

To verify the validity of the boundary conditions enforcement approach, presented in Section 4.4, a DeepONet model is trained without enforcing boundary conditions. Afterward, the trained model is integrated into the I-FENN framework and tested for the median load case, identified in Section 5.1.2. Figure B.4 presents the I-FENN results for the tested configuration in which temperature boundary conditions were not enforced through the DeepONet. The results show a good agreement with the true values; however, the zones with relative error exceeding 5% are much larger than those with enforced boundary conditions (see Fig. 10). This observation underscores the added value of enforcing boundary conditions and suggests that, despite the potential complexity, they should be applied whenever feasible to improve accuracy.

B.2 Thermoelasticity example: 3D Thick-walled tube with thermal surface load

B.2.1 Error in early stages

Error plots in Fig. 19 revealed a relatively high error in the early stages of the simulation. To investigate the cause, we plot the norms of true values ($y_{2}^{t}=\left\|y_{true}\right\|_{2}$) and norms of errors ($e_{2}^{t}=\left\|y_{true}-y_{pred}\right\|_{2}$), as shown in Fig. B.5. In the initial time steps, the true response norm $y_{2}^{t}$ is relatively small, while the error norm $e_{2}^{t}$ is relatively high compared to the later time steps. These two factors are participating in inflating the relative error $L_{2}^{t}$ in the early stages of the simulation (see Fig. 19), since relative error $L_{2}^{t}$ is computed as a ratio and becomes more sensitive when the denominator (true norm $y_{2}^{t}$) is small.

To investigate the extent of the error in initial time steps, I-FENN results for the median load case applied to the 3D tube are plotted at the 1st and 5th time steps in Fig. B.6 and Fig. B.7, respectively. It can be seen that at the 1st time step, the relative error values $\epsilon_{rel}$ at zones with maximum response are significantly lower than 5%. However, areas with response values close to zero show large boundaries exceeding 5%. Moving to the 5th time step Fig. B.7, the results show a consistent behavior to that shown at the 120th time step (see Fig. 16), indicating a stabilized behavior. In conclusion, the figures introduced in this section (Fig. B.5, Fig. B.6, and Fig. B.7) suggest that despite having initial high $L_{2}^{t}$ at early time steps, the discrepancies are mainly concentrated in areas with response close to zero. Moreover, after almost 5 time steps (5% of the whole simulation duration), the effects of early-stage errors introduced by the GRU appear to diminish significantly.

B.2.2 Training dataset size

One of the important parametric studies carried out during the development of the current model is the dataset size and its effect on the model’s generalizability. In this study, different dataset sizes (load case count) are considered for training the model, ranging from 100 to 900 load cases. A consistent set of 100 load cases is used for testing all the trained models. The validation loss during training and the final testing loss values are presented in Fig. B.8. Both plots exhibit a consistent trend: error values decrease as the dataset size increases, confirming that increased dataset size enhances the model’s learning.

All models are trained for the same number of epochs, which raises the question of whether the improved accuracy is primarily due to the increased number of load cases seen per epoch, rather than the diversity of those cases. To address this, a complementary study is conducted. In this study, the model trained with 200 load cases is trained four times the number of epochs used for the model with 800 load cases, and the model with 400 load cases is trained for twice as many epochs as the model trained with 800 load cases.

The validation loss during training and the final testing loss values for this study are presented in Fig. B.9. The results indicate that increasing the number of training epochs alone does not lead to the same improvement in model accuracy as increasing the dataset size. This suggests that the model benefits from the increasing load cases in the training data more than additional exposure to the same data, confirming the data diversity and its role in the learning process.

B.2.3 Training stabilization

One of the challenges that existed in the early stages of training the model was the training stability, where the loss curves showed high fluctuation along the epochs. To tackle this issue, a custom learning scheduler is adopted to stabilize training without impacting the training accuracy or the need for increased training epochs. Three different cases of learning rate options and their corresponding validation loss curves are shown in Fig. B.10. The results show that using a custom scheduler enhanced the learning stability, showing less noise than the high learning rate model and more accuracy than the low learning rate model. The custom scheduler option is adopted in developing the DeepONet model used for I-FENN integration.

B.2.4 Normalization

Various normalization strategies were evaluated, and three key approaches are highlighted here:

• •

(a) min-max normalization to scale the data between 0 and 1, expressed as:

[TABLE]

• •

(b) min-max normalization to scale the data between -1 and 1, expressed as:

[TABLE]

• •

(c) standardization using the mean and standard deviation of the data, expressed as:

[TABLE]

The validation losses vs epochs using different normalization techniques are provided in Fig. B.11a, showcasing training stability for all options. However, accuracy can not be inferred using loss values as they are computed during the training phase with different normalization techniques, which inherently shift the data range, affecting normalized error computation. To address this discrepancy, during the testing phase, all outputs are denormalized before computing the loss metric $L_{2}^{All\ LCs}$ described in Section 4.5. Given that the lowest testing error is given by normalizing between -1 and 1 (refer to Fig. B.11b), this approach is adopted for the DeepONet trained for I-FENN implementation for the 3D tube example.

B.2.5 Loss function

Different loss functions were tested in the initial stage of developing the DeepONet model. The reader is referred to Section 4.5 for definitions of those loss functions. Figure B.12 shows the validation and testing loss using different loss functions during the early stages of model training trials. The loss function $\mathcal{L}_{\text{SSE}}$ showed the least testing error and, hence, was selected for developing the final model.

During the development of the final model, and due to different parametric studies, many hyperparameters were changed, including learning rate, batch size, and dataset size. On refreshing the comparison between different loss functions for the latest training setup, we find that the previous conclusion holds no more, as shown in Fig. B.13. For the latest setup, the best-performing loss function is $\mathcal{L}_{L_{2}}$ compared to $\mathcal{L}_{\text{SSE}}$ in the initial training trials. This finding highlights the dynamic nature of NNs and how changes in hyperparameters can affect the effectiveness of other training hyperparameters.

B.2.6 Branch and trunk size

Among the parametric studies conducted for hyperparameter selection, different values for the hidden layer size $(N_{H})$ for the branch and trunk were investigated. Fig. B.14 shows the effect of branch hidden size on validation and testing errors. Alternatively, Fig. B.15 depicts the effect of branch hidden size on validation and testing errors. Results show that no significant accuracy improvement is detected on increasing the branch hidden size. However, increasing the trunk hidden size showed a substantial increase in accuracy, highlighting the complexity of geometry being embedded through the trunk. Given these results, hidden layer sizes of 64 and 256 are selected for the branch and trunk, respectively, as listed in Table 3.

B.2.7 I-FENN results for the 10th and 90th error-percentile load cases

I-FENN results for the 10th and 90th error-percentile load cases are depicted in Fig. B.16 and Fig. B.17, respectively. The accuracy behavior of the framework is consistent with the median load case results shown in Fig. 16, with much less error for the 10th error-percentile load case. Results demonstrate the reliability of the I-FENN framework across varying loading conditions.

B.3 Poroelasticity example: 2D Excavation problem with hydraulic boundary flux (Dewatering)

B.3.1 Training stabilization

Training stability was one of the main challenges encountered during the early stages of model development, as the loss curves exhibited significant fluctuations across epochs. A custom learning rate scheduler was implemented, aiming to stabilize training without compromising accuracy or requiring additional training epochs. Three different learning rate choices with their corresponding validation loss curves are depicted in Fig. B.18. The results highlight the importance of using a training scheduler that allows broad exploration in the early stages through a high learning rate, followed by a gradual reduction to enable precise fine-tuning as training progresses. The custom scheduler learning rate shown in Fig. B.18a is the adopted scheduler for the DeepONet model used in the I-FENN implementation, with results in Section 5.3.3

B.3.2 I-FENN results for the 10th and 90th error-percentile load cases

I-FENN results for the 10th and 90th error-percentile load cases of the 2D excavation example are presented in Fig. B.19 and Fig. B.20, respectively. A consistent accuracy is detected compared to the median load case presented in Fig. 22, highlighting the I-FENN’s reliability for different loading conditions.