Multiscale graph neural networks with adaptive mesh refinement for accelerating mesh-based simulations

In this work, we develop the multiscale mesh-based GNN framework with BSAMR for PF fracture problems. PF fracture problems are computationally demanding due to a) near singular operator due to near-zero modulus inside the crack field and b) the requirement of high mesh resolution near the crack tip. In the PF formulation of fracture, the sharp crack is regularized using a smooth scalar field $\phi(\bm{x})$ ranging from 0 to 1. An energy functional $\Pi$ is constructed using $\phi(\bm{x})$ which is then minimized to obtain a set of coupled partial differential equations, a vector equation corresponding to elastic equilibrium, and a scalar equation for crack equilibrium. The reader is referred to the extensive literature on PF fracture for details on the formulation [67]. In this work, we use the second-order energy functional as below [4].

Here, $\bm{u}$ is the displacement field, $\bm{\varepsilon}$ is the strain, $W$ is the strain energy density, $G_{c}$ is the fracture energy constant, and $d$ is the regularization length scale of the crack. Here, $\phi$ takes the value $0$ inside the crack, and the value $1$ in the bulk material. In this work, we use the PF open-source high-fidelity model with BSAMR capabilities from [4] to simulate various crack propagation problems. Following [43], we chose a linear elastic isotropic brittle material $(0.5m\times 0.5m)$, with Young’s modulus $E=210GPa$, Poisson’s ratio $\nu=0.3$, fracture energy density $G_{c}=2.7N/m^{2}$, and $d=0.0125m$. We fixed the bottom edge and applied displacement to the top edge in $y$-direction for tensile load and $x$-direction for shear load. We varied initial crack length, $C_{L}$, edge position, $C_{P}$, and crack angle, $C_{\theta}$ to gather the datasets of unique initial conditions.

As shown in Figure 2, we formulated the graph representation using the instantaneous refined mesh, $\mathcal{M}^{ref}:\langle{\mathcal{\mathbf{U}}},\mathbf{E}\rangle$, where $\mathcal{\mathbf{U}}$ includes the mesh vertices, and $\mathbf{E}$ the resulting mesh edges. The node features, $\xi_{s}$, include the positions of the mesh vertices, $\hat{P}_{s}$, their x- and y-displacement values, $\hat{D}_{s}$, their crack field values, ${\phi}_{s}$, and the applied displacement loading along the x- and y-directions, $\{{u}_{0_{s}},{v}_{0_{s}}\}$.

The applied displacement loading node feature, $\{{u}_{0_{s}},{v}_{0_{s}}\}$, allows the framework to distinguish between cases subjected to tensile loading versus cases under shear loading.

The edge features, $e_{sr}$, are defined using the binary connectivity value, $b_{sr}\in\{0,1\}$, where $s$ and $r$ indicate the “sender” and “receiver” nodes, respectively. This means that $b_{sr}=1$ for nodes that share an edge and for cases where $s=r$.

Lastly, we note that for each time-step in a simulation, we construct the initial graph from the instantaneous refined mesh containing all five levels of mesh resolution, $\mathcal{M}^{ref}$, shown in Figure 2(a). However, as shown in Figure 1, the framework removes each refinement level iteratively until obtaining the coarsest mesh at level 0 resolution. We define the corresponding graphs as (i) $\mathcal{M}^{ref}_{4}$ the initial instantaneous refined mesh with resolution levels 0-4 (shown in Figure 2(a)), (ii) $\mathcal{M}^{ref}_{3}$ the first downscaled mesh with resolution levels 0-3, (iii) $\mathcal{M}^{ref}_{2}$ the second downscaled mesh with resolution levels 0-2 (shown in Figure 2(b)), (iv) $\mathcal{M}^{ref}_{1}$ the third downscaled mesh with resolution levels 0-1, and (v) $\mathcal{M}^{ref}_{0}$ the coarsest downscaled mesh with resolution level 0 only (shown in Figure 2(c)).

The developed GNN framework integrates AMR to the multigrid approach where the initial refined mesh is coarsened by one mesh level at each step. The resulting coarsened mesh resolutions provide new connections between distant nodes to transfer information beyond the previous local neighborhoods. This approach reduces the required number of MP steps, thus, avoiding over-smoothing and significantly reducing computational costs while maintaining high accuracy.

The architecture of the developed multiscale GNN framework is depicted in Figure 1. The first step of the framework is to input the graph representation defined in Section 2.2 and shown in Figure 2(a), at a given time “$t$”, into an encoder network. We use an MLP model as the encoder network, denoted by $MLP_{\text{in}}$.

We input the generated latent-space embedding, $\xi_{s}^{{}^{\prime}}$, to a MP network, $GNN_{1}$, to learn relations within the local neighborhoods. We employ Graph Transformers as our MP models in this work [68]. Graph Transformers were first introduced for the tasks of sequence modeling and language translation [69]. Recently, this technique was extended for graph-based tasks by applying the multi-head and self-attention mechanisms to the vertex and edge features of neighboring vertices, thus, creating additional attention embeddings [68]. The Transformer MP models in the developed framework follow the same architecture, involving four attention heads and a hidden node dimension of 128.

The “GNN_D1” block depicted in Figure 1 includes the MP network, $GNN_{1}$, and the first Downscale operation, $Down_{1}$. Each Downscale operation removes the current finest level of mesh refinement. For instance, the first Downscale operation in Figure 1 removes refinement level 4 from graph $\mathcal{M}^{ref}_{4}$ (Figure 2(a)), resulting in graph $\mathcal{M}^{ref}_{3}$.

Here, $\xi_{s}^{4}$ is the resulting node embedding in refined mesh $\mathcal{M}^{ref}$ from the MP block GNN₁, $e_{sr}^{3}$ is the edge feature vector for the downscaled mesh $\mathcal{M}^{ref}_{3}$, and $\xi_{s}^{3}$ is the resulting downscaled node embedding from $\mathcal{M}^{ref}$ to $\mathcal{M}^{ref}_{3}$. Similarly, “GNN_D2” from Figure 1 includes a MP network, $GNN_{2}$, and the second Downscale operation, $Down_{2}$, which removes refinement level 3 from graph $\mathcal{M}^{ref}_{3}$, resulting in graph $\mathcal{M}^{ref}_{2}$ shown in Figure 2(b).

where $e_{sr}^{2}$ is the edge feature vector for the downscaled mesh $\mathcal{M}^{ref}_{2}$, and $\xi_{s}^{2}$ is the resulting downscaled node embedding from $\mathcal{M}^{ref}_{3}$ to $\mathcal{M}^{ref}_{2}$. We repeat this process using two additional MP blocks each followed by their coarsening operations to obtain the downscaled node embedding (from $\mathcal{M}^{ref}_{1}$ to $\mathcal{M}^{ref}_{0}$) for the coarsest mesh, $\xi_{s}^{0}$, shown in Figure 2(c). The framework then includes an additional MP network (denoted by “GNN_n” in Figure 1) which operates on the resulting coarsest mesh node embedding.

Here, $e_{sr}^{0}$ is the edge feature vector for the downscaled mesh $\mathcal{M}^{ref}_{0}$ shown in Figure 2(c), $Agg$ is an aggregation operation for the skip connectors (shown as dashed green lines in Figure 1), and $\xi_{s}^{0^{{}^{\prime}}}$ is the resulting node embedding for mesh $\mathcal{M}^{ref}_{0}$ from the aggregation of the MP block GNN_n and the skip connector.

Next, we perform a series of upscale steps to reconstruct the original refined mesh. For instance, “GNN_U1” includes the first Upscale operation, $Up_{1}$, to reconstruct mesh level 1 ($\mathcal{M}^{ref}_{1}$), followed by the skip connector from $\xi_{s}^{1}$ to the MP GNN, $GNN_{4}$.

In equation (8), $\xi_{s}^{1^{{}^{\prime}}}$ defines the resulting reconstructed node embedding from the upscale operation, $Up_{1}$, (i.e., from $\mathcal{M}^{ref}_{0}$ to $\mathcal{M}^{ref}_{1}$), and $\xi_{s}^{1^{{}^{\prime\prime}}}$ denotes the resulting node embedding for mesh $\mathcal{M}^{ref}_{1}$ from the aggregation of the MP block $GNN_{4}$ and the skip connector from $\xi_{s}^{1}$. Next, “${GNN}_{U2}$” applies the second Upscale operation, $Up_{2}$, to regenerate mesh level 2 ($\mathcal{M}^{ref}_{2}$ from Figure 2(b)), followed by the skip connector from $\xi_{s}^{2}$ to the MP GNN, ${GNN}_{3}$.

where $\xi_{s}^{2^{{}^{\prime}}}$ defines the resulting reconstructed node embedding from the upscale operation, $Up_{2}$, (i.e., from $\mathcal{M}^{ref}_{1}$ to $\mathcal{M}^{ref}_{2}$), and $\xi_{s}^{2^{{}^{\prime\prime}}}$ denotes the resulting node embedding for mesh $\mathcal{M}^{ref}_{2}$ from the aggregation of the MP block $GNN_{3}$ and the skip connector from $\xi_{s}^{2}$. We repeat this process using two additional upscale blocks (i.e., $GNN_{U3}$ and $GNN_{U4}$ in Figure 1) until we reconstruct the initial refined mesh $\mathcal{M}^{ref}$ with mesh levels 0-4 and obtain the resulting aggregated node embedding $\xi_{s}^{4^{{}^{\prime\prime}}}$. Lastly, we use the final generated reconstructed embedding, $\xi_{s}^{4^{\prime\prime}}$, as input to a Decoder MLP network, $MLP_{\text{out}}$, to transfer the embedding from the latent-space to the real-space and predict the displacements and scalar damage field at the next time-step, “$t+1$”.

While a higher number of downscale/upscale operations and MP GNNs may provide the framework with higher prediction accuracy, it comes with an increased computational cost. We evaluated the prediction accuracy and computational cost of three different architectures following the procedure described in the previous Section 2.3. For each architecture, we reduced the number of MP GNNs, upscale operations, and downscale operations. The first framework was a Four-Stage Refinement (FSR) GNN involving four downscale and four upscale operations as shown in Figure 1. The FSR framework architecture is described in detail in Section 2.3.

The second framework was a Two-Stage Refinement (TSR) GNN involving two downscale and two upscale operations. In the TSR framework, a single downscale/upscale operation accounts for two mesh resolution levels. For instance, given an initial instantaneous refined mesh $\mathcal{M}^{ref}_{4}$ (Figure 2(a)), the first downscale operation in TSR, $GNN_{D1}$, removes mesh resolution levels 4 and 3 to obtain the new mesh graph $\mathcal{M}^{ref}_{2}$ (i.e., Figure 2(b)), along with the resulting downscaled node embedding $\xi_{s}^{2}$. By removing four MP GNNs, and two downscale and upscale operations, we expect the TSR framework to be faster than FSR.

Lastly, the third framework was a Single-Stage GNN (SSR) involving a single downscale/upscale operation. In the SSR framework, we removed six MP GNNs and three downscale and upscale operations. For instance, SSR involves a single downscale operation, $GNN_{D1}$, which removes mesh resolution levels 1-4 directly from the initial node embedding, $\xi_{s}^{4}$, to obtain the coarsest mesh level shown in Figure 2(c). The SSR framework is significantly less computationally expensive compared to FSR and TSR as it does not require computing new graphs for intermediate mesh resolution levels (1-4), nor additional MP Transformer GNNs.

First, we trained the multiscale GNN framework for the left-edge notched system under tensile loading, shown in Figure 3(a). We obtained a dataset of 1100 PF fracture simulations involving single-edge notched systems under tensile loading from [43]. Each of these simulations consisted of a unique crack configuration where the initial length $C_{L}$, edge position $C_{P}$, and orientation of the crack $C_{\theta}$ were varied. Additional details for the simulation set-up can be found in [43].

To perform TL, we first trained the multiscale GNN framework on the large dataset of 1100 single-edge notched simulations gathered from [43]. From the trained GNN framework, we transferred the weights of the encoder MLP (i.e., $MLP_{\text{in}}$), and the first MP model (i.e., GNN_D1) to a new GNN framework with similar architecture. Using this approach, we performed a series of sequential TL update steps. For each TL update step, we transferred the resulting pretrained weights to the next case study. We considered the following 4 case studies.

• 1.

  Case 1: Left-edge cracks subjected to Mode I loading as shown in Figure 3(a).

• 2.

  Case 2: Center cracks subjected to Mode I loading as shown in Figure 3(b).

• 3.

  Case 3: Left-edge cracks subjected to Mode II loading as shown in Figure 3(c).

• 4.

  Case 4: Right-edge cracks subjected to Mode I loading as shown in Figure 3(d).

Unlike Case 1 where left-edge cracks are only allowed to propagate towards the right, in Case 2 center cracks can propagate in both the left and right directions, thus, increasing the framework’s understanding. For Cases 2-3, we leveraged the PF model from [4] to generate new datasets for center cracks and shear loading scenarios. For each of these cases, we generated a dataset of 30 simulations (i.e., 15 for training, and 15 for testing). For Case 4 we leveraged the symmetry of the left-edge crack case study (i.e., Case 1) to mirror 30 randomly chosen simulations from the training dataset used in Case 1. We emphasize the significant decrease in the required size of the training dataset from 1100 (Case 1) to 15 simulations using TL (i.e., approximately 70 times smaller). The implementation of TL allows the multiscale GNN framework to be extended to other crack problems with a fast simulation time and decreased computational cost.

Here, we compare the performance of FSR, TSR, and SSR architectures of the multiscale GNN framework. We obtained the predictions of the crack variable $\phi$, $x$-displacement, and $y$-displacement for the left-edge crack cases. Figures 4(a) - 4(c) depict the predicted $\phi$ values, $x$-displacement, and $y$-displacement for a randomly chosen simulation from the test dataset of left-edge crack cases. The simulation chosen shows a crack of positive orientation, positioned towards the top of the left edge of the domain. The crack can also be seen approaching the right edge of the domain for complete material failure. From Figure 4(a) all three architectures predicted $\phi$ with nearly identical results to the ground truth (i.e., left-most case). We obtained a similar result for x- and y-displacements shown in Figures 4(b) - 4(c), where qualitatively the FSR, TSR, and SSR architectures predicted displacements with high accuracy. This qualitative analysis demonstrates that the TSR and SSR were able to maintain their prediction accuracy despite the reduced downsampling and upsampling steps.

Next, we computed the errors corresponding to each simulation from the test dataset. For each test simulation, we first computed the average percent error across all mesh points in $\mathcal{M}^{ref}$ for each time step. We then computed the average of the resulting percent error across all time steps of the simulations. For instance, we computed the error for the crack field as

were $\phi_{s}^{pred}$ and $\phi_{s}^{true}$ denote the predicted and ground truth crack field value at each mesh point $s\in\mathcal{M}^{ref}$, respectively, $M$ is the total number of mesh points in $\mathcal{M}^{ref}$, $t$ is the first predicted time-step, and $T_{f}$ the final time-step. Using equation (11), we gathered the $\phi$ and displacement errors corresponding to each test simulation for the FSR, TSR, and SSR architectures. Figures 6(a) - 6(c) show the resulting crack field and displacement errors for the FSR architecture. Figures 6(d) - 6(f) depict the errors for the TSR architecture. Similarly, Figure 6(g) - 6(i) shows the resulting average errors for the SSR architecture.

Comparing $\phi$ predictions (i.e., left-most), we note that the FSR, TSR, and SSR architectures resulted in average percent errors below 0.3$\%$. The errors in $x$-displacement and $y$-displacement also remained below 0.35$\%$ for all architectures. These results confirm our previous observations from the qualitative analysis which indicated that TSR and SSR were capable of maintaining high prediction accuracy despite their reduced downsampling and upsampling steps.

Next, we computed the average errors across all testing simulations for each architecture 7(a). From Figure 7(a), we note that FSR architecture demonstrated the highest accuracy when predicting the crack field $\phi$. However, the SSR architecture showed lower errors in crack field predictions than the TSR architecture. The TSR architecture resulted in the lowest error for $x$-displacement prediction, while the FSR and SSR architectures showed similar errors at approximately 0.08$\%$. For $y$-displacement predictions, the FSR and TSR architectures showed similar low errors at approximately 0.08$\%$, while the SSR architecture showed the highest error close to 0.12$\%$. Finally, Figure 7(a) shows that reducing the number of refinement operations for TSR and SSR did not significantly increase prediction error. The highest prediction error of 0.12$\%$ for SSR on $y$-displacement predictions is still considerably low compared to previous work [43].

We compared the computational cost for FSR, TSR, and SSR architectures by calculating the time required to generate 30 simulations. In Section 3.1, we showed that the prediction errors did not significantly increase for the TSR and SSR architectures despite their lower number of downscaling and upscaling operations. Each operation increases computational costs because it requires storing the resulting mesh configurations and their node embeddings. Each operation also comes with added MP GNNs. Figure 7(b) depicts the total time in hours required for each framework to generate 30 randomly chosen simulations from the training dataset. As shown in Figure 7(b), the FSR GNN architecture required the longest simulation time of 6.13 hours. The TSR architecture was the second most computationally demanding, requiring 4.75 hours to simulate 30 cases. As expected, the fastest architecture with the lowest computational costs was the SSR, which required only 3.91 hours. In contrast, we note from [43] that the high-fidelity PF model required approximately 43.5 hours to generate 30 cases. Due to its lower simulation time and high prediction accuracy, we chose the SSR architecture for subsequent TL steps.

In Sections 3.1 - 3.2, we determined that the SSR architecture provided high prediction accuracy at significantly lower computational costs. We implemented TL to the multiscale GNN framework with SSR architecture for cases with center cracks as shown in Figure 3(b). We tested the extended SSR-based framework to predict the crack and displacement fields for a randomly chosen case from the testing data. For this analysis, we first obtained the initial time-step prediction ($t_{0}$) as shown in Figure 9(a) for the ground truth versus predicted crack field, Figure 9(c) for the ground truth versus predicted $x$-displacement, and Figure 9(e) for the ground truth versus predicted $y$-displacement. Next, we propagated the simulation in time and obtained the corresponding predictions of crack and displacement fields for a time-step approaching complete material failure ($t_{f}$).

Figures 9(b), 9(d), and 9(f) depict the ground truth versus predicted crack field, $x$-displacement, and $y$-displacement, respectively, at $t_{f}$. These results depict similar accuracies compared to the left-edge crack cases. During both time steps (i.e., $t_{0}$ and $t_{f}$) the SSR-based framework’s predictions are virtually indistinguishable from the high-fidelity PF model. These qualitative results demonstrate that the SSR-based framework was capable of predicting crack propagation for center crack cases through TL.

We also performed a quantitative analysis to verify these qualitative observations by generating the errors of each test simulation. Figure 10 depicts the computed average percent errors for center crack test cases obtained using equation (11). For crack field errors, the extended SSR framework maintained high prediction accuracy (less than 0.125$\%$ error) across all test cases. Similarly, the $x$-displacement errors remained below 0.25$\%$ error, and y-displacement errors below 0.20$\%$. These results emphasize the high prediction accuracy of the SSR framework by implementing TL using significantly smaller training datasets. Additionally, extending the SSR framework from left-edge cracks to center cracks increases the framework’s capabilities in predicting cases where cracks propagate in both the left and right directions.

In this step, we apply TL update to the extended SSR-based framework obtained in Section 3.3. As mentioned in Section 2.2, the node features $\{u_{0},v_{0}\}$ allow the framework to distinguish between tensile loading (e.g., $u_{0}=\Delta u$ and $v_{0}=0$) and shear loading (e.g., $u_{0}=0$ and $v_{0}=\Delta v$) cases. Following a similar approach as for center cracks, we evaluated the resulting framework for shear load cases to predict the crack field, and displacement fields at an initial time-step, $t_{0}$, and at a later time-step approaching material failure, $t_{f}$.

Figures 12(a), 12(c), and 12(e) compare GNN predictions and PF predictions for the crack field, $x$-displacement, and $y$-displacement fields, respectively for initial time steps. Figures 12(b), 12(d), and 12(f) compare GNN and PF predictions for the crack field, $x$-displacement, and $y$-displacement fields, respectively for the final time step. These figures show that the new extended SSR framework was able to capture cases involving shear loads and generate accurate predictions at both time steps with virtually identical results compared to the ground truth. For quantitative analysis, we computed the average errors for each test simulation following the approach described in Section 3.1. Figure 13 shows the resulting average percent errors for all test cases under shear loads. As shown, the SSR also predicted shear cases with high accuracy. Average percent errors for the crack field, $x$-displacement, and $y$-displacement resulted below 0.25$\%$, 1.20$\%$, and 0.25$\%$, respectively.

Lastly, we extended the framework for cases involving right-edge cracks. For this TL update step, we mirrored the dataset involving left-edge cracks subjected to tension as described in Section 2.5. Following the same approach as in previous sections, we first tested the resulting trained SSR framework to simulate a randomly chosen simulation from the test dataset. Figures 15(a) - 15(b) depict the prediction history of the crack field parameter from $t_{0}$ to $t_{f}$ for the chosen test simulation. The prediction history depicts qualitatively identical results compared to the PF predictions. This high prediction accuracy is verified in Figure 16(a), where the average error across all simulations remains below 0.10$\%$. For the x-displacement and y-displacement predictions, we show the resulting SSR simulation histories in Figures 15(c) - 15(f), respectively. The SSR framework also shows high prediction accuracy compared to the PF simulations. The obtained average percent errors for x-displacement are shown in Figures 16(b). The errors across all testing samples remained below 0.25$\%$. Also, for y-displacement errors shown in Figure 16(c), it may be noted that the average percent errors remained below 0.30$\%$. Ultimately, these results show that by implementing a series of sequential TL update steps, the SSR-based framework was successful in simulating multiple problem configurations with high accuracy.