Graph Neural Networks for Full Waveform Inversion

1. Introduction

Full waveform inversion (FWI) has its origins in seismic applications  [10,11]. In an inversion framework, waveform information is exploited by iteratively fitting simulated data from a numerical model to measurements. In recent years, the application of FWI has also been explored in ultrasonic testing, where it can be applied for internal damage detection in solid structures thanks to its ability to reconstruct internal voids and hidden boundaries. The method involves using piezoelectric transducers to generate waves that propagate through a solid structure. These waves reflect off internal boundaries, such as cracks or other impurities, before being measured by a collection of sensors. FWI then employs an optimization process to iteratively update the material distribution of a computational model, thereby minimizing the difference between simulated and measured data. Numerical studies [12,13,14,15] and experimental validation [16] have demonstrated the enhanced tomographic capabilities of FWI compared to non-iterative approaches such as the total focusing method (TFM) [17] and reverse time migration (RTM) [18].

While classical methods using polynomial representations of material fields still dominate [11,19,20], recent studies in seismic exploration suggest that neural networks can improve the reconstruction quality by serving as efficient parameterizations for material parameter fields [21,22,23,24]. Similar benefits have been observed in structural optimization [1,2,3], where the parametrization has been coined neural topology optimization [4]. Herrmann et al. [25] extended these benefits to ultrasonic testing, demonstrating that neural network-based FWI can accurately detect internal voids. A key advantage of using neural networks as parametrizations is their ability to suppress high-frequency artifacts in the reconstructed images, thereby improving tomographic accuracy. This advantage arises because neural networks tend to learn smoother, low-frequency spatial patterns before capturing high-frequency details [26], which mitigates common artifacts in classical FWI.

However, a significant challenge remains: the performance of these neural network-based methods is highly sensitive to the initial weights of the network. To address this, researchers have proposed various advanced strategies. Vantassel et al. [27] used neural networks to generate informed initial guesses for classical FWI. Müller et al. [28] introduced a transfer learning approach that uses a network pre-trained on sensor data prior to their use in inversion. In [29], the network is pre-trained on sensor data and the corresponding damage and subsequently used for the FWI. The authors of the paper at hand further developed this idea in [30] by creating a transfer learning framework that uses adjoint gradients derived from sensor measurements to train a neural network to predict the underlying material field. This has the advantage of creating an image-to-image mapping with identical dimensions, which is agnostic to sensor and source placement. Further advances have been made by not relying on labeled data through meta-learning, in which multiple optimizations are performed with the goal of finding a better initial starting point [5].

The methods named above stand in stark contrast to purely data-driven approaches [31,32,33,34,35], which attempt to learn mappings from sensor data to material properties without incorporating physical constraints. Purely data-driven techniques cannot deliver reliable results due to the absence of an error measure on the predictions, as opposed to neural network-based optimization methods. Although such optimization methods are computationally expensive, pre-training can accelerate the minimization procedure.

A major limitation of all the aforementioned methods, whether purely data-driven or based on stringent minimization procedures constrained by the wave equation, is their reliance on grid-dependent neural network architectures, such as Convolutional Neural Networks (CNNs). These networks struggle to generalize effectively to complex, unstructured geometries that significantly deviate from their training data. This drawback impedes the application of neural network-based methods in non-academic scenarios.

To address this limitation, we propose utilizing Graph Convolutional Networks (GCNs) [36]. GCNs operate on graph-structured data (Figure 1), where nodes represent spatial points and edges encode geometric or physical relationships between them. In contrast to grid-based CNNs, GCNs inherently manage irregular topologies, which allows them to adapt to geometries starkly different from those they were trained on. This adaptability makes GCNs especially well-suited for transfer learning in Full Waveform Inversion (FWI), as they readily facilitate generalization across various geometries, especially those different from the training set.

Due to their inherent flexibility [37], GCNs have been widely adopted in domains requiring generalization over graph structures. They have been successfully applied to learning mesh-based physical simulations [6,38] and have also been integrated with Physics-Informed Neural Networks (PINNs) [39,40,41]. It is worth noting that PINNs are comparatively inefficient relative to classical methods, particularly in the context of FWI [7]. Readers seeking a comprehensive review of GCNs are directed to [42].

In this work, we propose using GCNs to parameterize the material field, an approach inspired by our prior work that utilized CNNs [25,29,30]. To the best of our knowledge, this is the first application of graph neural networks along with FWI. In particular, we demonstrate the following:

• We pre-train a GCN to learn to map an adjoint gradient, calculated from the first iteration of classical FWI, to the true material distribution. Both the input (the adjoint gradient) and the output (the material distribution) share the same underlying graph structure, enabling a simple network architecture.
• Using GCNs along with input and output data that are defined on identical graphs allows the method to be independent of specific experimental conditions, such as the number of sensors, the material’s properties, or the excitation frequency.
• Once the GCN is pre-trained, it generalizes to a wide range of setups and material characteristics (see Section 3). The network’s output then serves as the initial model for the FWI process. The GCN’s weights are subsequently fine-tuned during the FWI iterations to produce the final, refined damage distribution.

Figure 1 depicts the proposed approach.

The remainder of the paper is organized as follows: the next section presents the mathematical background of the proposed method. Section 3 details the results using three distinct case studies. We then discuss the broader implications and open questions in Section 4, and conclude in Section 5.

2. Methodology

2.1. Physical Models

In the current work, we use the scalar wave equation to generate inexpensive pre-training datasets on a rectangular geometry, while, in the inversion, wave propagation is modeled using the elastic wave equation. As proposed in [12], void- or crack-like defects can be described by a degradation in the material’s density. Introducing a dimensionless scaling function $\gamma \left(\mathit{x}\right)$, and assuming a constant density $\rho$ and P-wave speed $v_{\mathrm{p}}$ of the background material, the 2D scalar wave equation in the spatial domain $\mathsf{\Omega }\subset {\mathbb{R}}^2$ and for time $\mathcal{T}\subset \mathbb{R}$ can be written as

$$\gamma \left(\mathit{x}\right)\rho \ddot{u}(\mathit{x},t)-\nabla ·\left(\gamma \left(\mathit{x}\right)\rho v_{\mathrm{p}}^2\nabla u(\mathit{x},t)\right)=f(\mathit{x},t),\mathit{x}=[x,y]\in \mathsf{\Omega },t\in \mathcal{T},$$

(1)

where $u:\mathsf{\Omega }\times \mathcal{T}\to \mathbb{R}$ is the scalar wave field, $\ddot{u}$ the corresponding acceleration, and f the external force. Equivalently, the 3D elastic wave equation for $\mathsf{\Omega }\subset {\mathbb{R}}^3$ is given by

$$\gamma \left(\mathit{x}\right)\rho \ddot{\mathit{u}}(\mathit{x},t)-\nabla ·\left(\gamma \left(\mathit{x}\right)\mathbf{C}:\nabla \mathit{u}(\mathit{x},t)\right)=f(\mathit{x},t),\mathit{x}=[x,y,z]\in \mathsf{\Omega },t\in \mathcal{T},$$

(2)

where the elastic wave field $\mathit{u}:\mathsf{\Omega }\times \mathcal{T}\to {\mathbb{R}}^3$ has three components, and the constant linear isotropic elastic constitutive tensor $\mathbf{C}$ of undamaged material. Parametrized with the density $\rho$ and P- and S- wave speeds $v_{\mathrm{p}}$ and $v_{\mathrm{s}}$, $\mathbf{C}$ is written as

$$C_{ijkl}=\rho (v_{\mathrm{p}}^2-2v_{\mathrm{s}}^2){\delta }_{ij}{\delta }_{kl}+\rho v_{\mathrm{s}}^2{\delta }_{ik}{\delta }_{jl}+\rho v_{\mathrm{s}}^2{\delta }_{il}z,{\delta }_{jk}.$$

(3)

Additionally, homogeneous initial conditions are set at $t=0$, and homogeneous Neumann boundary conditions are prescribed on $\partial \mathsf{\Omega }$.

2.2. Classical FWI

In classical FWI, the material field is commonly represented by a set of coefficients defined on the underlying computational grid or mesh. In this work, the vector $\widehat{\gamma }$ contains the scaling function values at the vertices of a finite element mesh. We seek to find an optimal vector of coefficients ${\widehat{\gamma }}^{*}$ that minimizes the $L_2$ misfit function $\mathcal{L}\left(\widehat{\gamma }\right)$:

$${\widehat{\gamma }}^{*}=\underset{\widehat{\gamma }}{\mathrm{argmin}}\mathcal{L}\left(\widehat{\gamma }\right).$$

(4)

The $L_2$ error is summed over all experiments and receiver positions. For the elastic wave equation, we obtain

$$\mathcal{L}(\widehat{\gamma };{\mathit{u}}^{\mathrm{o}})={\displaystyle \frac{1}{2}}{\int }_T{\int }_{\mathsf{\Omega }}\sum _{s=1}^{N_{\mathrm{s}}}\sum _{r=1}^{N_{\mathrm{r}}}\left[\left({\widehat{\mathit{u}}}_s(\widehat{\gamma };\mathit{x},t)-{\mathit{u}}_s^{\mathrm{o}}({\mathit{x}}^r,t)\right)·{\mathit{n}}^r{)}^2\delta (\mathit{x}-{\mathit{x}}^r)\right]\mathrm{d}\mathsf{\Omega }\mathrm{d}t,$$

(5)

in which ${\widehat{\mathit{u}}}_{\mathrm{s}}(\widehat{\gamma };x,t)$ represents the simulated wave field based on the coefficients $\widehat{\gamma }$, while ${\mathit{u}}_s^{\mathrm{o}}(x^r,t)$ represents the sparse measurements from the target domain (referred to as the “observed wavefield” in the subsequent text), and ${\mathit{n}}^r$ is a unit vector pointing in the measurement direction of the receiver. The complete set of measurements across all $N_{\mathrm{s}}$ experiments is denoted as ${\mathit{u}}^{\mathrm{o}}={\left\{{\mathit{u}}_s^{\mathrm{o}}\right\}}_{s=1}^{N_{\mathrm{s}}}$ and $N_{\mathrm{r}}$ is the number of receivers.

To calculate how the cost function changes with respect to material coefficients, we employ the adjoint method. For those interested in the complete mathematical derivation, see, e.g. [12]. The adjoint method allows us to compute the derivative of the misfit functional with respect to a perturbation of the density scaling function $\gamma \left(\mathit{x}\right)$ in the direction $\delta \gamma$

$${\nabla }_{\gamma }\mathcal{L}\delta \gamma ={\int }_{\mathsf{\Omega }}{K}_{\gamma }\left(\mathit{x}\right)\delta \gamma \mathrm{d}\mathsf{\Omega }.$$

(6)

For the elastic wave equation, the sensitivity kernel $K_{\gamma }$ is

$$K_{\gamma }\left(\mathrm{x}\right)={\int }_{\mathcal{T}}-\dot{\mathbf{u}}·{\mathbf{u}}^{†}+\nabla \mathbf{u}:\mathbf{C}:\nabla {\mathbf{u}}^{†}\mathrm{d}t,$$

(7)

where ${\mathbf{u}}^{†}$ is the solution of the adjoint wave equation. The gradient from all the vertex points can be directly applied in gradient-based optimization algorithms to find the optimal material coefficients. For the scalar wave equation, corresponding equations can be formulated, e.g., see [12]. In this work, the commercial software Salvus from Mondaic [43] has been used as the forward solver and for calculating the adjoint gradients.

2.3. Graph Convolutional Network-Based Transfer Learning FWI

2.3.1. Graph Convolutional Networks

GCNs are structured around the principle of message passing [8], where each node aggregates information from the nodes it is connected with to refine its representations over the layers of the network. The aggregation operations in GCNs used in the current work are based on vanilla GCNs [44], where node features are updated using weighted combinations of adjacent node attributes,

$$x_i^{\prime }=W_1{x}_i+\sum _{j\in {\mathcal{N}}_i}{W}_2·f\left(x_j\right),$$

(8)

where $W_1$ and $W_2$ are the trainable parameters of the layer, $x_i$ is the value of the node from the previous layer, and $x_j$ is the value of the neighboring nodes. All the connected neighboring nodes share the same weight matrix $W_2$, which is the key feature that allows the GCN to generalize to nodes with different numbers of neighbors, and thus to different graphs. Figure 2 shows the representation of a finite element mesh as a graph.

A key advantage of GCNs lies in their ability to process mesh-based data, which is inherently structured as a graph. Traditional grid-based methods struggle with irregular geometries, whereas GCNs are designed to operate on these structures. The PyTorch geometric library [45] is used, because it offers an excellent out-of-the-box implementation of many types of GCNs.

2.3.2. Pre-Training

The graph-based neural network discussed above is pre-trained using the adjoint gradients of the first iteration of a classical FWI. The same training data as used in our previous work [30], expressed as a graph, is used to pre-train the GCN. The dataset is generated by varying the semi-major and minor axes a and b, the center $x_c$ and $y_c$, and the rotation $\varphi$ of a single elliptical damage. The generated damaged density distributions ${\gamma }_k^{\mathrm{true}}$ serve as input to the wave equation to compute the corresponding wavefields $\mathit{u}$. These wavefields, evaluated at the sensor locations, act as the observed wavefields ${\mathit{u}}^{\mathrm{o}}$. Assuming an undamaged density distribution ${\widehat{\gamma }}_0$, the adjoint gradient ${\nabla }_{\gamma }\mathcal{L}({\widehat{\gamma }}_0;{\mathit{u}}^{\mathrm{o}})$ is calculated using the observed wavefields. These adjoint gradients ${\nabla }_{\gamma }\mathcal{L}({\widehat{\gamma }}_0;{\mathit{u}}^{\mathrm{o}})$ form the input to the GCN. The damaged density distributions (with ellipsoidal voids) are used to generate the observed wavefield from the target output ${\gamma }_k^{\mathrm{true}}$. For training the GCN, the input and the output data are converted into a graph as depicted in Figure 2. The features of the mesh nodes become the features of the nodes processed by the GCN. The GCN is trained in a supervised manner to minimize the standard mean-squared error loss function

$${\mathcal{L}}_D={\displaystyle \frac{1}{N}}\sum _{k=1}^N\left(\parallel {\gamma }_k^{\mathrm{true}}-{\widehat{\gamma }}_k{\parallel }_2^2\right)$$

(9)

between predicted ${\widehat{\gamma }}_k$ and ground truth target values ${\gamma }_k^{\mathrm{true}}$.

Neural network architecture:

The GCN has a symmetric encoder-decoder architecture designed for node-level feature transformation. The model comprises an encoder, a bottleneck layer, and a decoder, with each component constructed from a series of GCN layers, followed by Graph Normalization layer1 for feature normalization and a PReLU activation (except for the final output layer). The encoder progressively transforms the input node features through the layers, resulting in a latent embedding space. A single bottleneck layer connects the encoder and decoder. Each layer maintains the dimensionality of the output from the previous layer. The decoder reconstructs node features from the bottleneck’s output. It mirrors the encoder’s structure, with its initial layers using CuGraphSAGEConv2 with max aggregation. At the end of the decoder layer, there is a layer of SAGEConv with a Multi-Aggregation module, which combines the outputs of three parallel Softmax Aggregation modules. Each Softmax Aggregation module $k\in \{1,2,3\}$ is parameterized by a learnable temperature $t_k$, initialized to $t_1=0.01$, $t_2=1$, and $t_3=10$ respectively ( Table 4). For each module k, an aggregated feature vector ${\mathbf{h}}_{S_k}$ is computed element-wise, where the $j^{\mathrm{th}}$ dimension ${\left({\mathbf{h}}_{S_k}\right)}_j$ is given by

$${\left({\mathbf{h}}_{S_k}\right)}_j=\sum _{u\in \mathcal{N}\left(v\right)}\left({\displaystyle \frac{exp\left(t_k·{\mathbf{m}}_{u,j}\right)}{{\sum }_{p\in \mathcal{N}\left(v\right)}exp\left(t_k·{\mathbf{m}}_{p,j}\right)}}\right)·{\mathbf{m}}_{u,j}.$$

Here, ${\mathbf{m}}_{u,j}$ denotes the $j^{\mathrm{th}}$ feature of message ${\mathbf{m}}_u$. The parameter $t_k$ is learnable, allowing the model to adaptively fine-tune the aggregation behavior of each neighbor. Finally, three aggregated vectors are combined element-wise using the LogSumExp function of PyTorch. The final aggregated message ${\mathbf{H}}_{\mathrm{final}}$ for node v, with its $j^{\mathrm{th}}$ dimension ${\left({\mathbf{H}}_{\mathrm{final}}\right)}_j$, is computed by

$${\left({\mathbf{H}}_{\mathrm{final}}\right)}_j=log\left(exp\left({\left({\mathbf{h}}_{S_1}\right)}_j\right)+exp\left({\left({\mathbf{h}}_{S_2}\right)}_j\right)+exp\left({\left({\mathbf{h}}_{S_3}\right)}_j\right)\right).$$

This LogSumExp combination serves as a smooth, differentiable approximation of the maximum function, enabling the GCN to dynamically emphasize the most important aggregation terms for each feature. Furthermore, the smooth nature of the function also allows better gradient flow during training as compared to other commonly used aggregation functions like max or min. The final layer of the GCN is essential since it improves the expressivity of the network. This, in turn, improves the convergence during the downstream task of TL-FWI. Adding more layers did not significantly improve the results.

Since GCNs and CNNs operate on fundamentally different data structures, there is a significant difference in their training speeds. GCNs are designed to handle irregular data structures such as graphs, where nodes can have a varying number of connections. This lack of a uniform structure means the GPU cannot apply the same highly parallel, uniform operations as it does for CNNs. Given the computational challenges, the GCN is trained on only 200 training cases as opposed to the 800 cases in our previous work [30]. This training data is based on the scalar wave equation on a rectangular 2D geometry, with a 70:30 train-test split. The corresponding training time for 200 epochs is $\sim 65$ minutes on a Nvidia Quadro RTX 8000. The input, i.e., the initial adjoint gradient, is converted into graph data using the PyTorch Geometric library. Here, the gradient value at each node of the mesh is stored as a 1D vector, and another 2D vector is constructed to represent the connectivity between all nodes. For constructing the connectivity, the radius threshold method was used. This forms the input to the GCN. The density distribution is analogously converted into a vector (in the same order as the input), which yields the output. Further details can be inferred from the code published in [46].

Pre-processing adjoint gradient:

During pre-training, the model was trained on data that is in a uniform, rectangular domain where the wavefields provided consistent illumination. This means the energy from the sources and sensors is spread evenly throughout the domain. As a result, the distortions in the gradients caused by distortions in the density due to the damage (damage fingerprint) were strong and easily detectable by the neural network. The addition of a region of interest (ROI) further refined this process by excluding near-field artifacts from the sources and sensors. In the downstream task, we tested irregular geometries and used a simulation time that was not long enough to allow the wavefields to illuminate the domain uniformly. This lack of uniform illumination caused the damage fingerprints to be weaker and difficult to detect, making them almost “invisible” to the neural network. By scaling the adjoint gradients using the wavefield, we effectively normalized the data. This resulted in the damage fingerprints being easily recognized by the neural network, resulting in better initial guesses. The scaling factor is determined as

$${\mathbf{u}}_{\mathrm{norm},i}={\left(\sum _{t=1}^T\left|u_{t,x,i}\right|\right)}^2+{\left(\sum _{t=1}^T\left|u_{t,y,i}\right|\right)}^2+{\left(\sum _{t=1}^T\left|u_{t,z,i}\right|\right)}^2,$$

(10)

where u is the displacement field for all node points along the x, y, and z directions, t is the time steps used for checkpointing, and i is the node. The normalized gradient is calculated by scaling the quantity

$${\nabla }_{\widehat{{\gamma }_{\mathit{i}}}}\mathcal{L}{(\widehat{\gamma };{\mathit{u}}^{\mathrm{o}})}_{\mathrm{norm}}={\displaystyle \frac{{\nabla }_{\widehat{{\gamma }_{\mathit{i}}}}\mathcal{L}(\widehat{\gamma };{\mathit{u}}^{\mathrm{o}})}{u_{\mathrm{norm},i}}},$$

(11)

where ${\nabla }_{\widehat{{\gamma }_{\mathit{i}}}}\mathcal{L}(\widehat{\gamma };{\mathit{u}}^{\mathrm{o}})$ is the adjoint gradient from the first iteration of the classical FWI. This process amplifies the weak gradients from the damaged boundaries, making them more prominent and easier for the neural network to identify. This normalization directly led to a significant improvement in the quality of the initial guesses, as the network could now reliably detect and localize the damage despite the non-uniform illumination. Figure 3 shows the gradient pre-processing steps.

Generalization to arbitrary shape and experiment configuration:

Using an adjoint gradient for pre-training the GCN instead of raw wavefields has several advantages. Most importantly, the graph network learns to map the gradient to the density scaling function on the same graph. This renders the architecture of the GCN to be independent of the number of sensors, the time steps of the raw wavefield, and the type of wave equation used (scalar or elastic) as opposed to other transfer learning approaches [28,29]. This flexible formulation constitutes a key ingredient of the paper at hand, as it enables automatic generalization to other geometries, sensor configurations, excitation locations, source types, and types of wave equations. Consequently, the data necessary to pre-train the GCN is generated cheaply, using simple rectangular geometry simulated with the scalar wave equation. The generalization capabilities are strong enough that the pre-trained GCN can be applied to FWI based on the elastic wave equation for geometrically and topologically more complex scenarios in 3D. This is remarkable, as the training data only covers a single 2D rectangular geometry on the scalar wave equation.

2.3.3. Graph Convolutional Network-Based Transfer Learning (GCN-based TL-FWI)

In transfer learning-based FWI, as proposed in [30], a CNN is used for discretization of the density scaling function $\gamma$. The CNN is pre-trained to identify the damage directly from the adjoint gradient calculated from the first iteration of a classical adjoint method. Using the adjoint gradient instead of the raw signal allows us to simplify the network architecture, because the input and output shapes are identical. Once trained, the network is used for the downstream task of FWI. In the paper at hand, a further step is taken towards generalization by using a GCN to parameterize the density scaling function $\gamma$. The developed method is termed GCN-based TL-FWI (Figure 4) in the remainder of this work. The gradients of the cost function with respect to the weights of the pre-trained GCN are calculated using the chain rule [25]. The pseudo code of the GCN-based TL-FWI algorithm is also shown in Algorithm 1.

2.4. Classical FWI with Initial Guess

Classical FWI can similarly be accelerated using the pre-trained GCNs to provide an initial guess of the damage. The initial guess serves as the starting model for a classical FWI method, where the material field is defined by a set of coefficients defined on the underlying computational grid or mesh. In this work, the L-BFGS optimization method from the PyTorch library was used to solve the optimization task. The pseudo code for classical FWI algorithm is shown in 2. Comparing the GCN-based TL-FWI with a classical FWI (provided with the same initial guess) isolates the regularizing effect of the GCN on the FWI.

3. Results

3.1. Bridge to Previous Work

In order to set a point of reference of established work, the performance of a pre-trained GCN is compared to that of a CNN from [30], where a scalar wave equation has been used for carrying out FWI.

Figure 5 shows that the CNN provides a better initial guess as compared to the GCN. Additionally, the CNN recovers all voids with higher accuracy and slightly faster in the inversion than the GCN. Both the CNN-based and GCN-based TL-FWI require fewer than 15 epochs to reconstruct the damage shapes. The reconstruction of the CNN-based TL-FWI contains fewer artifacts in the undamaged domain, and the resolved damage shape is very sharp. GCN-based TL-FWI requires some additional epochs to achieve the same level of accuracy. There are various reasons for the better performance of a CNN as compared to a GNN. Section 4 contains a detailed discussion about the differences between a CNN and a GCN. The drawbacks of GCN are largely offset by its generalization capabilities, as will be demonstrated in the sequel.

Note that the double voids depicted in Figure 5 were also outside the training dataset; the present case depicts the performance of the GCN on an arbitrary damage scenario. While CNN-based FWI may achieve faster convergence in the simple rectangular scenario, GCNs offer greater flexibility for handling more complex and unstructured geometries, which is crucial in cases where the investigated domain deviates significantly from a simple rectangle.

3.2. Generalization to 2D T-Shaped Geometry and Elastic Wave Equation

The generalization capability of the GCN-based TL-FWI is presented step by step, starting from the simplest case of a T-shaped 2D geometry, which presents a moderate generalization challenge. As previously elaborated, GCNs allow generalizing to arbitrary meshes by converting the mesh (and the corresponding quantities defined on the mesh nodes) into a graph. The details can be found in Section 2.3.2. The performance of GCN-based TL-FWI is compared with classical FWI with and without an initial guess from the pre-trained GCN (as described in Section 2.4).

3.2.1. Experimental Configuration

Moving from the initial rectangular 2D training domain (which relied on the scalar wave equation), first, we assess the generalization to a T-shaped domain and to the elastic wave equation. Furthermore, different material properties and excitation frequencies are utilized. The T-shaped geometry and the corresponding mesh were generated using the Gmsh library in Python. It was then converted into Exodus format using the Pyexodus library and read into the Salvus [47] solver from the software framework published by Mondaic AG in Python. The pre-trained GCN outputs the density values at the corresponding vertices of the mesh. For the simulation of elastic wave propagation, point sources were utilized. Unlike the scalar case, which typically involves a single force component or pressure source, the elastic simulations utilize vector point sources capable of generating both horizontal and vertical force components. The source time function for these simulations is a Ricker wavelet. The details of the simulation and the material properties are provided in Table 1.

The initial guess is also used as the starting model for the classical FWI. L-BFGS was used for the optimization of the classical method, whereas Adam [9] is utilized for the GCN-based TL-FWI. As the L-BFGS optimizer relies on multiple function evaluations per epoch, a comparison based on the number of epochs is unfair. Instead, the results are compared at specific times. Overall, the 10 iterations of L-BFGS optimization are approximately equivalent to the 75 iterations of Adam in terms of runtime.

3.2.2. Results

Good initial guess:

The GCN-based TL-FWI is used for detecting a simple scenario consisting of a single damage. The initial guess in Figure 6 (a) shows that the network is able to correctly guess the location of the damage. As a reference, classical FWI with the initial guess is also used to detect the same damage. The convergence results are shown in Figure 6. It shows that the TL-FWI method is much faster than the classical methods. At $405seconds$, the TL-FWI method has already recovered the complete damage shape with negligible artifacts, while the classical FWI is unable to recover the shape accurately, even after $1240seconds$.3. This highlights the capability of the proposed GCN-based TL-FWI method in reconstructing the damage shape accurately and faster than the classical FWI method.

Bad initial guess:

The following case contains three voids presenting a greater challenge. The pre-trained GCN only predicts a small portion of the damage closest to the sensors (and is unable to predict the other two), see Figure 7 (b). Furthermore, incorrect voids are predicted in the bottom left corner.

During the optimization cycles of FWI, the GCN-based TL-FWI is able to correct the bad initial guess, see Figure 7 (c). The mean squared error of the final reconstruction is $\approx 26000$. The classical FWI with the initial guess (Figure 7 (d)) is also able to recover the voids. However, the recovered shapes are not as accurate (mean squared error $\approx 109,000$). Furthermore, the classical FWI required $\approx 3400seconds$ as compared to GCN-based TLFWI, which required $\approx 2400seconds$. In particular, this demonstrates that classical FWI suffers from substantial noise even if the same prediction is given as a starting point for the iteration, as in the case of the GCN, as well as requiring more time. Therefore, it can be concluded that GCN-based TL-FWI still works the best at reconstructing damage on a T-shaped geometry despite a bad initial guess.

3.3. Generalization to 3D Geometry: Flange

The next numerical example demonstrates the ability of the GCN-based TL-FWI to generalize to a 3D flange geometry, which is a common mechanical component. The focus is on the prediction of the pre-trained GCN, which has only been trained on 2D rectangular domains.

Table 2. Properties of the simulation for Flange

Property	Value
Excitation frequency	$3\times {10}^5\mathrm{Hz}$
Simulation time	$2\times {10}^{-5}$ s
$\rho$	$7850{\displaystyle \frac{\mathrm{K}\mathrm{g}}{{\mathrm{m}}^3}}$
$v_{\mathrm{p}}$	$5900{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}$
$v_{\mathrm{s}}$	$3200{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}$

Table 2. Properties of the simulation for Flange

Property	Value
Excitation frequency	$3\times {10}^5\mathrm{Hz}$
Simulation time	$2\times {10}^{-5}$ s
$\rho$	$7850{\displaystyle \frac{\mathrm{K}\mathrm{g}}{{\mathrm{m}}^3}}$
$v_{\mathrm{p}}$	$5900{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}$
$v_{\mathrm{s}}$	$3200{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}$

3.3.1. Numerical Configuration

The flange, with an outer radius of $0.3$ m and an inner radius of $0.15$ m, has a thickness of $0.005$ m. As shown in Figure 8(a), it contains five known holes and two ellipsoidal internal damages located near them. For data acquisition, five sensors and 25 receivers are positioned on the outer radius and at half the thickness of the flange. The GCN, which is pre-trained on a 2D rectangle (Figure 1), is now applied to a Flange’s complex geometry. Figure 8(b) shows the initial gradient, which is the input to the pre-trained GCN. The initial gradient was pre-processed in the same way as described in the previous section (Section 2.3.2.2).

3.3.2. Results

Figure 8(a) shows two damaged areas located near holes in a flange. The initial prediction provided by our GCN, shown in Figure 8(b), accurately predicts the location of these damaged regions. However, it also incorrectly identifies a damaged area at the bottom of the flange. This false positive is likely due to strong reflections from the hole boundaries, which can produce high gradients that the GCN interprets as damage.

Classical FWI successfully recovers both damaged areas, but with significant artifacts, particularly in the vicinity of the damage(Figure 8(d)). The recovered shapes are generally satisfactory, but the widespread artifacts are a common issue with classical FWI (mean-squared error $219,000$). The runtime of classical FWI was approximately 13.7 hours. The proposed GCN-based TL-FWI (Figure 8(c)) outperforms the classical method as it recovers the damaged areas with higher accuracy (mean-squared error $83,000$) and substantially less noise in about half the time (5.5 hours).

3.4. Generalization to 3D Geometry: Cylinderical Column

This numerical example demonstrates the ability of the GCN-based TL-FWI to generalize to another 3D geometry of a support column, a geometry that is quite common in civil engineering.

3.4.1. Numerical Configuration

The column has a radius of $1.25$ m in the center, a radius of $1.875$ m at the base, and a length of $3.75$ m (Figure 9). There are four sources and 32 receivers placed on the outer surface of the column at a height of $0.75$ m and $3.25$ m from the base. A significant challenge is presented by the stark difference in numerical configurations; the new geometry is considerably larger, and the excitation frequency is two orders of magnitude lower. The GCN, which is pre-trained on the 2D rectangle as shown in Figure 1 with much higher excitation frequency, is now applied to the column for reconstructing the scaling function (Figure 9(a)). The initial gradient was preprocessed in the same way as described in Section 2.3.2.2 before being input to the pre-trained GCN.

Table 3. Properties of the simulation for cylindrical column

Property	Value
Excitation frequency	$2000\mathrm{Hz}$
Simulation time	$2\times {10}^{-3}\mathrm{s}$
Density	$2400{\displaystyle \frac{\mathrm{K}\mathrm{g}}{{\mathrm{m}}^3}}$
$v_{\mathrm{p}}$	$3500{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}$
$v_{\mathrm{s}}$	$2000{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}$

Table 3. Properties of the simulation for cylindrical column

Property	Value
Excitation frequency	$2000\mathrm{Hz}$
Simulation time	$2\times {10}^{-3}\mathrm{s}$
Density	$2400{\displaystyle \frac{\mathrm{K}\mathrm{g}}{{\mathrm{m}}^3}}$
$v_{\mathrm{p}}$	$3500{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}$
$v_{\mathrm{s}}$	$2000{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}$

3.4.2. Results

Despite the substantial differences in the configuration of the pre-training, the GCN successfully predicts the correct location of the damage (Figure 9(b)) in the computation. While the initial prediction of the damage shape was found to be imperfect, its accurate localization is still very beneficial in accelerating the GCN-based TL-FWI. This preliminary damage estimate is then used as a starting model for the subsequent GCN-based TL-FWI process.

The final reconstructed scaling function obtained via GCN-based TL-FWI in (Figure 9(c)) and classical FWI with an initial guess is shown in (Figure 9(d)). The GCN-based TL-FWI method results in the most accurate scaling function reconstruction, closely matching the actual shape of the defect. The favorable difference in terms of noise becomes evident in the final reconstruction when comparing to the results obtained by classical FWI (Figure 9(d)). Although the computational time for both methods is comparable in this case, the mean-squared error (MSE) for GCN-based TL-FWI is $\approx 21,000$, which is significantly lower than the $\approx 46,600$ recorded for classical FWI with the initial guess. This difference demonstrates the GCN-based TL-FWI’s capability to generalize well and effectively suppress artifacts, even within complex 3D geometries.

4. Discussion

The proposed GCN-based TL-FWI is able to accurately predict the damage even in 3D cases after being trained on data generated from the 2D rectangle depicted in Figure 1. The GCN-based TL-FWI possesses a much stronger generalization in terms of geometry, material properties, sensor placement, and excitation frequencies as compared to a CNN as used in [30]. Even in cases when the adjoint gradient from the first iteration does not contain information to allow the GCN to provide a good initial guess, the pretrained GCN is still able to recover the damage much faster than classical FWI, given the identical initial guess.

However, in comparison to the U-net [30], GCNs are more difficult and slower to train, but offer much better generalization properties. The difference in the performance between CNNs and GCNs is attributed to their fundamentally distinct architectural designs. Structured, grid-based data is processed by CNNs through the application of fixed-size convolutional kernels. The weights in a CNN are specific to each position, enabling the detection of spatially localized patterns. While this approach is effective for uniform grids, such as those found in images, its adaptability to irregular or varying geometries is limited, as the learned weights are inherently tied to the spatial arrangement of the training data.

By contrast, a graph-based representation is processed by GCNs, wherein spatial points are represented as nodes and their connectivity is encoded as edges. For each new geometry, a graph must be constructed to reflect the specific topological relationships. This fundamentally differentiates GCNs from CNNs. This graph construction introduces a minor additional overhead, as the node-edge relationships suitable for each geometry must be defined.

Information from neighboring nodes is aggregated in GCNs using a shared, trainable weight, with the receptive field determined by the graph’s connectivity and number of layers. This adaptability is critical in TL-FWI, where the density field is parameterized by the GCN, and damage distributions are predicted across structures of various scales and experimental configurations.

Influence of the number of neighbors:

The number of neighboring nodes considered in a GCN’s aggregation process is a pivotal hyperparameter that governs the trade-off between reconstruction accuracy and computational feasibility. A larger number of neighbors expands the receptive field, providing richer contextual information that helps in noise suppression and improves the accuracy of the damage reconstruction using TL-FWI. As illustrated in Figure 10, increasing the number of neighbors accelerates convergence of FWI and reduces the mean squared error (MSE) by capturing wider spatial dependencies. However, this benefit is also accompanied by the risk of oversmoothing, where excessive neighbor aggregation dilutes fine-scale damage features by incorporating irrelevant information from distant nodes.

Limitations of GCN-based TL-FWI:

Memory constraints further complicate the trade-off between the number of neighbors and TL-FWI convergence. The memory cost of GCNs scales linearly with respect to the number of nodes and the number of neighbors, as each node’s feature update requires aggregating high-dimensional representations. Therefore, selecting an optimal number of neighbors presents a trade-off between achieving the desired predictive accuracy and speed, and satisfying the computational and memory requirements of the system. Graph connectivity is a hyperparameter that can be adjusted based on available computational resources without altering the network architecture. This is a key consideration that sets GCNs apart from CNNs. Furthermore, the memory requirement also increases with the number of nodes in the finite element mesh. In the presented case studies, the peak memory requirement reached the maximum capacity of the GPU (48 GB) for the flange, which had roughly $2\times {10}^6$ nodes. In other words, the increase in the speed that is possible due to GCN-based TL-FWI comes at the cost of increased peak memory usage.

However, various methods have already been proposed to tackle the high memory requirements of big GCNs, such as using mixed precision, where the forward and backward pass of the neural network is carried out in half precision instead of full precision, which halves the memory requirement. Another method exploits subgraph sub-sampling [48,49,50] where instead of training the whole graph, small graphs are sampled from the large graph and trained in mini-batches. Gradient checkpointing [51] can be used as well. Instead of storing all intermediate activations of the network, only a selected few are saved, and during backpropagation, the unsaved activations are recomputed on-the-fly. All these recent developments promise a reduction of the high memory requirement of GCN-based TL-FWI.

5. Conclusions

This paper improves the Transfer-Learning Full Wave Form Inversion (TL-FWI) method proposed in [30] by replacing the CNN with a GCN. The advantage of using a GCN is its ability to generalize to any geometry. Additionally, pre-training is applied [30] where the adjoint gradient computed from the first iteration of classical FWI is used as an input to train the GCN to guess the corresponding damage density distribution $\gamma$. Using pre-training along with a GCN enables the presented method to be agnostic to geometry, experimental configuration (such as excitation frequency and number of sensors), material properties, and even specific wave physics. As demonstrated by various numerical case studies, the proposed method enables training the GCN on data generated from a 2D rectangle, as depicted in Figure 1, using the scalar wave equation. The pre-trained network is then able to generalize to 3D geometries of entirely different scales, with varying experimental configurations, different material properties, and applied to the elastic wave equation instead of the scalar wave equation. This is of great advantage, since using CNNs, as in [30], requires rectangular domains and application-specific pre-training, which is too costly in engineering-relevant scenarios.

The proposed GCN-based TL-FWI provides a good initial guess based on the adjoint gradient. Furthermore, GCN-based TL-FWI is able to reconstruct the damage faster than classical FWI and, at the same time, exhibits fewer artifacts.

The primary drawback of the proposed method is its high memory requirements. For future work, new methods will be explored to optimize the GCN so that larger geometries can be tackled and data from real-world experiments may be applied.

GCN Architecture

The GCN has an encoder-decoder architecture. The graph nodes are kept the same, and the number of channels is increased from one 256 till the bottleneck layer and subsequently reduced back to one in the final layer. The encoder part of the network is frozen. It was observed that unfreezing the weights did not improve the overall performance.

Table 4. Graph convolutional NN architecture has an encoder-decoder type architecture. The network consists of an encoder layer, a bottleneck layer, and a decoder layer. The total number of parameters is $571,281$. However, the encoder layer’s weights are frozen, which results in 351180 trainable parameters.

Layer	Shape after layer	Learnable parameters
Input	$N\times 1$	0
Encoder layer
Graph Convolution with Max Aggr & GN & PReLU	$N\times 32$	$96+96+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 64$	$4160+192+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 128$	$16512+384+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 256$	$65792+768+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 256$	$131328+768+1$
Bottleneck layer
Graph Convolution with Max Aggr & GN & PReLU	$N\times 256$	$131328+768+1$
Decoder layer
Graph Convolution with Max Aggr & GN & PReLU	$N\times 256$	$131328+768+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 128$	$65664+384+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 64$	$16448+192+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 32$	$4128+96+1$
Graph Convolution with Multi Aggr & GN & Sigmoid	$N\times 1$	$68+3+1$

Table 4. Graph convolutional NN architecture has an encoder-decoder type architecture. The network consists of an encoder layer, a bottleneck layer, and a decoder layer. The total number of parameters is $571,281$. However, the encoder layer’s weights are frozen, which results in 351180 trainable parameters.

Layer	Shape after layer	Learnable parameters
Input	$N\times 1$	0
Encoder layer
Graph Convolution with Max Aggr & GN & PReLU	$N\times 32$	$96+96+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 64$	$4160+192+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 128$	$16512+384+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 256$	$65792+768+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 256$	$131328+768+1$
Bottleneck layer
Graph Convolution with Max Aggr & GN & PReLU	$N\times 256$	$131328+768+1$
Decoder layer
Graph Convolution with Max Aggr & GN & PReLU	$N\times 256$	$131328+768+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 128$	$65664+384+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 64$	$16448+192+1$
Graph Convolution with Max Aggr & GN & PReLU	$N\times 32$	$4128+96+1$
Graph Convolution with Multi Aggr & GN & Sigmoid	$N\times 1$	$68+3+1$

GCN-Based TLFWI Pseudocode

Algorithm 1: GCN-based Transfer Learning FWI

1: Input:Background Mesh $\mathit{M}$, Pre-trained GCN Model $\mathcal{G}$, Observed Data ${\mathit{u}}_{obs}$, FWI Hyperparameters ${\mathrm{P}}_{\mathrm{tlfwi}}$   2: Output: Inverted Density Model ${\rho }_{\mathrm{ver}}$   3: Function InitialGradient($\mathit{M},{\mathit{u}}_{\mathrm{obs}}$)   4:    Compute initial forward wavefield ${\mathit{u}}_{\mathrm{init}}$ on $\mathit{M}$ with ${\rho }_{\mathrm{bg}}$   5:    Compute initial adjoint source ${\mathit{s}}_{\mathrm{adj}}\leftarrow \mathrm{L}2\mathrm{Misfit}({\mathit{u}}_{\mathrm{init}},{\mathit{u}}_{\mathrm{obs}})$   6:    Compute initial gradient $\nabla {\mathcal{J}}_{\mathrm{elm}}\leftarrow \mathrm{SalvusAdjoint}({\rho }_{\mathrm{bg}},{\mathit{s}}_{\mathrm{adj}})$   7:    Compute wavefield energy ${\mathit{u}}_{\mathrm{norm}}\leftarrow {\sum }_t{\sum }_i{\parallel {\mathit{u}}_{\mathrm{init}}\parallel }^2$   8:    Clip and scale the initial gradient ${(\nabla {\mathcal{J}}_{\mathrm{elm}})}_{\mathrm{norm}}\leftarrow \mathrm{ScaleAndNormalize}(\nabla {\mathcal{J}}_{\mathrm{elm}},{\mathit{u}}_{norm})$   9:    return ${(\nabla {\mathcal{J}}_{\mathrm{elm}})}_{\mathrm{norm}}$ 10: End Function 11: $\nabla {\mathcal{J}}_{\mathrm{init}}\leftarrow \mathrm{InitialGradient}(\mathit{M},{\mathit{u}}_{\mathrm{obs}})$ 12: Map initial gradient to graph vertices ${(\nabla {\mathcal{J}}_{\mathrm{init}})}_{\mathrm{ver}}\leftarrow \mathrm{MeshToGraph}(\nabla {\mathcal{J}}_{\mathrm{init}})$ 13: Generate the graph connectivities $\mathit{E}$ 14: Load the pre-trained GCN network ${\theta }_t\leftarrow$ LoadWeights 15: $\mathrm{optimizer}\leftarrow \mathrm{Adam}(\left[{\theta }_{\mathrm{t}}\right],\mathrm{lr})$ 16: for$\mathrm{it}=1$toepochsdo 17:     $\mathrm{optimizer}.\mathrm{zero\_grad}\left(\right)$ 18:     Initialize ${\rho }_{ver}$ 19:     GCN prediction at mesh vertices ${\rho }_{\mathrm{ver}}\leftarrow \mathcal{G}({(\nabla {\mathcal{J}}_{\mathrm{init}})}_{\mathrm{ve}r},\mathit{E})*{\rho }_{\mathrm{bg}}$ 20:     Map density values to Elements ${\rho }_{\mathrm{elm}}\leftarrow \mathrm{GraphToMesh}\left({\rho }_{\mathrm{ver}}\right)$ 21:     Update mesh $\mathit{M}$ with current ${\rho }_{\mathrm{elm}}$ model 22:     ${\mathit{u}}_{\mathrm{syn}}\leftarrow \mathrm{SalvusForward}\left(\mathit{M}\right)$ 23:     $\mathcal{L}\leftarrow \mathrm{L}2\mathrm{Misfit}({\mathit{u}}_{\mathrm{syn}},{\mathit{u}}_{\mathrm{obs}})$ 24:     ${\mathit{s}}_{\mathrm{adj}}\leftarrow \mathrm{GenerateAdjointSource}\left(\mathcal{L}\right)$ 25:     $\nabla {\mathcal{J}}_{\mathrm{elm}}\leftarrow \mathrm{SalvusAdjoint}({\rho }_{\mathrm{elm}},{\mathit{s}}_{\mathrm{adj}})$ 26:     $\nabla {\mathcal{J}}_{\mathrm{elm}}\leftarrow \mathrm{Clip}(\nabla {\mathcal{J}}_{\mathrm{elm}},\text{clip\_fraction})$ 27:     Map gradient values to graph vertices $\nabla {\mathcal{J}}_{\mathrm{ver}}\leftarrow \mathrm{MeshToGraph}(\nabla {\mathcal{J}}_{\mathrm{elm}})$ 28:     Calculate the gradient of the weights ${\rho }_{\mathrm{ver}}.\mathrm{backward}(\nabla {\mathcal{J}}_{\mathrm{ver}})$ 29:     Update the weights optimizer.step() 30: end for

Classical FWI with Initial Guess

Algorithm 2: Classical FWI with initial guess FWI

1: Input:Background Mesh $\mathit{M}$, Pre-trained GCN Model $\mathcal{G}$, Observed Data ${\mathit{u}}_{obs}$, FWI Hyperparameters ${\mathrm{P}}_{\mathrm{tlfwi}}$   2: Output: Inverted Density Model ${\rho }_{\mathrm{ver}}$   3: Function InitialGradient($\mathit{M},{\mathit{u}}_{\mathrm{obs}}$)   4:    Compute initial forward wavefield ${\mathit{u}}_{\mathrm{init}}$ on $\mathit{M}$ with ${\rho }_{\mathrm{bg}}$   5:    Compute initial adjoint source ${\mathit{s}}_{\mathrm{adj}}\leftarrow \mathrm{L}2\mathrm{Misfit}({\mathit{u}}_{\mathrm{init}},{\mathit{u}}_{\mathrm{obs}})$   6:    Compute initial gradient $\nabla {\mathcal{J}}_{\mathrm{elm}}\leftarrow \mathrm{SalvusAdjoint}({\rho }_{\mathrm{bg}},{\mathit{s}}_{\mathrm{adj}})$   7:    Compute wavefield energy ${\mathit{u}}_{\mathrm{norm}}\leftarrow {\sum }_t{\sum }_i{\parallel {\mathit{u}}_{\mathrm{init}}\parallel }^2$   8:    Clip and scale the initial gradient ${(\nabla {\mathcal{J}}_{\mathrm{elm}})}_{\mathrm{norm}}\leftarrow \mathrm{ScaleAndNormalize}(\nabla {\mathcal{J}}_{\mathrm{elm}},{\mathit{u}}_{norm})$   9:    return ${(\nabla {\mathcal{J}}_{\mathrm{elm}})}_{\mathrm{norm}}$ 10: End Function 11: $\nabla {\mathcal{J}}_{\mathrm{init}}\leftarrow \mathrm{InitialGradient}(\mathit{M},{\mathit{u}}_{\mathrm{obs}})$ 12: Map initial gradient to graph vertices ${(\nabla {\mathcal{J}}_{\mathrm{init}})}_{\mathrm{ver}}\leftarrow \mathrm{MeshToGraph}(\nabla {\mathcal{J}}_{\mathrm{init}})$ 13: Generate the graph connectivities $\mathit{E}$ 14: Load the pre-trained GCN network ${\theta }_t\leftarrow$ LoadWeights 15: Initialize ${\rho }_{ver}$ as torch parameter 16: GCN prediction at mesh vertices ${\rho }_{\mathrm{ver}}\leftarrow \mathcal{G}({(\nabla {\mathcal{J}}_{\mathrm{init}})}_{\mathrm{ver}},\mathit{E}).\mathrm{detach}\left(\right)*{\rho }_{\mathrm{bg}}$ 17: $\mathrm{optimizer}\leftarrow \mathrm{Adam}(\left[{\rho }_{\mathrm{ver}}\right],\mathrm{lr})$ 18: for$\mathrm{it}=1$toepochsdo 19:     $\mathrm{optimizer}.\mathrm{zero\_grad}\left(\right)$ 20:     Map density values to Elements ${\rho }_{\mathrm{elm}}\leftarrow \mathrm{GraphToMesh}\left({\rho }_{\mathrm{ver}}\right)$ 21:     Update mesh $\mathit{M}$ with current ${\rho }_{\mathrm{elm}}$ model 22:     ${\mathit{u}}_{\mathrm{syn}}\leftarrow \mathrm{SalvusForward}\left(\mathit{M}\right)$ 23:     $\mathcal{L}\leftarrow \mathrm{L}2\mathrm{Misfit}({\mathit{u}}_{\mathrm{syn}},{\mathit{u}}_{\mathrm{obs}})$ 24:     ${\mathit{s}}_{\mathrm{adj}}\leftarrow \mathrm{GenerateAdjointSource}\left(\mathcal{L}\right)$ 25:     $\nabla {\mathcal{J}}_{\mathrm{elm}}\leftarrow \mathrm{SalvusAdjoint}({\rho }_{\mathrm{elm}},{\mathit{s}}_{\mathrm{adj}})$ 26:     $\nabla {\mathcal{J}}_{\mathrm{elm}}\leftarrow \mathrm{Clip}(\nabla {\mathcal{J}}_{\mathrm{elm}},\text{clip\_fraction})$ 27:     Map gradient values to graph vertices $\nabla {\mathcal{J}}_{\mathrm{ver}}\leftarrow \mathrm{MeshToGraph}(\nabla {\mathcal{J}}_{\mathrm{elm}})$ 28:     ${\rho }_{\mathrm{ver}}.\mathrm{grad}=\nabla {\mathcal{J}}_{\mathrm{ver}}$ 29:     optimizer.step() 30: end for