HI-MeshGraphNets: Efficient and Accurate Mesh-based Physics Learning with Hierarchical Multi-scale Graph Neural Networks

# Efficient and Accurate Mesh-based Physics Learning with Hierarchical Multi-scale Graph Neural Networks
Source: [https://arxiv.org/html/2608.13827](https://arxiv.org/html/2608.13827)
## HI\-MeshGraphNets: Efficient and Accurate Mesh\-based Physics Learning with Hierarchical Multi\-scale Graph Neural NetworksJournal:Elsevier

SiHun LeeEmail:[s\.hun\.lee@samsung\.com](mailto:[email protected])Affiliation:Mobile Experience, Digital Twin AX Lab\., Samsung Electronics Co\., 129, Samsung\-ro, Yeongtong\-gu, Suwon\-si, 16677, Gyeonggi\-do, Republic of KoreaDong\-Hyuk ParkAffiliation:Department of Aerospace system Engineering, Sejong University, 209, Neungdong\-ro, Gwangjin\-gu, Seoul, 05006, Republic of KoreaTaesoo BangAffiliation:Department of Aerospace system Engineering, Sejong University, 209, Neungdong\-ro, Gwangjin\-gu, Seoul, 05006, Republic of KoreaSeung\-Hoon KangEmail:[shkang@sejong\.ac\.kr](mailto:[email protected])Corresponding author:Corresponding authorAffiliation:Department of Aerospace system Engineering, Sejong University, 209, Neungdong\-ro, Gwangjin\-gu, Seoul, 05006, Republic of Korea

###### Abstract

Machine\-learned physical surrogate models have become promising alternatives to mesh\-based numerical solvers\. Among them, graph neural networks \(GNNs\) have been successful in representing simulation meshes as graphs and learning nodal state evolution via message passing\. However, flat message passing can become inefficient for large, high\-fidelity meshes or significant variations within a step\. Since information propagates only one hop per step, long\-range interactions require many hops, increasing computational burden, memory usage, and the risk of over\-smoothing\.

To mitigate the long\-range interaction problem, we propose Hierarchical Interpolating MeshGraphNets \(HI\-MGN\), a flexible multiscale extension of GNNs designed to improve long\-range communication on large unstructured meshes\. HI\-MGN replaces the flat processor with a hierarchical multiscale processor\. The hierarchical processor works by coarsening the graphs with farthest\-point sampling \(FPS\) and Voronoi partitioning\. Such coarsening preserves the underlying mesh topology without introducing artificial spatial connections or nodes\. These coarse graphs enable information to propagate over larger geometric distances with fewer message passing steps\. After message passing on the coarse graphs, the information is projected back onto the original mesh resolution via a learned graph interpolation network\.

Through numerical comparisons on three structural and fluid problems against existing GNN methods, such as MeshGraphNets and the Bi\-Stride Multi\-Scale GNN \(BSMS\-GNN\), we show that our method achieves improved accuracy with a comparable number of message passing blocks while reducing training time and peak memory usage\. These results demonstrate that topology\-aware hierarchical multiscale processing and learned coarse\-to\-fine interpolation provide an effective and practical route toward scalable mesh\-based physics surrogate modeling\.

###### Keywords:

Graph Neural Network , Mesh\-based Simulation , Multi\-scale , Physics AI , Scientific Machine Learning

††highlights:HI\-MGN improves long\-range communication through hierarchical message passing\. FPS–Voronoi coarsening builds coarse graphs from original mesh adjacency\. Learned interpolation reconstructs fine\-scale features after coarse processing\. HI\-MGN improves accuracy while reducing memory and training cost\.## 1Introduction

As the fields of artificial intelligence \(AI\) and machine learning \(ML\) have progressed rapidly, efforts to build ML\-based physics surrogate models have also increased\. Also known as scientific machine learning \(SciML\) or AI for computer aided engineering \(AI\-CAE\), earlier methods mainly focused on applying classical regressors such as Gaussian process regression and radial basis function regression to linearly reduced subspaces such as proper orthogonal decomposition\[[37](https://arxiv.org/html/2608.13827#bib.bib1),[26](https://arxiv.org/html/2608.13827#bib.bib2)\]\. Following these regressors, simple feed\-forward neural networks\[[8](https://arxiv.org/html/2608.13827#bib.bib3),[20](https://arxiv.org/html/2608.13827#bib.bib4),[13](https://arxiv.org/html/2608.13827#bib.bib5)\], long short\-term memory networks\[[25](https://arxiv.org/html/2608.13827#bib.bib18),[35](https://arxiv.org/html/2608.13827#bib.bib6),[6](https://arxiv.org/html/2608.13827#bib.bib7)\], variational autoencoders \(VAEs\)\[[16](https://arxiv.org/html/2608.13827#bib.bib9)\], and generative adversarial networks \(GANs\)\[[15](https://arxiv.org/html/2608.13827#bib.bib8),[10](https://arxiv.org/html/2608.13827#bib.bib14)\]were considered upon the linear subspace\. Then, as ML methods improved, fully data\-driven ML methods were used as physical surrogate models\. Models such as autoencoders\[[11](https://arxiv.org/html/2608.13827#bib.bib21),[39](https://arxiv.org/html/2608.13827#bib.bib20),[14](https://arxiv.org/html/2608.13827#bib.bib15)\], and VAEs\[[33](https://arxiv.org/html/2608.13827#bib.bib16),[12](https://arxiv.org/html/2608.13827#bib.bib17),[34](https://arxiv.org/html/2608.13827#bib.bib19),[17](https://arxiv.org/html/2608.13827#bib.bib13),[18](https://arxiv.org/html/2608.13827#bib.bib12)\]were investigated\. While these methods have been successful on structured grids or fixed discretizations, they can quickly lose effectiveness on unstructured meshes and geometry\-deforming problems\.

Recent progress in geometry\-aware learning has substantially changed this landscape by introducing mesh\-agnostic methods that enable the use of various meshes for training and inference\. Early mesh\-agnostic attempts emerged mainly from point\-cloud methods and soon extended to graph neural networks \(GNNs\) and discretization\-independent operator\-learning methods\. Point\-set networks such as PointNet, PointNet\+\+, and related point\-based convolutions were among the early methods explored for mesh\-agnostic surrogate models\[[30](https://arxiv.org/html/2608.13827#bib.bib22),[31](https://arxiv.org/html/2608.13827#bib.bib23),[38](https://arxiv.org/html/2608.13827#bib.bib24)\]\. However, point\-based convolutional networks require a standardized set of points distributed in a point\-cloud\. This can result in decreased performance in representing local geometry, especially when the original degrees of freedom are very large\. The PointNet architecture may also be prone to errors when the original dataset is complex and has highly diverse point densities\.

In parallel, neural operator learning has provided another important route toward mesh\-agnostic surrogate modeling\. Rather than learning a regression map tied to a particular mesh or discretization, operator\-learning methods aim to approximate the solution operator that maps input functions to output functions\. DeepONet was one of the earliest and most influential formulations of this idea along with subsequent geometry\-aware extensions, such as Geom\-DeepONet and Fusion\-DeepONet\[[23](https://arxiv.org/html/2608.13827#bib.bib25),[7](https://arxiv.org/html/2608.13827#bib.bib26),[28](https://arxiv.org/html/2608.13827#bib.bib27)\]\. Additionally, Fourier\- and geometry\-based neural operators further extended this operator\-learning direction\. Geo\-FNO generalized the Fourier Neural Operator to irregular domains by learning a deformation from physical geometry to a latent computational domain\[[21](https://arxiv.org/html/2608.13827#bib.bib28)\]\. Geometry\-informed neural operators introduced richer geometric encodings, including point\-cloud, signed\-distance, and graph\-based representations, to improve scalability for large three\-dimensional PDE problems\[[22](https://arxiv.org/html/2608.13827#bib.bib29)\]\. While these operator\-learning methods have been successful to some degree, they still face challenges when applied to industrial CAE settings, particularly in handling geometric, boundary condition variations, preserving fine\-scale geometric features, and scaling to large unstructured meshes with localized gradients or shock\-induced discontinuities\. It is also known that many of operator learning methods in practice, can suffer from out\-of\-bounds extrapolation\. More recently, transformer\-based PDE solvers such as Transolver and Transolver\+\+ have introduced physics\-attention mechanisms that communicate through learned physical slices instead of explicit connectivity\[[36](https://arxiv.org/html/2608.13827#bib.bib30),[24](https://arxiv.org/html/2608.13827#bib.bib32)\]\. While these methods have proven to be accurate, they are known to be sensitive to hyperparameters and may require substantial computational resources for training compared to existing methods\[[1](https://arxiv.org/html/2608.13827#bib.bib11)\]\.

As meshes in numerical simulations are naturally graphs, GNN\-based mesh\-agnostic methods that treat simulation meshes as graphs have also been investigated\. One of the earliest success with GNN\-based physics surrogate models was MeshGraphNets \(MGN\)\[[29](https://arxiv.org/html/2608.13827#bib.bib33)\]\. Following MGN, many GNN frameworks have used message passing processors to propagate information across the mesh, learning the evolution of physical states\. While GNN\-based frameworks have been successful, two main bottlenecks exist: excessive memory requirements as mesh size increases and limited long\-range interaction on large graphs\. The former can be mitigated by decreasing batch size and using better hardware, while the latter remains a more fundamental issue\. Counterintuitively, GNN\-based methods can show deteriorated performance as the number of degrees of freedom increases, contrary to the behavior often expected from classical numerical simulations\. This is caused by the nature of message passing, where each layer propagates information only between directly connected nodes\. In traditional message passing, a node’s receptive field grows by a single hop per layer\. Communicating between two physically distant regions therefore requires a number of message passing steps proportional to their geodesic distance on the graph\. Therefore, long\-range interactions in a mesh require deeper processor stacks, resulting in over\-smoothing and information attenuation\.

To alleviate limited long\-range interaction in GNNs, multiscale graph methods such as MultiScale MGN and Bi\-Stride Multi\-Scale GNN \(BSMS\-GNN\) have been proposed\. MultiScale MGN introduced coarse\-resolution message passing to improve communication on high\-resolution systems, while BSMS\-GNN proposed bi\-stride coarsening based on breadth\-first search \(BFS\), avoiding manually generated coarse meshes and reducing erroneous spatial\-proximity connections\[[5](https://arxiv.org/html/2608.13827#bib.bib34),[4](https://arxiv.org/html/2608.13827#bib.bib35)\]\. However, these methods still have limitations\. MultiScale MGN requires a pre\-generated set of coarse meshes in addition to the fine mesh for training and inference by design\. BSMS\-GNN uses repeated bi\-stride coarsening, which may require excessive hierarchical levels to achieve sufficiently coarse representations on very large meshes\. More recently, methods such as X\-MeshGraphNet introduced scalable METIS partitioning and halo\-region strategies for large engineering graphs\. It enables the handling of larger domains by decomposing the original graph into smaller subgraphs while preserving inter\-partition communication\[[27](https://arxiv.org/html/2608.13827#bib.bib36)\]\. Such an approach may be beneficial for memory efficiency or distributed scalability; however, partitioning alone does not resolve the long\-range interaction within each high\-resolution subgraph\. Moreover, long\-range information may be compressed through narrow graph bottlenecks, resulting in the over\-squashing problem\.

Recent global\-processing approaches such as MeshGraphNet\-Transformer \(MGN\-T\) address long\-range interactions using attention\-based global processors without hierarchical coarsening\[[9](https://arxiv.org/html/2608.13827#bib.bib10)\]\. Also, Physics\-Informed Ollivier–Ricci Flow \(PIORF\) introduced a graph\-rewiring approach for improving long\-range interaction in mesh\-based GNNs\[[40](https://arxiv.org/html/2608.13827#bib.bib40)\]\. PIORF uses Ollivier–Ricci curvature to identify graph bottlenecks and introduces physics\-informed artificial long\-range connections toward regions with large velocity gradients, thereby mitigating the over\-squashing problem\. Unlike hierarchical coarsening approaches, however, PIORF modifies the computational graph connectivity by introducing artificial additional edges\. While such rewiring can improve long\-range information propagation, the original mesh\-edge topology may not be strictly preserved\. This may require additional connectivity constraints in multi\-body or contact problems, where unintended connections between physically disconnected components can be undesirable\.

In this work, we proposeHierarchical Interpolating MeshGraphNets \(HI\-MGN\), a MeshGraphNets\-style surrogate model designed for efficient and accurate full\-field prediction with significantly enhanced long\-range interactions\. HI\-MGN preserves the standard encode–process–decode structure of MeshGraphNets, but replaces the flat processor with a configurable multi\-scale hierarchical processor\. The model first applies local message passing on the fine mesh, constructing latent state considering its neighbors\. Then the latent information are pooled onto a coarser graph\. On the coarser graph, the model performs another set of message passing operations to enlarge the receptive field and then reconstructs the fine latent representation from the coarse mesh\. Such hierarchical design enables long\-range interaction via the coarse graph and considerably alleviates the long\-range communication bottleneck with only a few message passing steps\.

The main contribution of this work is a self\-contained hierarchical graph construction and interpolation framework for mesh\-based GNN surrogates\. Unlike multiscale approaches that rely on prescribed multi\-resolution meshes or fixed recursive coarsening rules, HI\-MGN constructs its hierarchy directly from a single input mesh using farthest\-point sampling \(FPS\)–Voronoi clustering\. Coarse connectivity is induced from the original mesh adjacency, preventing geometric proximity alone from introducing connections between nearby but disconnected components\. The resulting hierarchy permits strongly reduced coarse representations with user\-specified node counts\. To recover fine\-scale information after coarse\-level processing, we introduce a learned geometry\-aware interpolation operator that combines neighboring coarse latent states, fine\-resolution skip features, and relative geometric positions\. Through these components, HI\-MGN aims to provide an efficient and mesh\-native surrogate modeling framework for high\-resolution problems involving complex geometries and long\-range physical interactions\. Throughout the manuscript, we evaluate HI\-MGN against the original MGN and our most closely related method, BSMS\-GNN\. Across three numerical examples, we show that HI\-MGN achieves improved accuracy while substantially reducing training time and peak memory usage\.

## 2Methodology

### 2\.1Model definition

Throughout the model, follow the basic structure of MeshGraphNets\[[29](https://arxiv.org/html/2608.13827#bib.bib33)\]\. We model atttemporal discretized numerical simulation as a sequence of graphs\{G0,G1,…,Gt\}\\\{G^\{0\},G^\{1\},\\dots,G^\{t\}\\\}whereGt=\(Vt,Et\)G^\{t\}=\(V^\{t\},E^\{t\}\)denotes the graph at time steptt\. Each node in vertexVVstores nodal information such as deformation, physical quantities of interest, and categorical tags including boundary conditions, part numbers, and material types\. The edge set consists of two types of bidirectional edges: mesh edgesEME^\{M\}and world edgesEW,tE^\{W,t\}, such thatEt=EM∪EW,tE^\{t\}=E^\{M\}\\cup E^\{W,t\}\.

Mesh edgesEME^\{M\}follow the connectivity of the original numerical mesh\. Each mesh edge stores a relative geometric vector and its Euclidean norm,\{d​x,d​y,d​z,∥\(d​x,d​y,d​z\)∥\}\\\{dx,dy,dz,\\lVert\(dx,dy,dz\)\\rVert\\\}between the connected nodes\. The geometric features are evaluated in both the reference configuration and the configuration at timett, yielding an eight\-dimensional features\. World edgesEW,tE^\{W,t\}are constructed dynamically at each time step, connecting nodes whose current separation in the Lagrangian description falls below a certain threshold\. These edges are intended to represent non\-mesh interactions such as contact in structural dynamics and are constructed only at the original mesh resolution\.

HI\-MGN learns the temporal evolution of a system: givenGtG^\{t\}, it predicts the increment of the nodal physical state fromtttot\+1t\+1\. Denoting the physical state of nodeiiat timettby𝐮it\\mathbf\{u\}\_\{i\}^\{t\}, the network predictsΔ​𝐮it=𝐮it\+1−𝐮it\\Delta\\mathbf\{u\}\_\{i\}^\{t\}=\\mathbf\{u\}\_\{i\}^\{t\+1\}\-\\mathbf\{u\}\_\{i\}^\{t\}\. In the static case, the network instead predicts the nodal state difference𝐮i−𝐮i0\\mathbf\{u\}\_\{i\}\-\\mathbf\{u\}\_\{i\}^\{0\}from the reference state givenG0G^\{0\}\.

### 2\.2Architecture

HI\-MGN uses the standard ’Encoder–Processor–Decoder’ structure adopted by previously developed GNN surrogate models\[[2](https://arxiv.org/html/2608.13827#bib.bib38),[32](https://arxiv.org/html/2608.13827#bib.bib37),[29](https://arxiv.org/html/2608.13827#bib.bib33)\]\. Throughout the model, multi\-layer perceptrons \(MLPs\) with SiLU activations and LayerNorm are used\. However, HI\-MGN replaces the conventional flat processor with a hierarchical multiscale processor\. The encoder and decoder follow the same formulation as those of MeshGraphNets and are summarized below\[[29](https://arxiv.org/html/2608.13827#bib.bib33)\]\.

A nodal encoder𝔼V\\mathbb\{E\}\_\{V\}maps each nodal feature, consisting of the normalized physical state concatenated with one\-hot encoded categorical tags, to a nodal latent vectorziNz^\{N\}\_\{i\},

ziN=𝔼V​\(vi\),z^\{N\}\_\{i\}=\\mathbb\{E\}\_\{V\}\(v\_\{i\}\),\(1\)whereviv\_\{i\}denotes nodal state ofii\-th node\. Separate mesh and world edge encoders \(𝔼M\\mathbb\{E\}\_\{M\}and𝔼W\\mathbb\{E\}\_\{W\}\) embed the edge features of the corresponding two edge types\. The corresponding latent vectors become

zi​jM=𝔼M​\(ei​jM\),zi​jW=𝔼W​\(ei​jW\),z^\{M\}\_\{ij\}=\\mathbb\{E\}\_\{M\}\\\!\\big\(e^\{M\}\_\{ij\}\\big\),\\qquad z^\{W\}\_\{ij\}=\\mathbb\{E\}\_\{W\}\\\!\\big\(e^\{W\}\_\{ij\}\\big\),\(2\)whereei​jMe^\{M\}\_\{ij\}andei​jWe^\{W\}\_\{ij\}denotes mesh and world edges connecting nodeiiandjj\. The geometric information associated with both the reference and current configurations is provided to the network\.

A message passing \(MP\) block first updates the edge latent states based on the latent states of their corresponding sender and receiver nodes, and then updates the node latent states using the aggregated incoming edge messages, as follows:

zi​jM\\displaystyle z^\{M\}\_\{ij\}←zi​jM\+ϕM​\(\[ziN,zjN,zi​jM\]\),\\displaystyle\\leftarrow z^\{M\}\_\{ij\}\+\\phi\_\{M\}\\\!\\big\(\[\\,z^\{N\}\_\{i\},\\,z^\{N\}\_\{j\},\\,z^\{M\}\_\{ij\}\\,\]\\big\),\(3\)ziN\\displaystyle z^\{N\}\_\{i\}←ziN\+ϕV​\(\[ziN,∑jzi​jM,∑jzi​jW\]\),\\displaystyle\\leftarrow z^\{N\}\_\{i\}\+\\phi\_\{V\}\\\!\\Big\(\\big\[\\,z^\{N\}\_\{i\},\\;\\textstyle\\sum\_\{j\}z^\{M\}\_\{ij\},\\;\\textstyle\\sum\_\{j\}z^\{W\}\_\{ij\}\\,\\big\]\\Big\),\(4\)whereϕM\\phi\_\{M\}andϕV\\phi\_\{V\}are independently parameterized MLPs\. For the world edges, if present, are updated byϕW\\phi\_\{W\}and enter the node update through the aggregated messages\. The decoder𝔻\\mathbb\{D\}maps the refined message\-passed nodal latent of each node to a normalized per\-node physical state update,

y^i=𝔻⁡\(zi\)\.\\hat\{y\}\_\{i\}=\\mathbb\{D\}\(z\_\{i\}\)\.\(5\)The predicted update is then denormalized and accumulated onto the current physical state to obtain the next state,

uit\+1=uit\+𝒩−1​\(y^i\),u^\{t\+1\}\_\{i\}=u^\{t\}\_\{i\}\+\\mathcal\{N\}^\{\-1\}\(\\hat\{y\}\_\{i\}\),\(6\)where𝒩\\mathcal\{N\}denotes normalizing operator\.

At inference, this update is applied autoregressively\. The mesh\-edge topology remains fixed, while geometry\-dependent mesh\-edge features and the world\-edge set are recomputed from the predicted configuration at each time step\. For static problems, the decoder output is instead interpreted as the physical\-state difference relative to the reference configuration or the initial condition\.

### 2\.3Hierarchical multi\-scale processor

Rather than stacking processor blocks at a single, original resolution, HI\-MGN arranges message passing in a hierarchical manner overHHcoarse graphs as shown in Fig\.[1](https://arxiv.org/html/2608.13827#S2.F1)\. Hierarchical multi\-scale message passing makes the modeling of long\-range graph interactions tractable\. Considering that two nodes separated byddedges requireO⁡\(d\)O\(d\)blocks to communicate in a conventional flat processor, long\-range interaction may require prohibitively deep network that will frequently result in over\-smoothing on large or high\-fidelity meshes\[[19](https://arxiv.org/html/2608.13827#bib.bib39)\]\. A multi\-scale processor shortens these communication paths by performing message passing on coarser graphs\. A few blocks at the coarsest level can therefore exchange information across a much larger portion of domain, enabling information to propagate farther with the same number of message passing steps\.

The hierarchical multi\-scale processor consists of2​H\+12H\{\+\}1message passing stages, with\{L0,…,L2​H\}\\\{L\_\{0\},\\dots,L\_\{2H\}\\\}denoting the number of message passing layers at each stage\. The processor first performsL0L\_\{0\}message passes at the original level after which the latent states are pooled to the coarser graph\. At the next coarse level, the latent vectors are refined again forL1L\_\{1\}passes and then pooled onto the next coarse graph\. After the latent vectors are refined at the coarsest level, they are upsampled by the learned, geometry\-aware interpolation GNN\. The upsampled latent vectors are then merged with the corresponding skip states and refined once more to generate the final latent state at the original resolution\.

![Refer to caption](https://arxiv.org/html/2608.13827v1/architecture.png)Figure 1:Architecture of HI\-MGN with hierarchical, multi\-scale message passing \(MP\)\.A farthest\-point sampling \(FPS\) with Voronoi clustering \(hereafter, FPS\-Voronoi\) strategy is considered for the graph coarsening\. The FPS\-Voronoi coarsening requires no externally supplied meshes\. First,kkseeds are chosen by greedy FPS using the reference coordinates\. A multi\-source breadth\-first search over the original mesh adjacency is used to assign nodes to graph\-based Voronoi clusters\.

The coarse nodes consist of the FPS\-selected seeds\. The coarse node latent features are generated by mean\-pooling the latent states of the nodes within each Voronoi cluster, while the coordinates are inherited from the original FPS seeds\. After theL0L\_\{0\}fine\-level blocks, each node’s latent state contains information from up to itsL0L\_\{0\}\-hop neighbors\. For the edges, two coarse nodes are connected if at least one fine\-mesh edge crosses between their corresponding Voronoi clusters\. This boundary\-induced edge construction and coarse node selection process preserves mesh topology and, unlike proximity\-based coarsening, introduces no artificial connections across nearby but unconnected components\. ForH\>1H\>1the procedure is applied recursively, yielding multiple hierarchical levels\. Algorithm[1](https://arxiv.org/html/2608.13827#alg1)shows the FPS\-Voronoi coarsening method used in HI\-MGN\.

Algorithm 1FPS–Voronoi coarsening strategy of HI\-MGN\.1:Fine mesh graph

\(V,E\)\(V,E\), node positions

\{pi\}i∈V\\\{p\_\{i\}\\\}\_\{i\\in V\}, number of coarse nodes

kk
2:FPS sampling

3:

s1=v0s\_\{1\}=v\_\{0\},

S=\{s1\}S=\\\{s\_\{1\}\\\}⊳\\trianglerightSS: set of FPS seeds

4:

di=∥pi−ps1∥,∀i∈Vd\_\{i\}=\\lVert p\_\{i\}\-p\_\{s\_\{1\}\}\\rVert,~\\forall i\\in V
5:for

j=2,…,kj=2,\\dots,kdo

6:

sj=arg⁡maxi∈V∖S⁡dis\_\{j\}=\\arg\\max\_\{i\\in V\\setminus S\}d\_\{i\}
7:

S=S∪\{sj\}S=S\\cup\\\{s\_\{j\}\\\}
8:

di=min⁡\(di,∥pi−psj∥\),∀i∈V∖Sd\_\{i\}=\\min\\\!\\left\(d\_\{i\},\\,\\lVert p\_\{i\}\-p\_\{s\_\{j\}\}\\rVert\\right\),~\\forall i\\in V\\setminus S
9:endfor

10:Voronoi partitioning

11:

C⁡\(i\)=arg⁡mina∈\{1,…,k\}⁡hop⁡\(i,sa\),∀i∈VC\(i\)=\\arg\\min\_\{a\\in\\\{1,\\dots,k\\\}\}\\mathrm\{hop\}\(i,s\_\{a\}\),~\\forall i\\in V⊳\\trianglerightCC: Voronoi assignment map

12:⊳\\trianglerighthop⁡\(i,j\)\\mathrm\{hop\}\(i,j\): shortest path between nodesiiandjjin mesh edges

13:

𝒞a=\{i∈V:C⁡\(i\)=a\},∀a∈\{1,…,k\}\\mathcal\{C\}\_\{a\}=\\\{\\,i\\in V:C\(i\)=a\\,\\\},~\\forall a\\in\\\{1,\\dots,k\\\}
14:Coarse graph construction

15:

Vc=SV\_\{c\}=S⊳\\trianglerightVcV\_\{c\}: Coarse nodes

16:

Ec=\{\(C\(i\),C\(j\)\)\|\(i,j\)∈E,C\(i\)≠C\(j\)\}E\_\{c\}=\\left\\\{\(C\(i\),C\(j\)\)\\;\\middle\|\\;\(i,j\)\\in E,\\ C\(i\)\\neq C\(j\)\\right\\\}⊳\\trianglerightEcE\_\{c\}: Coarse edges

17:return

C,S⁡\(Vc,Ec\)C,\\ S\\ \(V\_\{c\},E\_\{c\}\)

Upsampling lifts the coarse latents back to the fine resolution using a learned, geometry\-aware GNN\. A coarse\-to\-fine graph connects each fine nodeViV\_\{i\}to its own coarse seed the neighboring coarse nodes of that seed, forming a set of coarse source nodes𝒞⁡\(i\)\\mathcal\{C\}\(i\)\. Each sourcec∈𝒞⁡\(i\)c\\in\\mathcal\{C\}\(i\)sends a message built from the coarse latentzcz\_\{c\}, the stored fine\-scale skip connectionsis\_\{i\}, and the reference\-configuration offset ofiifrom the seed anchor,

μi​c=ψ⁡\(\[zc,si,pi−psc\]\),\\mu\_\{ic\}=\\psi\(\[z\_\{c\},s\_\{i\},p\_\{i\}\-p\_\{s\_\{c\}\}\]\),\(7\)and the fine latent state is reconstructed from the aggregated messages together with the skip connection,

ziup=M​L​P​\(\[si,∑c∈𝒞⁡\(i\)μi​c\]\)\.z^\{\\mathrm\{up\}\}\_\{i\}=MLP\(\[s\_\{i\},\\sum\_\{c\\in\\mathcal\{C\}\(i\)\}\\mu\_\{ic\}\]\)\.\(8\)The coarsening, pooling, and interpolation operators are illustrated in Fig\.[2](https://arxiv.org/html/2608.13827#S2.F2), and the forward pass is summarized in Algorithm[2](https://arxiv.org/html/2608.13827#alg2)\.

![Refer to caption](https://arxiv.org/html/2608.13827v1/Downsample_Upsample.png)Figure 2:Coarsening and GNN interpolation in HI\-MGN\.Algorithm 2Hierarchical message\-passing processor of HI\-MGN\.1:Fine latent features

zz; graph edges

\{E\(ℓ\)\}ℓ=0H\\\{E^\{\(\\ell\)\}\\\}\_\{\\ell=0\}^\{H\}; coarsening maps

\{C\(ℓ\),S\(ℓ\)\}ℓ=0H−1\\\{C^\{\(\\ell\)\},S^\{\(\\ell\)\}\\\}\_\{\\ell=0\}^\{H\-1\}; message\-passing depths

\{Ls\}s=02​H\\\{L\_\{s\}\\\}\_\{s=0\}^\{2H\}
2:Downsampling

3:for

ℓ=0,…,H−1\\ell=0,\\dots,H\-1do

4:

z=MPLℓ​\(z,E\(ℓ\)\)z=\\mathrm\{MP\}^\{L\_\{\\ell\}\}\\bigl\(z,E^\{\(\\ell\)\}\\bigr\)
5:

zskip\(ℓ\)=zz\_\{\\mathrm\{skip\}\}^\{\(\\ell\)\}=z
6:

z=M​e​a​n​P​o​o​l​\(z,C\(ℓ\)\)z=MeanPool\(z,C^\{\(\\ell\)\}\)
7:endfor

8:Coarsest\-level propagation

9:

z=MPLH​\(z,E\(H\)\)z=\\mathrm\{MP\}^\{L\_\{H\}\}\\bigl\(z,E^\{\(H\)\}\\bigr\)
10:Upsampling

11:for

ℓ=H−1,…,0\\ell=H\-1,\\dots,0do

12:

z=Interpolate⁡\(z,zskip\(ℓ\),C\(ℓ\),S\(ℓ\)\)z=\\mathrm\{Interpolate\}\\bigl\(z,\\,z\_\{\\mathrm\{skip\}\}^\{\(\\ell\)\},\\,C^\{\(\\ell\)\},\\,S^\{\(\\ell\)\}\\bigr\)
13:

z=MPL2​H−ℓ​\(z,E\(ℓ\)\)z=\\mathrm\{MP\}^\{L\_\{2H\-\\ell\}\}\\bigl\(z,E^\{\(\\ell\)\}\\bigr\)
14:endfor

15:return

zz

### 2\.4Loss and training

All input and output channels are standardized to zero mean and unit variance using statistics computed from the training split\. The network minimizes a channel\-weighted Huber loss on the normalized physical\-state updates,

ℒ=1N​∑i∑qwq​Huberη​\(y^i,q−yi,q\),\\mathcal\{L\}=\\frac\{1\}\{N\}\\sum\_\{i\}\\sum\_\{q\}w\_\{q\}\\,\\mathrm\{Huber\}\_\{\\eta\}\\\!\\big\(\\hat\{y\}\_\{i,q\}\-y\_\{i,q\}\\big\),\(9\)with normalized feature weightswqw\_\{q\}summing to one for all of the output channelsqq\. The Huber loss is less sensitive to large residuals than the mean squared error because of its linear growth for residuals beyond the transition threshold, which can be beneficial for localized large errors near stress concentrations or contact regions\. Following the noise\-injection strategy used in MeshGraphNets\[[29](https://arxiv.org/html/2608.13827#bib.bib33)\], zero\-mean Gaussian noise \(σ=0\.01\\sigma=0\.01\) is added to the temporal physical\-state inputs during training\. The corresponding geometry\-dependent edge features are then recomputed from the perturbed nodal states, and the prediction targets are corrected consistently with the applied perturbation to improve autoregressive rollout robustness\. We further apply data augmentations, including random rotations and reflections when the governing physics and boundary conditions are preserved under the corresponding transformations\. During training, we use the AdamW optimizer with gradient clipping, learning\-rate warmup followed by cosine annealing, 16\-bit mixed\-precision training, and an exponential moving average \(EMA\) of the model weights\.

## 3Numerical examples

In this section, we compare the results of HI\-MGN with those of the original MeshGraphNets and BSMS\-GNN\[[29](https://arxiv.org/html/2608.13827#bib.bib33),[4](https://arxiv.org/html/2608.13827#bib.bib35)\]\. The number of total message passing blocks for is kept identical for HI\-MGN and MeshGraphNets\. BSMS\-GNN multiscale levels are set to values between 7–9, as the cases reported by the authors\. The hyperparameters of MeshGraphNets and BSMS\-GNN are also set to the default values set by the original authors except for the number of message passing blocks\. Through the three benchmarks, we show that our method performs consistently across the considered settings and outperforms the previously developed methods\. HI\-MGN not only achieves improved accuracy but also requires competent training time and peak memory usage\.

### 3\.12D static thermoelastic analysis

We consider a two\-dimensional multi\-material static thermoelastic problem\. Figure[3](https://arxiv.org/html/2608.13827#S3.F3)illustrates the configuration\. The domain comprises two parts, a backgroundAAand an inclusionBBwith different thermoelastic material properties, and is discretized using plane\-stress three\-node triangular elements\. The background is a rectangle subjected to a fixed displacement boundary condition and a prescribed temperature \(0​K0\\,\\mathrm\{K\}\) on its left edge, whereas a heat flux of100​W/m2100\\,\\mathrm\{W/m^\{2\}\}is imposed on its right edge\. The background dimensions vary across samples, characterized with a width of2​a2aand a height of2​b2b, wherea∈\[1000,3000\]​mma\\in\[1000,3000\]\\,\\mathrm\{mm\}andb∈\[1000,3000\]​mmb\\in\[1000,3000\]\\,\\mathrm\{mm\}\. The inclusion assumes two distinct geometries, a triangle and a pentagon\. The parameterized inclusions are governed by the parametersc∈\[200,700\]​mmc\\in\[200,700\]\\,\\mathrm\{mm\},d∈\[200,700\]​mmd\\in\[200,700\]\\,\\mathrm\{mm\},γc∈\[0\.2,0\.8\]\\gamma\_\{c\}\\in\[0\.2,0\.8\], andγd∈\[0\.2,0\.8\]\\gamma\_\{d\}\\in\[0\.2,0\.8\]\. For each inclusion geometry, 50 parameter sets are generated using Latin hypercube sampling\. The material properties are summarized in Table[1](https://arxiv.org/html/2608.13827#S3.T1)\.

![Refer to caption](https://arxiv.org/html/2608.13827v1/Figure_3.jpg)Figure 3:Schematics for 2D thermoelastic exampleTable 1:Material properties of the background and the inclusion\.For training, we consider a three\-level hierarchy\. LettingNNas the number of nodes in the training mesh, HI\-MGN is configured with the hierarchy,\[N,5000,100\]\[N,5000,100\], while BSMS\-GNN usesL=7L=7hierarchy and MeshGraphNets uses a single level,\[N\]\[N\]\. We set the number of message passing blocks to be identical across HI\-MGN and MeshGraphNets, so that HI\-MGN block configuration become \[4, 6, 8, 6, 4\], and that of MeshGraphNets is \[28\]\. The message passing blocks of BSMS\-GNN is 13 blocks spread out seven hierarchies\. The remaining hyperparameters are kept identical across the models: 5000 epochs, a batch size of 1, a latent dimension of 128, the use of an EMA model, and 16\-bit mixed\-precision training\. All models are trained on NVIDIA GeForce RTX 3090 GPU and an Intel Xeon Gold 5218R CPU\.

For inference, we consider geometric extrapolation\. Based on various training samples with triangular and pentagon inclusion, the considered GNN models are tested on a rectangular inclusion which is not included in the training geometries\. The schematic of the rectangular inclusion is shown in Fig\.[4](https://arxiv.org/html/2608.13827#S3.F4)and the corresponding inference results are shown in Fig\.[5](https://arxiv.org/html/2608.13827#S3.F5)\. TheR2R^\{2\}values between the ground truth and the considered GNNs are also shown in Table[2](https://arxiv.org/html/2608.13827#S3.T2)\.

![Refer to caption](https://arxiv.org/html/2608.13827v1/Figure_4.jpg)Figure 4:Schematic and finite element discretization for the rectangular inclusion inference, witha=b=2,000​mma=b=2,000\\,\\mathrm\{mm\}andc=d=450​mmc=d=450\\,\\mathrm\{mm\}\.![Refer to caption](https://arxiv.org/html/2608.13827v1/ex1/ex1_inference.png)Figure 5:Inferred results for geometric extrapolation with rectangular inclusion\. All units are in Pa \(stress\) and mm \(displacement\)\.Table 2:R2R^\{2\}between the ground truth and rectangular inclusion inference results of HI\-MGN, BSMS\-GNN, and MeshGraphNets\.The results show that HI\-MGN outperforms BSMS\-GNN and MeshGraphNets in terms of accuracy\. In particular, Fig\.[5](https://arxiv.org/html/2608.13827#S3.F5)shows that HI\-MGN shows considerably smoother and accurate prediction compared to BSMS\-GNN and MeshGraphNets\. Especially, MeshGraphNets accurately captures the responses near the boundaries while poorly estimating the central region\. Table[3](https://arxiv.org/html/2608.13827#S3.T3)shows computational cost breakdown of the considered GNNs\. Setting the original MeshGraphNets as the baseline, HI\-MGN is found to reduce peak VRAM usage by 38\.6% while reducing training time by 59\.9%\. BSMS\-GNN required smaller VRAM and training time as it required 9\.93GB for VRAM and 6\.48 hours for training time\.

Table 3:Training time and peak VRAM usage of HI\-MGN, BSMS\-GNN, and MeshGraphNets\.
### 3\.23D nonlinear contact analysis

Next, we consider a three\-dimensional nonlinear contact problem\. Figure[6](https://arxiv.org/html/2608.13827#S3.F6)illustrates the configuration, where the indenter \(AA\) moves downward to press the rectangular parallelepiped foundation block \(BB\)\. With the coordinate origin placed at the center of the bottom face of the foundation block, the configuration is symmetric with respect to theX​YXY\- andY​ZYZ\-planes\. The top surface of the indenter is prescribed with a displacement ofδy=−160​mm\\delta\_\{y\}=\-160\\,\\mathrm\{mm\}, which is applied over 50 steps, while the bottom surface of the foundation block is fixed\. The indenter and foundation block, designated as the slave and master bodies, are discretized with 4\-node tetrahedral elements and 8\-node hexahedral elements, respectively\. The indenter dimensions are parameterized and vary across samples, with a width of2​a2aand a contact surface radius ofR=γa​aR=\\gamma\_\{a\}a, wherea∈\[160,320\]​mma\\in\[160,320\]\\,\\mathrm\{mm\}andγa∈\[0\.125,0\.875\]\\gamma\_\{a\}\\in\[0\.125,0\.875\], which are selected via Latin hypercube sampling with 50 samples\. Both bodies consist of a compressible neo\-Hookean hyperelastic material, whose properties are summarized in Table[4](https://arxiv.org/html/2608.13827#S3.T4)\.

![Refer to caption](https://arxiv.org/html/2608.13827v1/Figure_6.jpg)Figure 6:Schematic of the 3D contact example\. All unspecified rounding fillet radii are30​mm30\\,\\mathrm\{mm\}, and all units are in mm\. The finite element discretization shown in isometric view corresponds to the example case ofa=240​mma=240\\,\\mathrm\{mm\}andγa=0\.5\\gamma\_\{a\}=0\.5\.Table 4:Material properties of the indenter and the foundation block\.To train the dataset, we consider the same three\-level hierarchy for HI\-MGN\. The message passing block configuration and coarse\-graph construction are all kept the same as those used in the 2D static thermoelastic analysis\. However, because the original MeshGraphNets study reported that 15 message\-passing blocks were effective for time\-stepping problems, we additionally evaluate a MeshGraphNet model with 15 message\-passing blocks\[[29](https://arxiv.org/html/2608.13827#bib.bib33)\]\. We train the models for 100 epochs and keep the remaining hyperparameters the same as those used in the 2D static thermoelastic analysis\.

For inference, we also consider geometric extrapolation\. A large indenter with minimal rounding for numerical stability is selected to test the models on stress concentration with large spatial gradient in the physics field\. A geometric extrapolation test case witha=340a=340mm andγa=0\.05\\gamma\_\{a\}=0\.05is selected as seen in Fig\.[7](https://arxiv.org/html/2608.13827#S3.F7)\. The inference results for this geometric\-extrapolation case are shown in Fig\.[8](https://arxiv.org/html/2608.13827#S3.F8)and Table[5](https://arxiv.org/html/2608.13827#S3.T5)\.

![Refer to caption](https://arxiv.org/html/2608.13827v1/Figure_7.jpg)Figure 7:Schematic and finite element discretization for the extrapolation sample of 3D contact problem, witha=340​mma=340\\,\\mathrm\{mm\}andγa=0\.05\\gamma\_\{a\}=0\.05\.![Refer to caption](https://arxiv.org/html/2608.13827v1/ex2/ex2_inference.png)Figure 8:Inferred results for geometric extrapolation\.Table 5:Time\-averagedR2R^\{2\}between the ground truth and inference results of HI\-MGN, BSMS\-GNN, and MeshGraphNets\.![Refer to caption](https://arxiv.org/html/2608.13827v1/ex2/ex2_R2.png)Figure 9:R2R^\{2\}results with respect to time\.Table 6:Training time and peak VRAM usage of HI\-MGN, BSMS\-GNN, and MeshGraphNets\.The results show that HI\-MGN also outperforms the existing methods in terms of prediction accuracy\. While the time\-averagedR2R^\{2\}shown in Table[5](https://arxiv.org/html/2608.13827#S3.T5)may appear relatively modest, but they are strongly influenced by the early stages of contact as seen in Fig\.[9](https://arxiv.org/html/2608.13827#S3.F9), where the magnitudes and spatial variations of the predicted physical fields are small, but greatly influencesR2R^\{2\}\. More importantly, shortly after the onset of contact \(t=12t=12\), both the 15\- and 28\-block MeshGraphNets exhibit prediction errors consistent with limited long\-range information propagation\. The contours in Fig\.[8](https://arxiv.org/html/2608.13827#S3.F8)illustrate the limitation of conventional flat message passing, whose receptive field is restricted by the number of message\-passing blocks\. BSMS\-GNN also exhibits noticeable errors in the bottom\-right region despite using 28 message\-passing blocks, whereas HI\-MGN maintains substantially better agreement with the ground\-truth solution\. Even with the later stage of rollout \(t=49t=49\), BSMS\-GNN and both MeshGraphNets variants continue to exhibit localized errors consistent with insufficient long\-range information propagation\.

In terms of computational efficiency, HI\-MGN also achieves the best overall performance\. Table[6](https://arxiv.org/html/2608.13827#S3.T6)summarizes the computational costs of the considered GNN models\. HI\-MGN reduces peak VRAM usage by up to approximately 48\.3% and training time by up to approximately 63\.0% compared with the other GNN models\. The reason for large computational burden in BSMS\-GNN is due to its bi\-stride contraction and edge building\. BSMS\-GNN halves the number of nodes at every coarsening level, but reconnects any two surviving nodes that were within two hops of each other, so the connection radius doubles at every level\. Whether this stays cheap depends on the mesh dimension d: halving the nodes spreads them apart by only a factor of2\(1/d\)2^\{\(1/d\)\}, so on a 2D surface the doubled radius is roughly matched by the wider spacing and the coarse graphs stay sparse, while on a 3D volume the added connectivity radius grows faster than the spacing and each coarse node reaches ever more neighbors — the coarse graphs become denser instead of smaller\. Our hierarchy avoids this entirely: coarse edges are created only by merging fine edges, so their number can never grow

### 3\.33D steady aerodynamic analysis

Lastly, we consider a 3D steady aerodynamic analysis of the NASA Common Research Model \(CRM\)\. The NASA CRM is a publicly available benchmark configuration representing a scaled transport aircraft and has been widely recognized in the aerodynamics community\. Here, we adopt the simulation results generated by DLR using the DLR TAU Reynolds\-Averaged Navier–Stokes \(RANS\) solver\[[3](https://arxiv.org/html/2608.13827#bib.bib31)\]\. The dataset consists of steady\-state pressure and surface friction results under various flight conditions of the NASA CRM\. Computational fluid dynamics \(CFD\) results are provided for six varying input parameters, including the freestream Mach number \(M​aMa\), angle of attack \(AoA\), inboard and outboard aileron deflections, elevator deflection, and horizontal tail plane angle\. The original dataset consists of 454,404 surface nodes with the pressure coefficient \(CpC\_\{p\}\), and surface friction coefficients inxx,yy, andzzdirections \(Cfx,Cfy,CfzC\_\{f\_\{x\}\},C\_\{f\_\{y\}\},C\_\{f\_\{z\}\}\)\. Further information regarding the NASA CRM dataset can be found in the work of Bekemeyer et al\.\[[3](https://arxiv.org/html/2608.13827#bib.bib31)\]\.

In this section, we train on a subsampled version of the original DLR dataset due to the large memory requirements of MeshGraphNets\. Based on the DLR results, the original surface nodal data of 454,404 nodes are subsampled to 122,778 nodes with reconstructed edges\. The dataset contains 105 flight conditions for the training and 44 conditions for testing\. We train the models using the same hyperparameters as in the previous numerical examples, except for the BSMS\-GNN hierarchical level set toL=9L=9and the training epochs set to 1,000\. As in the previous tests, HI\-MGN and BSMS\-GNN are trained using three\-level hierarchies, while MeshGraphNets is evaluated using 15 and 28 message\-passing blocks\.

For inference, the considered GNN methods are evaluated on the 44 test cases\. Results for the first sample in the test set are shown in Fig\.[10](https://arxiv.org/html/2608.13827#S3.F10), corresponding toM​aMa: 0\.72, AoA:3\.45∘3\.45^\{\\circ\}, inboard ailerons:−12\.41∘\-12\.41^\{\\circ\}, outboard ailerons:1\.36∘1\.36^\{\\circ\}, elevator:−2\.89∘\-2\.89^\{\\circ\}, and horizontal tailplane:0\.95∘0\.95^\{\\circ\}\. TheR2R^\{2\}values averaged over all test samples are also shown in Table[7](https://arxiv.org/html/2608.13827#S3.T7)\.

![Refer to caption](https://arxiv.org/html/2608.13827v1/ex3/ex3_Cp.png)

![Refer to caption](https://arxiv.org/html/2608.13827v1/ex3/ex3_Cp_wing.png)

Figure 10:InferredCpC\_\{p\}results for the first sample in the test set\.Table 7:Sample\-averagedR2R^\{2\}between the ground truth and inference results of HI\-MGN, BSMS\-GNN, and MeshGraphNets\.Table 8:Training time and peak VRAM usage of HI\-MGN, BSMS\-GNN, and MeshGraphNets\.The results shown in Fig\.[10](https://arxiv.org/html/2608.13827#S3.F10)and Table[7](https://arxiv.org/html/2608.13827#S3.T7)indicate that all considered models achieve high prediction accuracy with only minor differences among them\. HI\-MGN achieves the highestR2R^\{2\}values forCp,Cfx,CfyC\_\{p\},C\_\{f\_\{x\}\},C\_\{f\_\{y\}\}although the improvements over the other methods are relatively small\. ForCfzC\_\{f\_\{z\}\}, MeshGraphNets with 28 message\-passing blocks marginally outperform HI\-MGN\. These results suggest that all four models provide comparable accuracy for this steady aerodynamic problem, with HI\-MGN showing a small overall advantage across most predicted quantities\. In contrast, the computational\-cost differences are more pronounced\. As summarized in Table[8](https://arxiv.org/html/2608.13827#S3.T8), HI\-MGN requires the lowest peak VRAM usage while BSMS\-GNN required shorter training time among the considered models\. Such phenomenon occurs as the dataset consists of only the 3D surface mesh, making the dimension similar to 2D mesh\. Compared with the other GNN models, HI\-MGN reduces peak VRAM usage by up to 51\.3% and training time by up to approximately 58\.5%\.

## 4Conclusions

In this paper, we present a new hierarchical extension of MeshGraphNets\. The model, HI\-MGN, is designed to alleviate limited long\-range information propagation in graph neural networks\. The developed HI\-MGN is a self\-contained graph neural network that performs hierarchical message passing across multiple graph resolutions\. The model selects coarse nodes using FPS, assigns fine nodes through graph\-based Voronoi partitioning, and constructs coarse edges from the connectivity of the original mesh\. This construction preserves the underlying mesh topology while avoiding proximity\-based connections between disconnected components\. Upsampling from the coarse to the fine level is performed by a learned geometry\-aware GNN interpolator that combines neighboring coarse latent states, fine\-scale skip features, and relative geometric information to reconstruct the fine\-resolution representation\.

HI\-MGN is compared against the original MeshGraphNets and its multiscale variant, BSMS\-GNN, across static structural, transient contact, and steady aerodynamic problems\. The numerical results show that the advantages of hierarchical message passing are particularly pronounced when the prediction requires information to propagate over large graph distances\. In the static thermoelastic problem, the solution is globally coupled to spatially localized boundary conditions and must be reconstructed within a single network evaluation\. Consequently, the limited receptive field of flat message passing leads to larger errors away from the boundaries, whereas the coarse graphs of HI\-MGN provide substantially shorter communication paths\. In the contact problem, the performance differences are smaller but remain visible, particularly shortly after the onset of contact, when newly generated local contact information must propagate rapidly through the domain\. For the steady aerodynamic problem, all considered models achieve high prediction accuracy, while HI\-MGN maintains comparable or slightly improved accuracy for most quantities\.

Across the three benchmarks, HI\-MGN consistently provides improved accuracy accross all scenarios\. For the computational burden, HI\-MGN showed competitive computational requirements in 2D and surface meshes while reducing significant amount of memory and training time for 3D meshes\. The aggressive hierarchical coarsening allows a large portion of the message passing to be performed on substantially smaller graphs, reducing the number of high\-resolution latent node and edge states that must be retained during training\. At the same time, topology\-aware coarse\-edge construction and learned coarse\-to\-fine interpolation allow the model to retain fine\-scale geometric information while improving communication over long graph distances\. In the considered numerical examples, these properties reduce peak VRAM usage by up to 51\.3% and training time by 63\.0% relative to the more computationally demanding baselines\.

Overall, the results demonstrate that hierarchical topology\-aware message passing provides an effective approach for mitigating the long\-range communication and computational and memory limitations of conventional mesh\-based GNNs\. Rather than increasing the depth of a flat processor, HI\-MGN reduces graph communication distances through aggressive coarsening and reconstructs the fine\-scale representation through learned geometric interpolation\. These characteristics make HI\-MGN a practical framework for unstructured engineering meshes involving complex geometries, multiple components, deformation, and long\-range physical interactions\.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper\.

## Data availability

The source code for HI\-MGN, including the graph\-coarsening and interpolation modules, training configurations, benchmark preprocessing scripts, and scripts used to reproduce the main results of this study, is publicly available at https://github\.com/leesihun/MeshGraphNets\.

## Appendix AComputational cost evaluation

To distinguish model complexity from implementation dependent wall\-clock cost, we report the number of trainable parameters, end\-to\-end preprocessing time, and inference time in addition to peak GPU memory and total training time\. For hierarchical models, preprocessing time includes the construction of the coarse graphs and inter\-level mappings\. We separately report model\-only latency, for which all graph hierarchies are precomputed\.

Table 9:Additional information regarding the 2D static thermoelastic analysisTable 10:Additional information regarding the 3D nonlinear contact analysisTable 11:Additional information regarding the 3D steady aerodynamic analysis
## Appendix BParameter matching BSMS\-GNN

The computational cost shown in Appendix A shows that number of parameters are considerably smaller for BSMS\-GNN\. While this is because the original authors intended to use BSMS\-GNN with only one message\-passing for each hierarchy\. For that, we also test BSMS\-GNN with matching number of parameters for the 2D static thermoelastic analysis\. BSMS\-GNN is trained with 192 latent dimension with other hyperparameters fixed as stated in the manuscript\. The number of parameters, peak VRAM, training time, andR2R^\{2\}results are shown in Table[12](https://arxiv.org/html/2608.13827#A2.T12)\.

Table 12:Numerical results of parameter matching BSMS\-GNN on 2D static thermoelastic analysis\.While the stress and displacment accuracy has enhanced compared to 128 latent dimension BSMS\-GNN, HI\-MGN still performs better even with less number of parameters\. Instead, VRAM and train time has risen for BSMS\-GNN, making HI\-MGN more efficient in both the train time and peak memory usage\. Here, we show that the number of parameters is not the cause of accuracy discrepancy of HI\-MGN and BSMS\-GNN\.

## Appendix CAblation study

For the sake of clarity and readability of the paper, we perform the ablation study on 2D static thermoelastic analysis, the case in which the largest variety in accuracy is observed\. All variants share the same encoder, decoder, training schedule, data split, and loss; one design axis is varied at a time\.R2R^\{2\}is reported as the average over the two output channels, displacement and stress\. The baseline is a two\-stage hierarchy with\[5,000,100\]\[5,000,100\]clusters, coarse\-centric message passing blocks of\[4,6,8,6,4\]\[4,6,8,6,4\], FPS\-Voronoi latent mean pooling and learned GNN interpolator for upsampling used throughout the manuscript; each block of rows varies one axis of it\.

Table 13:Ablation study on 2D static thermoelastic analysis\.The ablation results show that HI\-MGN is relatively insensitive to the detailed choice of message\-passing parameters\. Changing the distribution of message\-passing blocks from the baseline \(\[4,6,8,6,4\]\) to either a flat \(\[5,6,6,6,5\]\) or fine\-centric \(\[7,5,4,5,7\]\) configuration results in only moderate changes in accuracy\. Similarly, reducing the total number of message\-passing blocks from 28 to 14 retains anR2R^\{2\}of 0\.9920, while increasing it to 38 does not improve the prediction accuracy\. These results indicate that simply increasing the number of message\-passing operations, or carefully tuning their distribution across resolutions, has a relatively limited effect on the overall performance\.

In contrast, larger differences are observed when the hierarchical representation itself is modified\. Extending the hierarchy to three stages reduces theR2R^\{2\}to 0\.9862, while replacing the latent mean pooling with direct FPS\-seed latent inheritance results in anR2R^\{2\}of 0\.9874\. Replacing the learned GNN interpolation with linear interpolation also decreases the accuracy\. Overall, these results indicate that HI\-MGN does not rely on extensive tuning of the message\-passing parameters to achieve high accuracy\. The performance remains relatively stable over substantial changes in both the number and distribution of message\-passing blocks, whereas modifications to the hierarchical representation and the inter\-level transfer operators lead to more noticeable degradation\. This suggests that the primary source of accuracy is the hierarchical multi\-scale structure itself, which enables effective information propagation across different graph resolutions\. The detailed processor configuration plays a secondary role once this hierarchical communication path is established\.

## References

- \[1\]\(2025\)AB\-UPT: scaling neural CFD surrogates for high\-fidelity automotive aerodynamics simulations via anchored\-branched universal physics transformers\.arXiv preprint arXiv:2502\.09692\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p3.1)\.
- \[2\]P\. W\. Battaglia, J\. B\. Hamrick, V\. Bapst, A\. Sanchez\-Gonzalez, V\. Zambaldi, M\. Malinowski, A\. Tacchetti, D\. Raposo, A\. Santoro, R\. Faulkner,et al\.\(2018\)Relational inductive biases, deep learning, and graph networks\.arXiv preprint arXiv:1806\.01261\.Cited by:[§2\.2](https://arxiv.org/html/2608.13827#S2.SS2.p1.1)\.
- \[3\]P\. Bekemeyer, N\. Hariharan, A\. M\. Wissink, and J\. Cornelius\(2025\)Introduction of applied aerodynamics surrogate modeling benchmark cases\.InAIAA Scitech 2025 Forum,External Links:[Document](https://dx.doi.org/10.2514/6.2025-0036)Cited by:[§3\.3](https://arxiv.org/html/2608.13827#S3.SS3.p1.1)\.
- \[4\]Y\. Cao, M\. Chai, M\. Li, and C\. Jiang\(2023\)Efficient learning of mesh\-based physical simulation with bi\-stride multi\-scale graph neural network\.InProceedings of the 40th International Conference on Machine Learning,A\. Krause, E\. Brunskill, K\. Cho, B\. Engelhardt, S\. Sabato, and J\. Scarlett \(Eds\.\),Proceedings of Machine Learning Research, Vol\.202,pp\. 3541–3558\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p5.1),[§3](https://arxiv.org/html/2608.13827#S3.p1.1)\.
- \[5\]M\. Fortunato, T\. Pfaff, P\. Wirnsberger, A\. Pritzel, and P\. Battaglia\(2022\)Multiscale meshgraphnets\.arXiv preprint arXiv:2210\.00612\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p5.1)\.
- \[6\]F\. J\. Gonzalez and M\. Balajewicz\(2018\)Deep convolutional recurrent autoencoders for learning low\-dimensional feature dynamics of fluid systems\.arXiv preprint arXiv:1808\.01346\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[7\]J\. He, S\. Koric, D\. Abueidda, A\. Najafi, and I\. Jasiuk\(2024\)Geom\-DeepONet: a point\-cloud\-based deep operator network for field predictions on 3D parameterized geometries\.Comput\. Methods Appl\. Mech\. Eng\.429,pp\. 117130\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p3.1)\.
- \[8\]J\. S\. Hesthaven and S\. Ubbiali\(2018\)Non\-intrusive reduced order modeling of nonlinear problems using neural networks\.J\. Comput\. Phys\.363,pp\. 55–78\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2018.02.037)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[9\]M\. M\. Iparraguirre, I\. Alfaro, D\. Gonzalez, and E\. Cueto\(2026\)MeshGraphNet\-transformer: scalable mesh\-based learned simulation for solid mechanics\.arXiv preprint arXiv:2601\.23177\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p6.1)\.
- \[10\]T\. Kadeethum, F\. Ballarin, Y\. Choi, D\. O’Malley, H\. Yoon, and N\. Bouklas\(2022\)Non\-intrusive reduced order modeling of natural convection in porous media using convolutional autoencoders: comparison with linear subspace techniques\.Adv\. Water Resour\.160,pp\. 104098\.External Links:[Document](https://dx.doi.org/10.1016/j.advwatres.2021.104098)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[11\]T\. Kadeethum, F\. Ballarin, D\. O’Malley, Y\. Choi, N\. Bouklas, and H\. Yoon\(2022\)Reduced order modeling with barlow twins self\-supervised learning: navigating the space between linear and nonlinear solution manifolds\.arXiv preprint arXiv:2202\.05460\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[12\]Y\. Kang, S\. Yang, and K\. Yee\(2022\)Physics\-aware reduced\-order modeling of transonic flow viaβ\\beta\-variational autoencoder\.Phys\. Fluids34\(7\),pp\. 076103\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[13\]J\. Kneifl, D\. Grunert, and J\. Fehr\(2021\)A nonintrusive nonlinear model reduction method for structural dynamical problems based on machine learning\.Int\. J\. Numer\. Methods Eng\.122\(17\),pp\. 4774–4786\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[14\]K\. Lee and K\. T\. Carlberg\(2020\)Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders\.J\. Comput\. Phys404,pp\. 108973\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[15\]S\. Lee, K\. Jang, H\. Cho, H\. Kim, and S\. J\. Shin\(2021\)Parametric non\-intrusive model order reduction for flow\-fields using unsupervised machine learning\.Comput\. Methods Appl\. Mech\. Eng\.384,pp\. 113999\.External Links:[Document](https://dx.doi.org/10.1016/j.cma.2021.113999)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[16\]S\. Lee, K\. Jang, S\. Lee, H\. Cho, and S\. J\. Shin\(2024\)Parametric model order reduction by machine learning for fluid–structure interaction analysis\.Eng\. Comput\.40,pp\. 45–60\.External Links:[Document](https://dx.doi.org/10.1007/s00366-023-01782-2)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[17\]S\. Lee, S\. Lee, K\. Jang, H\. Cho, and S\. J\. Shin\(2024\)Data\-driven nonlinear parametric model order reduction framework using deep hierarchical variational autoencoder\.Eng\. Comput\.40,pp\. 2385–2400\.External Links:[Document](https://dx.doi.org/10.1007/s00366-023-01916-6)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[18\]S\. Lee, S\. Kang, S\. Lee, and S\. Shin\(2025\)Physics\-aware neural network\-based parametric model\-order reduction of the electromagnetic analysis for a coated component\.Eng\. Comput\.41,pp\. 785–799\.External Links:[Document](https://dx.doi.org/10.1007/s00366-024-02056-1)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[19\]Q\. Li, Z\. Han, and X\. Wu\(2018\)Deeper insights into graph convolutional networks for semi\-supervised learning\.InProceedings of the AAAI conference on artificial intelligence,Vol\.32\.Cited by:[§2\.3](https://arxiv.org/html/2608.13827#S2.SS3.p1.1)\.
- \[20\]T\. Li, S\. Deng, K\. Zhang, H\. Wei, R\. Wang, J\. Fan, J\. Xin, and J\. Yao\(2021\)A nonintrusive parametrized reduced\-order model for periodic flows based on extended proper orthogonal decomposition\.Int\. J\. Comput\. Methods18\(9\),pp\. 2150035\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[21\]Z\. Li, D\. Z\. Huang, B\. Liu, and A\. Anandkumar\(2023\)Fourier neural operator with learned deformations for pdes on general geometries\.J\. Mach\. Learn\. Res\.24\(388\),pp\. 1–26\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p3.1)\.
- \[22\]Z\. Li, N\. Kovachki, C\. Choy, B\. Li, J\. Kossaifi, S\. Otta, M\. A\. Nabian, M\. Stadler, C\. Hundt, K\. Azizzadenesheli,et al\.\(2023\)Geometry\-informed neural operator for large\-scale 3D PDEs\.InAdvances in Neural Information Processing Systems,Vol\.36,pp\. 35836–35854\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p3.1)\.
- \[23\]L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, and G\. E\. Karniadakis\(2021\)Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators\.Nat\. Mach\. Intell\.3,pp\. 218–229\.External Links:[Document](https://dx.doi.org/10.1038/s42256-021-00302-5)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p3.1)\.
- \[24\]H\. Luo, H\. Wu, H\. Zhou, L\. Xing, Y\. Di, J\. Wang, and M\. Long\(2025\)Transolver\+\+: an accurate neural solver for pdes on million\-scale geometries\.arXiv preprint arXiv:2502\.02414\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p3.1)\.
- \[25\]A\. T\. Mohan and D\. V\. Gaitonde\(2018\)A deep learning based approach to reduced order modeling for turbulent flow control using lstm neural networks\.arXiv preprint arXiv:1804\.09269\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[26\]A\. Moosavi, R\. Ştefănescu, and A\. Sandu\(2018\)Multivariate predictions of local reduced\-order\-model errors and dimensions\.Int\. J\. Numer\. Methods Eng\.113\(3\),pp\. 512–533\.External Links:[Document](https://dx.doi.org/10.1002/nme.5624)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[27\]M\. A\. Nabian, C\. Liu, R\. Ranade, and S\. Choudhry\(2024\)X\-meshgraphnet: scalable multi\-scale graph neural networks for physics simulation\.arXiv preprint arXiv:2411\.17164\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p5.1)\.
- \[28\]A\. Peyvan, V\. Kumar, and G\. E\. Karniadakis\(2026\)Fusion\-DeepONet: a data\-efficient neural operator for geometry\-dependent hypersonic and supersonic flows\.J\. Comput\. Phys544,pp\. 114432\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p3.1)\.
- \[29\]T\. Pfaff, M\. Fortunato, A\. Sanchez\-Gonzalez, and P\. W\. Battaglia\(2020\)Learning mesh\-based simulation with graph networks\.arXiv preprint arXiv:2010\.03409\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p4.1),[§2\.1](https://arxiv.org/html/2608.13827#S2.SS1.p1.1),[§2\.2](https://arxiv.org/html/2608.13827#S2.SS2.p1.1),[§2\.4](https://arxiv.org/html/2608.13827#S2.SS4.p1.2),[§3\.2](https://arxiv.org/html/2608.13827#S3.SS2.p2.1),[§3](https://arxiv.org/html/2608.13827#S3.p1.1)\.
- \[30\]C\. R\. Qi, H\. Su, K\. Mo, and L\. J\. Guibas\(2017\)PointNet: deep learning on point sets for 3d classification and segmentation\.In2017 IEEE Conference on Computer Vision and Pattern Recognition \(CVPR\),pp\. 77–85\.External Links:[Document](https://dx.doi.org/10.1109/CVPR.2017.16)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p2.1)\.
- \[31\]C\. R\. Qi, L\. Yi, H\. Su, and L\. Guibas\(2017\)PointNet\+\+: deep hierarchical feature learning on point sets in a metric space\.InAdvances in Neural Information Processing Systems,I\. Guyon, U\. V\. Luxburg, S\. Bengio, H\. Wallach, R\. Fergus, S\. Vishwanathan, and R\. Garnett \(Eds\.\),Vol\.30\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p2.1)\.
- \[32\]A\. Sanchez\-Gonzalez, J\. Godwin, T\. Pfaff, R\. Ying, J\. Leskovec, and P\. Battaglia\(2020\)Learning to simulate complex physics with graph networks\.InProceedings of the 37th International Conference on Machine Learning,H\. D\. III and A\. Singh \(Eds\.\),Proceedings of Machine Learning Research, Vol\.119,pp\. 8459–8468\.Cited by:[§2\.2](https://arxiv.org/html/2608.13827#S2.SS2.p1.1)\.
- \[33\]A\. Solera\-Rico, C\. S\. Vila, M\. Gómez, Y\. Wang, A\. Almashjary, S\. Dawson, and R\. Vinuesa\(2023\)β\\beta\-Variational autoencoders and transformers for reduced\-order modelling of fluid flows\.arXiv preprint arXiv:2304\.03571\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[34\]J\. Wang, C\. He, R\. Li, H\. Chen, C\. Zhai, and M\. Zhang\(2021\)Flow field prediction of supercritical airfoils via variational autoencoder based deep learning framework\.Phys\. Fluids33\(8\),pp\. 086108\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[35\]S\. Wiewel, M\. Becher, and N\. Thuerey\(2019\)Latent space physics: towards learning the temporal evolution of fluid flow\.Comput\. Graph\. Forum38\(2\),pp\. 71–82\.External Links:[Document](https://dx.doi.org/10.1111/cgf.13620)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[36\]H\. Wu, H\. Luo, H\. Wang, J\. Wang, and M\. Long\(2024\)Transolver: a fast transformer solver for PDEs on general geometries\.arXiv preprint arXiv:2402\.02366\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p3.1)\.
- \[37\]D\. Xiao\(2019\)Error estimation of the parametric non\-intrusive reduced order model using machine learning\.Comput\. Methods Appl\. Mech\. Eng\.355,pp\. 513–534\.External Links:[Document](https://dx.doi.org/10.1016/j.cma.2019.06.018)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[38\]F\. Xiong, L\. Zhang, X\. Hu, and C\. Ren\(2023\)A point cloud deep neural network metamodel method for aerodynamic prediction\.Chin\. J\. Aeronaut\.36\(4\),pp\. 92–103\.External Links:[Document](https://dx.doi.org/10.1016/j.cja.2022.11.025)Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p2.1)\.
- \[39\]J\. Xu and K\. Duraisamy\(2020\)Multi\-level convolutional autoencoder networks for parametric prediction of spatio\-temporal dynamics\.Comput\. Methods Appl\. Mech\. Eng\.372,pp\. 113379\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p1.1)\.
- \[40\]Y\. Yu, J\. Choi, J\. Park, K\. Lee, and N\. Park\(2025\)PIORF: physics\-informed Ollivier\-Ricci flow for long\-range interactions in mesh graph neural networks\.arXiv preprint arXiv:2504\.04052\.Cited by:[§1](https://arxiv.org/html/2608.13827#S1.p6.1)\.