Adaptive Swarm Mesh Refinement using Deep Reinforcement Learning with Local Rewards

Niklas Freymuth^1correspondence to
Philipp Dahlinger^1
Tobias Würth^2
Simon Reisch^1

Luise Kärger^2
Gerhard Neumann^1

^1Autonomous Learning Robots,
Karlsruhe Institute of Technology,

^2Institute of Vehicle Systems Technology,
Karlsruhe Institute of Technology, Karlsruhe

Simulating physical systems is essential in engineering, but analytical solutions are limited to straightforward problems.
Consequently, numerical methods like the Finite Element Method (FEM) are widely used.
However, the FEM becomes computationally expensive as problem complexity and accuracy demands increase.
Adaptive Mesh Refinement (AMR) improves the FEM by dynamically allocating mesh elements on the domain, balancing computational speed and accuracy.
Classical AMR depends on heuristics or expensive error estimators, limiting its use in complex simulations.
While learning-based AMR methods are promising, they currently only scale to simple problems.
In this work, we formulate AMR as a system of collaborating, homogeneous agents that iteratively split into multiple new agents.
This agent-wise perspective enables a spatial reward formulation focused on reducing the maximum mesh element error.
Our approach, Adaptive Swarm Mesh Refinement (ASMR), offers efficient, stable optimization and generates highly adaptive meshes at user-defined resolution during inference.
Extensive experiments, including volumetric meshes and Neumann boundary conditions, demonstrate that ASMR exceeds heuristic approaches and learned baselines, matching the performance of expensive error-based oracle AMR strategies.
ASMR additionally generalizes to different domains during inference, and produces meshes that simulate up to $2$ orders of magnitude faster than uniform refinements in more demanding settings.

§ INTRODUCTION

Across industries and disciplines, carefully tuning or evaluating a technical system is essential for of most engineering tasks.
While this procedure historically relied on costly prototype testing, advances in virtual simulations over the last decades have provided a cost-effective alternative that closely mimics real systems.
These simulations utilize fundamental principles like mass, momentum, and energy conservation, which are often expressed through complex pde.
Since analytical solutions of such pde are limited to simple cases, numerical approximations, particularly the fem are commonly employed [Brenner and Scott, 2008, Reddy, 2019, Anderson et al., 2021].
The fem partitions the continuous problem domain into a mesh consisting of smaller, finite elements, allowing for an efficient numerical solution whose accuracy depends on the number used elements.
However, as the physics becomes more complex, accurate simulations demand significantly more elements, quickly rendering simulations prohibitively expensive.

A popular extension of the fem is amr, which dynamically allocates more mesh elements to regions of high solution variability, striking a favorable balance between computational efficiency and accuracy [Plewa et al., 2005, Huang and Russell, 2010, Fidkowski and Darmofal, 2011].
As amr allows the creation of a simulation-specific mesh, it has become increasingly important for complex problems across domains like fluid dynamics [Berger and Colella, 1989, Baker, 1997, Zhang et al., 2020, Wallwork et al., 2022], structural mechanics [Ortiz and Quigley Iv, 1991, Stein, 2007, Gibert et al., 2019], and astrophysics [Cunningham et al., 2009, Bryan et al., 2014, Guillet et al., 2019].
Yet, general-purpose amr methods based on heuristics [Zienkiewicz and Zhu, 1992], or computationally expensive error estimates [Bangerth and Rannacher, 2013] are generally limited in their efficiency and adaptability
[Mukherjee, 1996, Kita and Kamiya, 2001, Yano and Darmofal, 2012, Cerveny et al., 2019, Wallwork, 2021], complicating their effective use in practical scenarios.

As an emerging alternative, a body of work employs rl [Sutton and Barto, 2018] to frame amr as a sequential decision-making process [Yang et al., 2023, Yang et al., 2023, Foucart et al., 2023].
These rlamr methods encode the state of the current simulation, including an intermediate mesh, as an observation that is fed to one or many rl agents, who then determine which mesh elements to refine.
However, previous work currently only scales to simple problems due to an expensive inference process [Yang et al., 2023], misaligned objectives and high variance in the state transitions [Foucart et al., 2023], or noisy reward signals during training [Yang et al., 2023].
We instead formulate amr as a Swarm rl [Šošić et al., 2017, Hüttenrauch et al., 2019] problem, extending existing frameworks to a shared observation space and per-agent, spatial rewards.
Crucially, we allow for agents to split into new agents over time to model element sub-division in the refinement process, mapping the per-agent reward signal over refinement steps to assign credits for individual agents throughout an episode.
We use a mpn [Sanchez-Gonzalez et al., 2020], a class of gnn [Scarselli et al., 2009, Bronstein et al., 2021], for our policy due to their effectiveness in learning physical simulations [Pfaff et al., 2021, Brandstetter et al., 2022].
Our method, Adaptive Swarm Mesh Refinement (ASMR), consistently produces highly efficient mesh refinements with thousands of elements and can be applied to arbitrary pde.

This paper extends previously published work on ASMR [Freymuth et al., 2023].
While the ASMR policies could produce meshes at a pre-defined granularity, we now use an adaptive element penalty which is also provided as context information for the policy to  allow a single policy to produce meshes of different granularity.
Additional extensions include an improved mapping between agents over time that acts as a regularizer and leads to better results, an improved network architecture and an improved reward formulation based on the reduction of the maximum local error of the mesh.
While this reward formulation was already introduced and ablated in our prior work with inconclusive results, we now present a large scale study that shows its improved performance compared to the previously used reward.
We refer to the improved version as method, and to the previous version as asmr.
Experimentally, we find that the proposed changes substantially improve the method on more difficult task setups.
We further provide additional visualizations of our method, including a schematic of the agent mapping and a qualitative error comparison to a uniform mesh.
Finally, we add two challenging tasks that utilize Neumann boundary conditions and a $3$-dimensional domain to showcase the applicability of our method to a wider range of applications.
Figure <ref> provides a schematic overview of a single method refinement step.

We validate the effectiveness of method on a wide range of pde that require complex refinement strategies.
The experiments use static meshes with conforming triangular or tetrahedral elements and corresponding h-adaptive refinements [Arnold et al., 2000, Stevenson, 2008], i.e., refinements via element subdivision, due to their prevalence in engineering [Ho-Le, 1988, Jones and Plassmann, 1997, Nagarajan and Soghrati, 2018].
We consider several recent rlamr methods as baselines [Yang et al., 2023, Yang et al., 2023, Foucart et al., 2023].
As these methods have been shown to work for shallow mesh refinement and coarsening on dynamic tasks, we adapt them to our setting of static meshes that require multiple precise refinement steps and thousands of elements.
We further compare to a threshold-based refinement heuristic that uses the popular zz estimate [Zienkiewicz and Zhu, 1992] and different oracle error estimates based on a fine uniform mesh to mark elements for refinement.
We implement all tasks as OpenAI gym [Brockman et al., 2016] environments
and publish our code, including all approaches and tasks presented in this paper, at <https://github.com/NiklasFreymuth/ASMRplusplus>.

Our experimental results show that method learns to predict stable and consistent refinements across a wide range of tasks, while the baselines rlamr methods struggle to produce high-quality meshes.
Our method generally outperforms the heuristic when using the zz error estimate, and even compares well when using oracle error estimates while being $2$ to $100$ times faster to compute than the uniform mesh needed for the oracle information.
A series of additional experiments demonstrates that our method robustly generalizes to different domains and initial conditions, and that it can up-scale to larger domains and meshes during inference when using data-augmentation during training.
We show the importance of our per-agent reward function and its combination with our proposed mapping between agents over time on a series of ablation experiments.

§ RELATED WORK

§.§ Learned Physics Simulation and Graph Network Simulators.

Classical simulators are often highly problem-specific and quickly become prohibitively expensive for complex problems, especially when considering iterative optimization or parametric studies [Wanner and Hairer, 1996, Zimmerling et al., 2022, Brandstetter et al., 2022].
As such, a recent trend in machine learning aims to directly learn the simulation of physical systems from data using deep neural networks.
These approaches usually require an underlying ground-truth simulator to generate large amounts of training data, and use this data to predict quantities of interest for a given system of equations.

Early work on learned simulators uses simple feed-forward [Um et al., 2018, Zimmerling et al., 2019] or convolutional neural networks [Guo et al., 2016, Zimmerling et al., 2019, Zimmerling et al., 2022], which are ill-suited for, e.g., mesh-based representations in the fem.
Here, gns [Sanchez-Gonzalez et al., 2020, Pfaff et al., 2021] have emerged as an alternative architecture acting on graphs.
gns encode the simulation state as a graph, and utilize mpn, a type of gnn [Scarselli et al., 2009, Bronstein et al., 2021] to compute per-node quantities.
Compared to their classical counterparts, learned gns are orders of magnitude faster and fully differentiable [Battaglia et al., 2018, Pfaff et al., 2021], making them applicable to inverse design problems [Allen et al., 2022] and real-world planning of robotic deformable object manipulation [Shi et al., 2023].
They have recently been combined with a block-structured amr strategy that uses the intermediate refinement steps as coarser graph representations for the latent gns layers [Perera and Agrawal, 2024] to facilitate learning multiscale interactions.

As an alternative to these data-driven models, physics-informed neural networks [Raissi et al., 2019, Würth et al., 2023] directly optimize a neural network to predict simulations that satisfy the governing equations of a physical system.
While usually mesh-free, recent extensions combine these approaches with gns and meshes to facilitate generalization to new domains during inference [Würth et al., 2024].

All of the above methods share a common goal of using recent advances in deep learning to speed up the simulation of complex physical systems.
Yet, they do so by directly approximating the simulation, causing any prediction error of the learned model to directly affect the simulated quantities.
In this work, we instead learn an amr strategy to propose an efficient mesh for a classical solver, implicitly speeding up the simulation by reducing the size of the system of equations that must be solved.
This approach is more robust and risk-averse than directly learning the simulation, especially when considering complex physics problems that demand high precision and accuracy [Zhang et al., 2020, Huang et al., 2021, Wallwork et al., 2022].

§.§ Supervised Learning and Adaptive Mesh Refinement.

Several recent methods apply supervised learning to improve existing or devise new amr strategies.
Examples include using an mlp [Zhang et al., 2020] to directly calculate an error per element on a coarse mesh, subsequently refining or re-generating the mesh based on this error.
Extensions instead use a two-grid paradigm self-adaptive amr algorithm [Rachowicz et al., 1989, Rachowicz et al., 2006] to collect refinement labels for an mlp-based classifier acting on individual mesh elements [Służalec et al., 2023].
Another approach uses a convolutional neural network to predict local mesh densities from an image-based domain representation for fluid dynamics problems [Huang et al., 2021].
Other work uses recurrent networks to find optimal marking strategies for second-order elliptical pde [Bohn and Feischl, 2021], or employs neural networks to solve the strong form of the adjoint problem and use hand-crafted features to compute error estimates directly [Roth et al., 2022, Wallwork et al., 2022].
Finally, a class of methods speeds up the computation of Dual Weighted Residual [Becker and Rannacher, 1996, Becker and Rannacher, 2001] error estimators by learning a metric tensor from solution information [Fidkowski and Chen, 2021, Chen and Fidkowski, 2021] that can then be used in existing refinement procedures [Yano and Darmofal, 2012].

These supervised methods usually predict intermediary metrics for specifics facets of amr, typically proposing a greedy next-step refinement instead of focusing on long-term optimization.
We instead use rl to directly find optimal refinement strategies, which in our case consists of the error reduction for a given refinement sequence minus the cost of creating new mesh elements.

Reinforcement Learning for Adaptive Mesh Refinement.
Multiple strategies for rlamr have been developed, all of which are motivated by the potential of rl to optimize non-differentiable, long-term rewards.
Among these, one technique focuses on the generation of quadrilateral meshes through a process of iterative element extraction from the given problem domain [Pan et al., 2023].
Another approach utilizes a global threshold prediction for heuristic-based refinement, leveraging existing error estimates on the entire mesh [Gillette et al., 2022].
A similar method learns to optimize the polynomial order of mesh elements by evaluating the difference in solutions between the current polynomial degree and its immediate higher degree [Huergo et al., 2024].
These methods either learn a global, per-mesh quantity, or fully act on local mesh elements.
We instead directly manipulate the mesh elements themselves, iteratively considering for all mesh elements simultaneously whether they should be sub-divided.

When considering mesh refinement via element sub-division, a defining feature is the change in mesh topology through iterative refinements.
As the mesh grows with each refinement, the rl policy needs to be able to handle varying state, and potentially observation and action spaces.
Some work circumvents this issue by only considering a single refinement level each for the element size and its polynomial order [Dzanic et al., 2024].
Another [Wu et al., 2023] jointly learns a gns and an edge-based rlamr strategy on the simulated mesh to automatically adapt the mesh resolution during the learned simulation.
The method deals with the potentially complex action space by sampling a fixed number of edges to coarsen and refine in each step.
Instead of optimizing the mesh with respect to any concrete quality metric, they reward the reduction in downstream simulation error of the learned gns compared to not using any refinement.

More similar to our method, a class of approaches considers iterative element sub-division with the goal of optimizing the trade-off between speed and accuracy for a classical simulator [Yang et al., 2023, Yang et al., 2023, Foucart et al., 2023].
These techniques are predominantly tailored for non-stationary pde, thus incorporating mesh coarsening operations essential for addressing time-variant problems.
While our approach can readily be adapted to time-dependent problems and mesh coarsening, we instead focus on static meshes and mostly stationary problems due to their prevalence in engineering [Nagarajan and Soghrati, 2018, Ho-Le, 1988].
In the context of these stationary problems, the challenge is to find multiple levels of precise refinement on the full mesh instead of either shallow or local refinements and coarsening operations.

Of these approaches, Single Agent [Yang et al., 2023] feeds an encoding of the full mesh into a gnn-based policy to provide a categorical action that selects a single element to refine.
Single Agent refines a single element at each time step, requiring to solve the resulting system of equations after each refinement.
Thus, the method suffers from inference runtime and has to learn from comparatively little information per environment sample.
Partially motivated by this, Sweep [Foucart et al., 2023] decides on a refinement for all mesh elements in parallel during inference.
To achieve this, the method is trained by iteratively selecting a random element of the mesh and determining its refinement based on an mlp that processes local and global features.
Yet, this random element selection during training leads to high variance in the state transitions and thus a more complex learning problem, especially since the prediction mostly considers local element features.
Additionally, the increased inference speed comes as the cost of a misalignment in the environment transition between training and inference [Yang et al., 2023].

We circumvent this issue by refining all mesh elements in parallel during both training and inference, viewing the mesh as a collaborative multi-agent system, with one agent per element.
Crucially, this paradigm implies that the number of agents changes and that individual agents may vanish between refinement steps, making it unclear how to assign credit for the individual agents' actions over time [Cohen et al., 2021].
Here, vdgn [Yang et al., 2023] sets a maximum refinement depth to upper bound the number of agents, and employs a value decomposition network [Sunehag et al., 2017] to learn the credit assignment by decomposing a shared Q-function.
While this method theoretically allows for efficient training and inference, this learned value decomposition quickly becomes challenging for larger meshes.
As the optimization of the individual agents depends on the quality of this decomposition, training quickly becomes unstable or diverges entirely.
We instead formulate amr as a Swarm rl problem to naturally integrate changing observation and action spaces.
We additionally utilize the spatial nature of amr to define an agent-wise reward, providing a strong and stable feedback signal to all agents during training.

§ ADAPTIVE SWARM MESH REFINEMENT++

method treats elements of a mesh as a swarm of homogeneous agents that collaborate to find an optimal refinement.
For this, each agent's state and observation are defined through its position in the mesh and includes local features of both the mesh and the pde.
At each step, all agents decide whether to mark their respective element for refinement.
Agents that refine their element receive a local reward based on how much their refinement has improved the quality of the underlying simulation.
Crucially, each refinement subdivides the refined elements, introducing new agents in the process.
We introduce a asmdp that features a mapping of agents over time, allowing us to optimize over multiple refinement steps by propagating reward from agents of later time steps back to related earlier agents.
Figure <ref> provides a schematic overview of method, and the following sections describe the individual aspects in more detail.

§.§ Adaptive Swarm Markov Decision Process.

We view mesh refinement as a collaborative multi-agent rl problem with a changing number of homogeneous agents and agent-wise rewards.
For this, we adopt a swarm rl view, adapting the SwarMDP framework [Šošić et al., 2017, Hüttenrauch et al., 2019] to vector-valued rewards and changing state, action and observation spaces.
Formally, let an asmdp be a tuple $\langle \mathbb{S}, \mathbb{O}, \mathbb{A}, P, \mathbf{r}, \xi, \phi \rangle$.
The state space, observation space and action space are given by $\mathbb{S}$, $\mathbb{O}$, and $\mathbb{A}$ respectively.
Since the number of agents changes over time, we denote subsets of the state, observation, and action spaces with $N$ agents as $\mathcal{S}^N\subset\mathbb{S}$, $\mathcal{O}^N\subset\mathbb{O}$, $\mathcal{A}^N\subset\mathbb{A}$.
The transition function $P : \mathcal{S}^N \times \mathcal{A}^N \rightarrow \mathcal{S}^M$ takes an action over $N$ agents and leads to a new state with $M$ agents.
The reward function $\mathbf{r} : \mathcal{S}^N \times \mathcal{A}^N \rightarrow \mathbb{R}^N$ takes the full state into account to produce a scalar reward for each agent.
Similarly, the observation function $\xi: \mathcal{S}^N \rightarrow \mathcal{O}^N$ calculates local observations for each agent from the global state.

We consider a finite-horizon setting with $T$ time steps.
In the scalar reward case, i.e., with $r(\mathbf{s},\mathbf{a})\in\mathbb{R}$, an optimal rl policy $\pi: \mathbb{O} \times \mathbb{A} \rightarrow [ 0, 1 ]$ generally maximizes the return, i.e., the expected cumulative future reward

\begin{equation}
\label{eq:scalar_return}
J^t:=\mathbb{E}_{\pi(\mathbf{a}|\xi(\mathbf{s}))}\left[\sum_{k=0}^{T} r(\mathbf{s}^{t+k},\mathbf{a}^{t+k})\right]\text{.}
\end{equation}

A key challenge of using rl for amr is the changing number of agents introduced by splitting agents that are marked to be refined.
To adapt the Markov Decision Process to changing numbers of agents within a single episode, we thus introduce an agent mapping $\phi^t \in \mathbb{R}^{N\times M}$.
Intuitively, each entry $\phi^t_{ij}$ describes the responsibility of agent $i$ at step $t$ for agent $j$ at step $t+1$.
The responsibilities of agents at step $t$ for agents at step $t'>t$ can then be computed via the matrix multiplication

Given this mapping, we can propagate the rewards of all future agents for which agent $i$ is responsible for at step $t$ to calculate a return

\begin{equation}
\label{eq:mapped_return}
\end{equation}

For training, e.g., a value function $V(\mathbf{s})$, this leads to a TD error [Sutton and Barto, 2018]

\begin{equation}
\label{eq:td_error}
\delta^t= \mathbf{r}(\mathbf{s}^t, \mathbf{a}^t)_i+\sum_j \phi^t_{ij} V_j(\mathbf{s}^{t+1}) - V_i(\mathbf{s}^t)
\end{equation}

of actions .
We derive the targets for $Q$-functions analogously.

The asmdp framework allows for a collaborative optimization of changing numbers of agents, each of which is equipped with its own reward.
Compared to, e.g., a centralized partially observable mdp [Yang et al., 2023, Foucart et al., 2023], it provides a more efficient learning environment, since every environment step provides a learning signal for each agent.
Similarly, the asmdp supports a multi-agent setting with changing numbers of homogeneous agents, circumventing the posthumous credit assignment problem [Cohen et al., 2021] without having to resort to dummy states and learned value decompositions [Yang et al., 2023].
The agent mapping additionally allows for an optimization of the return, i.e., the reward over time, which would not be possible with a purely local per-agent view.
While this work focuses on amr, the asmdp framework can readily be applied to other scenarios that are characterized by localized decisions of changing numbers of entities.

§.§ Actions and Agent Mappings.

In the context of amr, a system state $s\in\mathbb{S}$ is defined by the problem domain $\Omega$, a mesh discretization $\Omega$, and a pde formulated on $\Omega$ to be solved under specific conditions.
method views every element $\Omega^t_i$ of a mesh $\Omega^t:=\{\Omega^t_i \subseteq \Omega |~ \dot{\bigcup}_i~\Omega^t_i = \Omega \}$ as one of many homogeneous agents.
Each agent has a binary action space, deciding on whether it wants to mark its element for refinement or not.
These markings are provided to a remesher, yielding a finer mesh $\Omega^{t+1} = \{\Omega_{j}^{t+1}\}_j$.
The remesher subdivides each refined element $\Omega_i^t$ into a set of smaller elements $\{\Omega_{j}^{t+1}\}_j$, such that the disjoint union of all $\Omega_{j}^{t+1}$ reconstructs $\Omega_i^t$, i.e., $\dot{\bigcup}_j \Omega_{j}^{t+1} = \Omega_i^t$.
It may also refine unmarked elements to close the mesh and thus assert a conforming solution [Arnold et al., 2000], i.e., to make sure that elements of the mesh align with each other at the boundaries to ensure continuity of solution variables between adjacent elements.

To map between elements and their successors, asmr uses the indicator function $\phi^t_{ij} := \mathbb{I}(\Omega^{t+1}_j\subseteq \Omega^t_i)$.
This mapping assigns each element to all elements that it subdivides into, or equivalently, lets each agent be responsible for all new agents that it spawns.
Here, we extend the mapping to include a normalization factor

\begin{equation}
\label{eq:agent_mapping}
\phi^t_{ij} := \frac{|\Omega^t|}{|\Omega^{t+1}|}\mathbb{I}(\Omega^{t+1}_j\subseteq \Omega^t_i)\text{.}
\end{equation}

Compared to asmr, the normalization factor $\frac{|\Omega^t|}{|\Omega^{t+1}|}$ scales the responsibilities such that $\sum \phi^t=\sum \phi^{t'}$.
This procedure intuitively ensures that the total `mass' of agents is preserved over time.
We find that this factor regularizes the mapping of rewards within the agent swarm, akin to, e.g., batch normalization [Ioffe and Szegedy, 2015].

Conceptually, we may represent each mesh in a series of refinements as a layer in a hierarchical graph, with nodes corresponding to mesh elements.
These nodes connect to their predecessor in the prior layer and to all successors in the subsequent layer.
Applying an agent mapping equates to navigating through the respective level of hierarchy of this graph.
An example refinement procedure and its corresponding refinement graph for an initial mesh with a single element and $2$ refinement steps can be seen in Figure <ref>

This work focuses on adaptive mesh refinement via iterative mesh subdivision.
For mesh coarsening, the above mapping can be extended to respect element coarsening as, e.g.,
In this case, the hierarchical graph is no longer a tree but instead a directed acyclic graph, since multiple elements can coarsen into the same new element.

§.§ Observations and Policy Architecture.

Given a mesh $\Omega^t$ and its corresponding system of equations, we construct an observation as the bidirectional graph $\gG_{\Omega^t} = \gG = (\mathcal{V}, \mathcal{E}, \mathbf{X}_\mathcal{V}, \mathbf{X}_\mathcal{E})\in\mathcal{O}$ with mesh elements as nodes $\mathcal{V}$ and their neighborhood relation as edges $\mathcal{E} \subseteq \mathcal{V}\times \mathcal{V}$.
Mesh- and task-dependent node and edge features of dimensions $d_\mathcal{V}$ and $d_\mathcal{E}$ are provided as $\mathbf{X}_\mathcal{V}: \mathcal{V}\rightarrow \mathbb{R}^{d_\mathcal{V}}$ and $ \mathbf{X}_\mathcal{E}: \mathcal{E}\rightarrow \mathbb{R}^{d_\mathcal{E}}$.

This observation graph is consumed by an mpn policy.
mpn [Pfaff et al., 2021, Linkerhägner et al., 2023] are a popular gnn architecture for mesh-based physical simulation [Pfaff et al., 2021, Linkerhägner et al., 2023, Würth et al., 2024] as they encompass the function class of several classical pde solvers [Brandstetter et al., 2022].
Given an input graph, they iteratively update latent features on this graph for $L$ consecutive message passing steps.
Each step $l$ receives the output of the previous step and updates the features $\mathbf{X}_\mathcal{V}$ , $\mathbf{X}_\mathcal{E}$ for all nodes $v\in \mathcal{V}$ and edges $e\in\mathcal{E}$.
Using linear embeddings $\mathbf{x}_v^0$ and $\mathbf{x}_e^0$ of the initial node and edge features, the $l$-th step is given as

\mathbf{x}^{l+1}_{e} = f^{l}_{\mathcal{E}}(\mathbf{x}^{l}_v, \mathbf{x}^{l}_u, \mathbf{x}^{l}_{e}), \textrm{ with } e = (u, v)\text{,}

\mathbf{x}^{l+1}_{v} = f^{l}_{\mathcal{V}}(\mathbf{x}^{l}_{v}, \bigoplus_{e=(v,u)\in \mathcal{E}} \mathbf{x}^{l+1}_{e})\text{.}

The operator $\oplus$ is a permutation-invariant aggregation such as a sum, mean or maximum over all aggregated elements.
Each $f^l_\cdot$ is a learned function, usually parameterized as a simple mlp.
The mpn's final output is a learned representation $\mathbf{x}^L_v$ for each node $v\in\mathcal{V}$.
This learned representation can be fed into a policy head to compute an action per agent, yielding a joint action vector $\pi(\mathbf{a}|\mathcal{G}_{\Omega^t})$, or into a value head to compute a per-element value function estimate.
Both are parameterized as mlp that are shared between the different agents.
Since mpn are permutation-equivariant by design [Bronstein et al., 2021], this policy architecture ensures that the amount and ordering of mesh elements does not matter, i.e., that the elements are only defined by their features and local position in the mesh.

§.§ Reward Definition

amr aims to generate meshes with enhanced local resolution in areas that are relevant for the underlying system of equations.
Its objective is to achieve an optimal balance between the solution's accuracy on the mesh $\Omega^t$, and its total element count $\Omega^t_i \in \Omega^t$.

Since closed-form solutions are generally not available for most systems of equations, we compute an approximate ground truth solution $u_{\Omega^*}$ using a fine-grained uniform reference mesh $\Omega^*$ [Yang et al., 2023].
Calculating an accurate ground truth solution on a uniform mesh is expensive, because the mesh requires a high resolution and thus a large number of elements.
Given $u_{\Omega^*}$, a refined mesh $\Omega^t$ and its solution $u_{\Omega^t}$, we define the error per element as the maximum difference between the solutions on this element.
We use the maximum difference over solutions as a measure of simulation quality that is independent of, e.g., the size of the used elements.
Given reference elements $\Omega^*_m$ with midpoints $p_{\Omega_m^*} \in \Omega^*_m$ as sampling points to query the solutions on, this yields an error

\begin{equation}
\label{eq:err_per_element}
\hat{\text{err}}(\Omega_i^t)\approx \max_{\Omega_m^*\subseteq \Omega_i^t} \Big|u_{\Omega^*}(p_{\Omega_m^*})-u_{\Omega^t}(p_{\Omega_m^*})\Big|\text{.}
\end{equation}

We efficiently calculate the assignment $\Omega^*_m\in\Omega^t_i$ required for this reward using a $k$-d tree [Bentley, 1975].
If a queried point $p_{\Omega_m^*}$ exactly lies on the edge between two elements in $\Omega^t$, which is possible for some refinement algorithms, we assign it to both of these elements with a weight of $0.5$.
We normalize this error estimate with the initial mesh error to get the relative improvement of the mesh-error, i.e.,
This estimate is consistent across geometries and process conditions.
As the size of the solution $u_{\Omega^*}$ and subsequently that of the error $\hat{\text{err}}(\Omega_i^t)$ may differ significantly between systems of equations due to varied process conditions, this normalization prevents systems of equations with large solution quantities from dominating the learning process.
Given the normalized error estimate, we next construct an element-wise local reward function that estimates the benefit of a refinement and compares it to the cost of adding additional elements.
We introduce a variable element penalty $\alpha$ that penalizes each added element, and compare this to the decrease in error within each element when it is refined, yielding the reward

\begin{equation}
\label{eq:asmr_local_reward}
\mathbf{r}(\Omega^t_i) := \left(\text{err}(\Omega^t_i)
-\alpha\left(\sum_j\mathbb{I}(\Omega^{t+1}_j\subseteq \Omega^t_i)-1\right)\text{.}
\end{equation}

The local reward for each element is $0$ if no refinement is made.
If an element is refined, its reward is positive if and only if the refinement decreases the highest error in this element by more than the cost of adding new elements, which is scaled with the element penalty $\alpha$.
In comparison, asmr [Freymuth et al., 2023] rewards the integrated average change in error for each mesh element instead of its maximum change in error.
This average change in error is integrated over the element volume, meaning that small elements, i.e., those that have been previously refined, quickly become comparatively insignificant for the reward.
To counteract this, asmr scales each reward by the inverse of the respective element's volume, leading to widely varying and potentially unbound reward values that do not directly match any optimization criteria on the full mesh.
We instead directly reward the change in maximum error within the volume of each element.
This maximum reward variant has been proposed as an ablation study in Freymuth et al., 2023, but only been explored for a single task.
We find that the maximum reward results in more well-behaved optimization, presumably because it directly targets refining elements where the reduction in error justifies adding new elements.
We provide a more detailed comparison to the asmr reward in Appendix <ref>.
Since Equation <ref> is based on an element's decrease in maximum error, it is independent of the element size, ensuring that rewards of elements of different sizes are on the same scale without the need for explicit scaling.

In systems of equations with multiple quantities of interest, such as, e.g., elasticity problems where both the element displacement and stresses are important, the mesh must be optimized for all the quantities of interest.
Here, we calculate rewards independently for each solution dimension and use some application-dependent convex function to compute a scalar reward per element.
Equation <ref> rewards each agent for its local refinement decision for each time step.
While this direct and explicit local credit assignment leads to a temporally and spatially dense reward, it disregards non-local effects of mesh refinement, which are common for elliptical pde.
To encourage global optimization, we therefore extend the return $J_i^t$ of Equation <ref> to include a global term as

\begin{equation}
\label{eq:half_half_return}
{J_i^t}'=\frac{1}{2} J^t_i+\frac{1}{2} J^t\text{,}
\end{equation}

where $J^t$ is the global return of Equation <ref> using the average reward $r=\frac{1}{N}\sum_j \mathbf{r}_j$.

§.§ Element Penalty

We add the element penalty $\alpha$ in Equation <ref> an input to the policy to inform it about the penalty of adding new elements.
During training, we draw $\alpha$ from a log-uniform distribution at the start of each episode, allowing for environment samples on a wide range of different penalties.
Similar to vdgn [Yang et al., 2023], we then specify a value of $\alpha$ during inference to control the refinement level of the created mesh.
Compared to asmr, which trains a policy on a fixed element penalty value, we thus condition the policy on $\alpha$ for each rollout.
This conditional policy allows for different mesh granularities during inference, whereas asmr needs to train a new policy when the desired number of mesh elements changes.

Balancing the extra computational cost of adding mesh elements against the benefit of reducing simulation errors can be seen as an instance of multi-objective rl [Van Moffaert and Nowé, 2014, Hayes et al., 2022].
The cost of adding new elements increases linearly with a penalty factor $\alpha$, and there exists some approximate ranking of the effectiveness of different mesh refinements at reducing simulation errors.
This setup creates a roughly convex relationship between $\alpha$ and the policy's behavior.
In other words, lowering $\alpha$ generally leads to a finer mesh that includes all refinements made with a higher $\alpha$, while raising $\alpha$ results in fewer refinements throughout the mesh.
Figure <ref> provides an example of resulting method meshes for the same trained policy and Neumann boundary task for different values of the element penalty $\alpha$.
Regardless of the penalty, the parts of the mesh with the highest potential error reduction, in this case the areas near the load function and the sinusoidal Neumann boundaries, are refined the most.

§ EXPERIMENTS

§.§ Systems of Equations.

The physical behavior of a technical system can be described by a set of pde.
Simulations predict the physical states of a system by solving those pde.
For non-trivial problems, the solving method is usually based on numerical approximation methods, such as the fem, which considers the weak formulation of the pde.
Multiplying the strong pde formulation with a test function $v(\textbf{x}) \in V$ and integrating by parts yields the weak formulation [Larson and Bengzon, 2013], which we denote as $\Phi(u,v)$.
We denote $V$ as the space spanned by a set of basis functions $\phi_v(\textbf{x}) \in V$ and the sought solution as $u(\textbf{x})$.
In the Galerkin fem, the solution $u(\textbf{x}, t)$ is approximated by a weighted sum of these basis functions $u(\textbf{x}) = \sum_{v=1}^{N_\varphi} u_v(t) \phi_v(\textbf{x})$.
Additionally, the same basis function are used to define the test function $v(\textbf{x}) = \sum_{v=1}^{N_\varphi} \phi_v(\textbf{x})$.
In this work we consider linear elliptical or parabolic pde.
For elliptical pde, we seek for a space-dependent solution $u=u(\textbf{x})$.
This results in $N_\varphi$ equations for $N_\varphi$ unknown degrees of freedom $u_v$.
Due to the linearity of the pde and the discretization of the fem, we can assemble a linear system of equations $\textbf{A}\textbf{u}=\textbf{b}$ out of the weak formulation $\Phi(\phi_u, u_v, \phi_v)$.
The time-dependent parabolic pde requires a discretization in time, e.g., using forward or backward Euler schemes.
Subsequently, a system of equation can be assembled for each time step.
We can efficiently solve the obtained linear system of equations with numerical solvers.
For each pde, we create a task that consists of the weak formulations, the domain and varying process conditions, and aim to generate meshes that minimize the error of the fem when compared to a fine-grained reference mesh while using few mesh elements.

Figure <ref> shows examples for each task and visualizes final meshes created using an method policy.
Each task is associated with a family of domains, including L-shapes, rectangles with a square hole or multiple rhomboid holes, convex polygons, star shapes with different numbers of peaks, and a $3$-dimensional plate.
From left to right and top to bottom, these tasks are briefly characterized as follows.
Laplace's equation task considers a rectangular heat source at the inner boundary and features high solution gradients near this source.
Poisson's equation uses a multi-modal load function comprised of a Gaussian mixture model, which in turn causes multiple distinct regions of interest in the solution.
The Stokes flow task solves for the velocity of a fluid streaming into the domain from inlet on the left side.
It uses more complex shape functions and requires high precision near the corners of the rhomboid obstacles.
The linear elasticity task simulates the deformation and the resulting stress of a metal plate.
The heat diffusion task simulates a time-dependent heat source and is evaluated at the final simulation step.
Adding sinusoidal Neumann Boundaries to Poisson's Equation models scenarios where the boundary flux is specified, resulting in high solution variability near the domain edges.
The $3$d Poisson task considers a $3$-dimensional variant of Poisson's equation, using tetrahedral meshes and a different underlying refinement algorithm.

We implement the fem, all systems of equations, meshes and refinements in scikit-fem [Gustafsson and Mcbain, 2020].
All tasks use conforming triangular meshes and linear elements unless mentioned otherwise and are provided as OpenAI gym [Brockman et al., 2016] environments.
We specify the systems of equations and their process conditions in Appendix <ref>.

§.§ Adaptive Mesh Refinement.

For each system of equations, we generate an initial coarse mesh $\Omega^0$, which is then refined during an episode through iterative element subdivision.
In $2$-dimensional domains, we employ the red-green-blue refinement method [Carstensen, 2004], subdividing each marked element into $4$ smaller elements and then closing the mesh, as illustrated in Figure <ref>.
For $3$-dimensional meshes, we apply longest edge bisection [Rivara, 1984, Suárez et al., 2005] to halve each marked element.
In $2$ dimensions, we uniformly refine the initial coarse mesh $6$ times to create the reference mesh $\Omega^*$, and use the same amount of $6$ refinement steps per episode.
We additionally evaluate a simpler task setup with $4$ refinements for the reference and amr algorithm to emulate the task complexity of existing work [Yang et al., 2023, Fortunato et al., 2022, Yang et al., 2023].
In both cases, always refining all elements precisely reconstructs the reference mesh.
In $3$-dimensional domains, we uniformly refine three times, resulting in $8^3$ times more elements in $\Omega^*$ and employ $7$ longest edge bisection amr steps, yielding refinements that are always topologically different from the reference mesh.

§.§ Training and Evaluation.

We train all rlamr methods on $100$ randomly generated systems of equations, including randomized domains and process conditions.
This choice reduces the number of expensive reference meshes $\Omega^*$ that must be computed during training, and is further explored in Appendix <ref>.
We evaluate each converged rlamr policy and the heuristics on $100$ randomly sampled but fixed evaluation pde, training each rlamr algorithm for $10$ random seeds.
The systems of equations used for evaluation are disjoint from those used for training.

For all methods, we experiment with a range of target mesh resolutions to produce meshes of different refinement levels, as detailed in Appendix <ref>.
Appendix <ref> provides additional details on the experimental setup and additionally discusses the required computational budget.
We train the rlamr methods using both ppo [Schulman et al., 2017] with discrete actions as an on-policy algorithm and dqn [Mnih et al., 2013, Mnih et al., 2015] as an off-policy variant.

We evaluate mesh quality by calculating the squared error between the solution on this mesh and the solution on the fine-grained reference $\Omega^*$.
This squared error over the domain is numerically approximated over evaluations at the midpoints $p_{\Omega_m^*}$ of the reference mesh as

\begin{equation}
\label{eq:squared_integrated_difference}
\sum_{\Omega_m^*\in \Omega^*} \text{Volume}(\Omega^*_m)\left(u_{\Omega^*}(p_{\Omega_m^*})-u_{\Omega^t}(p_{\Omega_m^*})\right)^2\text{.}
\end{equation}

We normalize the resulting value by that of the initial mesh for comparability across pde, and call the resulting metric the remaining squared mesh error.
This squared error metric captures the overall error of the solution on the refined mesh, while punishing outliers.
It differs from our optimization objective in Equation <ref>, which instead rewards minimizing the maximum error of the solution on the refined mesh.
Appendix <ref> discusses and evaluates a mean error and an approximation of the maximum error as alternate metrics.
We list further algorithm and network hyperparameters in Appendix <ref>.

§.§ RL-AMR Baselines.

We compare to asmr [Freymuth et al., 2023] and several recent non-stationary rlamr methods [Yang et al., 2023, Yang et al., 2023, Foucart et al., 2023] that have been adapted to our application, focusing on stationary refinements.
Following asmr, we use the rewards as described in the respective papers for each method, employing the absolute integrated difference in solution compared to the reference mesh, i.e., Equation <ref> without the square as an error estimate.

Single Agent [Yang et al., 2023] uses the difference in error norm as a global reward and marks a single element for refinement in each step via a categorical action over the full mesh.
Sweep [Foucart et al., 2023] randomly samples a single mesh element and decides on its refinement based on local features and a global resource budget during training.
The method uses the logarithm of the change in error between consecutive steps as a reward.
During inference, each timestep considers all mesh elements in parallel.
vdgn [Yang et al., 2023] employs value decomposition networks [Sunehag et al., 2017] to estimate a global Q-function as the sum of element-wise local Q-functions.
For a vdgn-like ppo variant, we instead analogously decompose the value function.
As our method is an extension of asmr [Freymuth et al., 2023], this baselines largely follows our method. However, it instead uses the volume-scaled reward in Equation <ref>, omits the normalization factor in the agent mapping in Equation <ref>, considers an infinite horizon rl setting, and uses less regularization in the neural network architecture.

We use an mpn policy for all rlamr methods except for Sweep, which utilizes a simple mlp.
The rlamr methods handle the trade-off between mesh resolution and solution accuracy slightly differently.
Single Agent simply refines for a number of steps, Sweep uses a target number of elements, and the other methods penalize each added element.
Here, asmr trains a new policy for each penalty, while only method and the vdgn-like baseline use a conditional policy that can produce meshes of different resolutions during inference.
Appendix <ref> lists all baseline-specific hyperparameters, while Appendix <ref> provides details on the mesh resolution parameters for all methods.

§.§ Heuristic Baselines.

In addition to the rlamr methods above, we compare to a traditional error-based heuristics acting on a refinement threshold $\theta$ [Binev et al., 2004, Bangerth et al., 2012, Foucart et al., 2023].
Given an error estimate $\text{err}(\Omega^t_i)$ for element $i$, the heuristic marks all elements for refinement for which
$\text{err}(\Omega^t_i)>\theta\cdot\max_j \text{err}(\Omega^t_j)$.
We consider an Oracle Error Heuristic that uses the absolute integrated difference in solution as its error metric, as well as a Maximum Oracle Error Heuristic, for which the error estimate of Equation <ref> is used.
Both variants require the reference mesh $\Omega^*$ to estimate an error, which is expensive to compute and thus may be unavailable during inference.
We additionally consider the commonly used zz, which instead uses a superconvergent patch recovery process for its error estimate [Zienkiewicz and Zhu, 1992] and does not require access to $\Omega^*$.
As the recovery process averages over neighboring mesh elements, the zz generally produces smooth error estimates.
This property can lead to more coherent refinements when compared to the purely greedy oracle heuristics, especially when considering the closing operations needed to ensure a conforming mesh after every refinement step.

All heuristics greedily refine based on local error estimates, implicitly assuming that refining elements with high error is the best way to minimize the error on the mesh.
This assumption does not always hold, as, e.g., elliptical pde often feature globally propagating errors [Strauss, 2007, Foucart et al., 2023].
The rlamr methods instead learn to optimize an objective function, which rewards the error reduction over time.
This difference between greedily marking elements with a high local error, and learning to anticipate the long-term error reduction for any given refinement allows the rlamr methods to find more global strategies that take multiple refinement steps and non-local information of the mesh elements into account.

§.§ Observation Graph.

As node features for the observation graph of Section <ref>, we use the mean and standard deviation of the fem solution on the element's vertices, the element volume and the current environment timestep.
For method and the vdgn-like baseline, we also add the current element penality to condition the refinement process of the policy on a given target resolution.
Some tasks employ additional node features to encode their process conditions, which we describe in Appendix <ref>.
We use the Euclidean distance between the element midpoints as the single edge feature.
Compared to using either a distance vector between element midpoints, or providing absolute positions, this ensures that the predicted sizing fields are invariant under Euclidean group operations [Pfaff et al., 2021, Bronstein et al., 2021].

§ RESULTS

§.§ Qualitative Results.

Figure <ref> visualizes a full rollout of method on an instance of the heat diffusion task, including the per-element refinement markings after every policy step.
The method is trained to produce a sequence of refinements that improve the solution accuracy enough to justify the additional mesh elements, which here corresponds to more refinements along the path of the heat source.
Figure <ref> shows final refined meshes for the same Poisson problem for a fixed method policy that is conditioned on different values of the element penalty $\alpha$.
Here, decreasing the element penalty results in more mesh elements and thus a more accurate solution.
Our method is able to focus on the interesting parts of the mesh regardless of the final mesh resolution, thus providing a favorable trade-off between solution accuracy and computational cost.
Figure <ref> provides final refinements of method on exemplary systems of equations for all considered tasks.
The refinement strategy individually adapts to each task, choosing to refine elements that decrease the simulation error regardless of the underlying system of equations.
We further explore this error reduction property in Figure <ref>.
The visualization shows that our method achieves a much more uniform distribution of errors across the mesh when compared to simple uniform refinement, while also acquiring a much smaller simulation error for the same number of mesh elements.

Appendix <ref> visualizes all methods on the Stokes flow task to showcase the different refinement behaviors  for different target mesh resolution parameters.
Appendix <ref> provides additional method visualizations for the remaining tasks, including different camera angles for the $3$-dimensional variant of Poisson's equation that uses tetrahedral meshes.

§.§ Quantitative Results.

Quantitatively, we measure the mesh quality using a Pareto plot of the number of elements and the normalized error of the mesh.
We evaluate each rlamr policy on $100$ evaluation environments per mesh resolution parameter.
We then calculate the interquartile mean [Agarwal et al., 2021] of the resulting $100$ used elements and normalized errors of the final produced meshes, and plot the resulting values.
This process is done for $10$ independently trained policy seeds and repeated across a range of mesh resolution parameters.
For methods that require a fixed target mesh resolution, we use $10$ different parameters, for the adaptive methods we instead use $20$ parameters, as detailed in Appendix <ref>.
We further provide a log-log quadratic regression over the aggregated results of each method as a general trend-line.
To improve clarity and focus on representative behavior,  we omit data points with abnormally high element counts, which are sometimes produced as outliers by the vdgn-like baseline.

We experiment with ppo as an on-policy algorithm and dqn as an off-policy variant of the rl backbone for all rlamr methods on the Stokes flow task in Appendix <ref>.
All rlamr methods, including method, perform better with ppo, indicating that an on-policy algorithm is favorable when dealing with action and observation spaces of varying size.
Further, we compare our mpn to the gat-like [Velickovic et al., 2018] architecture proposed by vdgn [Yang et al., 2023] in Appendix <ref>.
For both our method and the vdgn-like baseline, the mpn shows minor but consistent benefits in performance.
We thus use ppo and mpn in all other experiments.

Appendix <ref> compares the performance of the zz Heuristic when applied directly to the initial mesh to variants that instead uniformly refine the mesh once or twice before using the estimator.
We find that the method benefits from $2$ uniform initial refinements, and use this variant for all experiments.
This tuning of the initial meshes in not required by the rlamr methods and prevents coarse refinements.

Using these initial results, we compare the different approaches on Laplace's equation.
On the left of Figure <ref>, we evaluate a simplified setup that uses $4$ refinement steps for both the reference mesh and the amr algorithms, except for Single Agent, which instead refines for up to $100$ time steps.
All methods clearly outperform uniform refinements.
The zz Heuristic performs well for larger meshes, but can not produce meshes with few elements as it starts with $2$ uniform initial refinements.

The right of Figure <ref> increases the complexity to $6$ refinement steps, and $400$ time steps for Single Agent.
Similarly, Figure <ref> provides results on the remaining $6$ tasks.
Across all tasks, only asmr and method can handle larger meshes while the other rlamr methods fail to provide better-than-uniform meshes or exhibit high variance in the quality of the proposed refinements.
Our method additionally outperforms all heuristics, including the Oracle Error and Maximum Oracle Error Heuristics, in several instances.
These results indicate that our optimization method facilitates learning non-local long-horizon refinement strategies that minimize the mesh error instead of simple targeting elements that currently have a high error.
method further improves over asmr in most settings, especially in the complex Stokes flow task where it compares well to the different Heuristics.
In the tasks where method does not perform better than asmr, both perform similar to or better than all heuristics, potentially indicating that both methods offer near-optimal refinements for these tasks.
On the $3$d Poisson's Equation, the different methods quickly converge to a relatively high final mesh error, likely due to the difference in mesh topology between the refined meshes and the uniform reference that is caused by the longest edge bisection refinement strategy.

Appendix <ref> presents additional results using a mean error and an adapted maximum error metric.
On both metrics, method is generally competitive with or better than the Heuristics and significantly outperforms the Single Agent, Sweep and vdgn-like baselines.
Comparing the different metrics, the Maximum Oracle Error Heuristic improves over the Oracle Error Heuristic on the maximum error metric and vice versa.
Likewise, asmr, which uses a scaled variant of the mean error instead of the maximum objective of Equation <ref>, outperforms method on half of the tasks in terms of the remaining mean error, while method is superior for every setup when evaluating the maximum error metric.

§.§ Ablations and Additional Experiments.

We conduct several ablations and additional experiments on the challenging Stokes flow task to determine what makes our approach uniquely effective.
These experiments individually change a single aspect of method while leaving the rest of the approach unchanged.

We compare the local maximum error reward of Equation <ref> to different local and global variants.
For the global variants, we obtain a scalar reward function by averaging over the local reward functions. Here, we compute the return using Equation <ref>, i.e., without mapping between agents over time.
In addition to the maximum reward of Equation <ref>, we compare to the volume-scaled reward proposed by asmr in Equation <ref>, and a variant that simple rewards a decrease in integrated error, i.e., Equation <ref> without the volume scaling term.

Figure <ref> shows the results.
Using any type of global reward function leads to unstable refinements, especially for large meshes, likely because the credit for each reward cannot be properly assigned to the large number of agents in the system.
For the local reward variants, the local maximum reward performs best, closely followed by the volume-scaled variant.
Simply integrating the error of each mesh element and rewarding a reduction in this integrated error performs poorly for finer meshes.
This result is likely the case because smaller elements generally contain less integration points of the reference mesh, and thus have less integrated error to reduce.
Optimizing this reward produces meshes that have an even distribution of integrated element errors, rather than refinements that effectively decrease the average error on the mesh.

Return and Agent Mapping.

Equation <ref> enables the computation of a TD error across agents at successive time steps by mapping each agent at time step $t$ to every agent it creates at time step $t+1$.
In this setup, each spawned agent is fully attributed to its predecessor, meaning the originating agent is considered wholly responsible for any agent it generates.
Compared to asmr, we additionally apply a regularization to this mapping through the ratio $\frac{\Omega^t}{\Omega^{t+1}}$.
An alternative approach would be to instead limit the total responsibility per agent to $1$, achieved by averaging the mapping as

\phi'^t_{ij} :=\frac{\mathbb{I}(\Omega^{t+1}_j\subseteq \Omega^t_i)}{\sum_{j'}\mathbb{I}(\Omega^{t+1}_{j'}\subseteq \Omega^t_i)}\text{.}

We assess these two methods in Figure <ref>, finding that the normalization factor $\frac{\Omega^t}{\Omega^{t+1}}$ significantly boosts performance, likely due to its role as a regularizer during training.
Furthermore, using a sum mapping, i.e., assigning each new agent a weight of $1$ from its creator, proves more effective than the mean mapping, which averages the total mapping weights per old agent.
We additionally explore optimizing a return that omits the global term in Equation <ref>, i.e., that only optimizes the local return of each agent.
Here, we find that the partially global objective improves performance, likely because it aids in global decision making.

Target Mesh Resolutions
Each rlamr method uses a single parameter to control the number of target elements of the final refined mesh.
method and vdgn condition their policy on an element penalty $\alpha$, while asmr specifies a fixed element penalty as a hyperparameter.
Similarly, Sweep uses a fixed budget $N_{\text{max}}$ per policy, and Single Agent trains on a specified number of rollout steps $T$.
Additional details are provided in Appendix <ref> and Table <ref>.

Figure <ref> compares training method on a fixed element penalty, as done in asmr [Freymuth et al., 2023], to training on an adaptive penalty that the policy is conditioned on for each rollout.
It additionally evaluates sampling the penalty uniformly instead of log-uniformly during training.
Interestingly, training method on a fixed element penalty does not improve performance, suggesting that the method policy is able to learn to condition on a concrete $\alpha$ value.
These results imply that method can successfully learn optimal tradeoffs between element cost and error reduction and apply it during inference.
Thus, we can train a single policy for a range of mesh granularities, omitting the need for asmr's expensive re-training for each target granularity.
If we train on uniformly instead of log-uniformly sampled values of $\alpha$, the method predominantly provides meshes with low numbers of elements, likely because there is a logarithmic relationship between the element penalty and the number of mesh elements.

Figure <ref> visualizes the methods' performance for different target resolutions.
Here, method provides meshes with consistent numbers of elements for a given target resolution, while the other rlamr baselines produce poorly refined meshes or inconsistent refinement granularities over target resolutions.

Additional Ablations
Appendix <ref> provides additional ablations over changes in the network architecture, varying numbers of training pde and different scales of the element penalty parameter in the reward.
The results indicate that as few as $10$ training pde are sufficient to provide accurate refinements, with performance increasing for up to $100$ pde.
In terms of architecture, edge dropout [Rong et al., 2019] of $0.1$, i.e., randomly omitting $10\,\%$ of edges in the mesh topology during training, slightly improves performance.
Linearly increasing or decreasing the range of the element penalty $\alpha$ leads to meshes with fewer or more elements, respectively, yet method always provides stable and efficient refinements for the respective element penalty ranges.

§.§ Runtime Comparison.

Figure <ref> shows the runtime comparison between our method and direct computation of the uniform reference $\Omega^*$ across all tasks.
Our method's timing includes the creation of the initial mesh, and iteratively solving the system of equations on the mesh, creating an observation graph, querying the policy for a refinement strategy, and finally refining a total of $6$ times.
For the uniform mesh, we simply measure the time it takes to refine the initial mesh $6$ times and to subsequently solve the problem on the resulting mesh.
We use a single 8-Core AMD Ryzen 7 3700X Processor for all measurements.

Our approach is always significantly faster than computing the fine-grained reference despite the comparatively large computational overhead.
A higher squared error requires less computation time, showing that method can adapt to task-specific computational budgets.
The Stokes flow task uses comparatively expensive $P_1/P_2$ Taylor-Hood-elements [John et al., 2016], causing our method to be more than $30$ times faster than $\Omega^*$ even for highly refined and accurate final meshes.

The other rlamr methods feature a similar iterative refinement procedure and thus runtime, but produce meshes of worse quality when compared to method.
In contrast, the Oracle Error Heuristic and Maximum Oracle Error Heuristic produce high-quality meshes, but require the calculation of $\Omega^*$, which dominates their runtime.
The zz Heuristic does not require access to $\Omega^*$ and is thus comparatively fast, but needs to tune the initial mesh and varies in solution quality depending on the considered task.

§.§ Generalization Capabilities.

We further experiment with the generalization capabilities of method to unseen and larger domains during inference.
Appendix <ref> shows the generalization capabilities of our approach when evaluated for different domains and load functions for Poisson's equation, indicating that it strongly generalizes to arbitrary domains and process conditions during inference.

We additionally train on a variant of Poisson's equation that may sample load functions modes outside of the considered domain and randomizes solution values on the domain's boundaries.
These changes effectively modify the training pde to mimic small patches of a larger mesh, acting as a data augmentation technique that facilitates up-scaling to larger domains and meshes during inference.
An exemplary refinement on a $20\times20$ from a policy trained on $1\times1$ domains can be seen in Figure <ref>.
The refined mesh has more than $50\,000$ elements, which is significantly larger than any meshes seen during training, yet accurately captures the interesting parts of the solution.
The full refinement procedure is roughly $100$ times faster than solving the fine-grained reference $\Omega^*$.

Figure <ref> shows the mesh error and inference speed for method policies trained on the augmented setup for Poisson's equation when evaluated on domains of different sizes.
For the $9\times9$ domain we used $10$ instead of $100$ evaluation pde due to their increased runtime and memory requirements.
While the number of elements for uniform refinements is linear in the domain's volume, meshes created by our method grow less quickly to achieve the same error threshold.
These results suggest that there are fewer elements with a significant error in larger domains.
Accordingly, evaluating the reference mesh $\Omega^*$ quickly becomes expensive as the domain size gets larger, while method provides fast and efficient refinements.
For large meshes, our method is about $100$ times faster than evaluating $\Omega^*$ while maintaining a mesh error of less than $0.001$.
Training details, as well as additional visualizations for adapted training pde and final refined meshes for different domain sizes are provided in Appendix <ref>.

The generalization capabilities of our method are likely a result of the agent-wise optimization and the utilized mpn, which facilitate refinement strategies that are largely based on local element neighborhoods.
Especially the ability of method to efficiently up-scale to larger meshes during inference opens up promising applications for practical engineering applications.
A policy can be trained on small and cheap environments, and then used in larger and more challenging setups during inference.

§ CONCLUSION

We introduce method, an adaptive mesh refinement method leveraging swarm reinforcement learning that increases the efficiency of numerical simulations with the finite element method.
Our method iteratively evaluates a system of equations on a mesh, generates an observation graph from the mesh and solution, processes this graph with a message passing graph neural network, and subsequently marks mesh elements for subdivision.
We treat each mesh as a homogeneous agent in a collaborative multi-agent system, training all agents with a shared policy.
Crucially, we provide a per-agent reward that balances the local reduction in simulation error against the cost of adding additional mesh elements.
This cost can be tuned during inference to generate meshes of different resolution using the same policy.
Accommodating different mesh sizes within an episode due to subdivisions, we propose a location-based mapping that assigns each new mesh element to the element in the previous mesh that created it.

Experimental results on static, conforming $2$d and $3$d meshes demonstrate that this agent-wise approach creates high-quality refinements for meshes with thousands of elements and without requiring any expensive error estimate during inference.
method surpasses both existing reinforcement learning-based methods and traditional refinement techniques, achieving a mesh quality comparable to expensive oracle-based error heuristics.
We conduct additional experiments that highlight the benefits of coupling our per-agent reward function with our temporal agent mapping.
Further, the method directly generalizes to different process conditions and domains during inference, and can be combined with simulation-specific data augmentation in the training environment to scale to significantly larger domains and meshes during inference.
For the considered systems of equations, our method is up to $30$ times faster than using a uniform refinement in small domains and can exceed a speed-up of factor $100$ when scaling to larger evaluation setups during inference.

Limitations and Future Work
Our method requires solving the system of equations on the mesh following each refinement, which quickly becomes expensive for finer meshes.
To mitigate this, we aim to integrate auxiliary physics-based losses or policy distillation methods to implicitly learn to predict which regions to refine from just the problem's process conditions.
Additionally, our reward function is currently based on an expensive reference mesh and solution.
While this reference is only required during training and method generalizes well from as few as $10$ different training systems of equations, the maximum mesh resolution is limited by this reference.
Lastly, we currently focus on static problems.
In future work, we will extend our method to time-dependent systems of equations and thus include both mesh refinement and coarsening operations.

§ ACKNOWLEDGEMENTS AND DISCLOSURE OF FUNING

NF was supported by the BMBF project Davis (Datengetriebene Vernetzung für die ingenieurtechnische Simulation).
This work is also part of the DFG AI Resarch Unit 5339 regarding the combination of physics-based simulation with AI-based methodologies for the fast maturation of manufacturing processes. The financial support by German Research Foundation (DFG, Deutsche Forschungsgemeinschaft) is gratefully acknowledged. The authors acknowledge support by the state of Baden-Württemberg through bwHPC, as well as the HoreKa supercomputer funded by the Ministry of Science, Research and the Arts Baden-Württemberg and by the German Federal Ministry of Education and Research.

§ INFERENCE-TIME GENERALIZATION

§.§ Same-size Generalization Capabilities

Figure <ref> visualizes final method meshes on Poisson's equation on different $1\times1$ domains and load functions that are not seen during training.
We experiment on the $5$ different domain classes that are used for the different tasks, plus a rectangular $\Omega=(0,1)^2$ domain.
For each class, we randomly sample $3$ domains and apply a Gaussian mixture model load function with $1$, $3$, and $5$ components, respectively.
We find that our approach generalizes well to novel domains and load functions, likely due to the strong local inductive bias of the Swarm rl objective and the mpn policy.

§.§ Domain Size Up-Scaling

We additionally experiment with the abilities of asmr to up-scale to larger domains during inference.
This up-scaling is crucial for practical scenarios, as it allows applying policies that are trained on small and comparatively cheap domains to apply to complex, large-scale domains where training would otherwise be too expensive.

Using Poisson's equation as an example, we augment the training pde to cover a larger space of potential simulations and mimic larger mesh segments by altering boundary conditions and load functions.
We train our policies on the square hole domains of, e.g., Figure <ref>.
We sample the means of the load function from a centered Gaussian with a standard deviation of $0.2$, allowing components outside the mesh, and additionally apply random Gaussian loads to selected boundary edges to emulate a larger domain outside of the simulated mesh.
As these changes lead to training pde of varying complexity, we use $10\cdot \hat{\text{err}}(\Omega)$ instead of $\text{err}(\Omega)$ in Equation <ref>, i.e., we omit the per-domain normalization, and use an unlimited number of training pde to further increase variety.
These changes only affect the training environments, and do not infer with our training algorithm.
Figure <ref> visualizes examplary training pde and policy refinements.

Figure <ref> compares method trained on these augmented training pde with the setup used throughout the paper.
The results show that the augmented training setup leads to a slight decrease in performance relative to method trained under standard conditions.
This decrease in performance is likely caused by the more challenging optimization problem.
Yet, the augmented training yields high-quality refinements.

Once trained, we evaluate the policy on spiral-shaped domains of increasing sizes, ensuring the volume of the initial mesh elements remains constant.
To align the complexity of the load function with the enlarged mesh sizes, we adapt the number of Gaussian mixture model modes, using more modes for larger domains.
Specifically, we use $16$ components which are first placed on a uniform grid and then randomly perturbed, allowing for them to be positioned outside of the mesh.
Figure <ref> shows an exemplary method refinement for a $5\times5$ domain.
Figure <ref> shows a final refined mesh on a $20\times20$ domain with $100$ components, including detailed views of specific sections.
Both figures use the same policy with a fixed element penalty.

§ BASELINE VISUALIZATIONS

§.§ RL-AMR visualizations

We visualize the final meshes and corresponding pde solutions of all rlamr methods on the Stokes Flow task for $6$ different refinement granularities on the same randomly selected system of equations.
For method and vdgn, which condition the policy on the element penalty parameter, we use the same trained model for the $6$ mesh refinement granularities.
Figure <ref> visualizes the results.

We find that method provides accurate refinements for different values of $\alpha$ using a single policy.
The produced meshes are of higher quality than those created by asmr, especially near the corners of the rhomboid obstacles in the domain.
All other rl-amr methods are unable to produce consistent meshes.
Single Agent tends to over-refine uninteresting regions, likely due to its iterative refinement procedure.
In contrast, Sweep often collapses to uniform or mostly uniform refinements, which may be a by-product of the difference in training on individual elements and sweeping over all elements at once during inference.
vdgn produces high-quality refinements in some cases, but is inconsistent across seeds and in some cases fails to optimize for different element penalties with the same model, which could be because its value decomposition objective may scale poorly to large numbers of agents.
These qualitative findings are consistent with the error measurements of Figure <ref>.

§.§ Heuristic visualizations

We visualize exemplary refinements for the error estimation-based threshold Heuristics in Figure <ref>.
All Heuristics greedily refine the elements with their respective largest error estimates, regardless of the resulting decrease in error.
This behavior is fully local, which may cause issues for global dependencies [Strauss, 2007] and conforming refinements.
The zz Heuristic uses a smoother error estimate than the oracle heuristics, which leads to potentially sub-optimal but generally more consistent refinements.

§ ADDITIONAL ASMR++ VISUALIZATIONS

We visualize meshes created by method policies for all considered tasks.
For the $2$-dimension tasks, we follow Section <ref> and use a fixed policy to create refinements with $6$ different granularities for a random pde for each task, only varying the element penalty $\alpha$.
For the $3$-dimensional variant of Poisson's equation, we instead visualize the same random pde and $\alpha$ value for $6$ different camera angles.
Figure <ref> displays the final refined meshes, as well as the solutions of the pde on the meshes.
On all tasks, method provides highly accurate refinements across target mesh refinement levels.

[Agarwal et al., 2021]
Rishabh Agarwal, Max Schwarzer, Pablo Samuel Castro, Aaron C Courville, and Marc Bellemare.
Deep reinforcement learning at the edge of the statistical precipice.
Advances in neural information processing systems, 34:0 29304–29320, 2021.

[Allen et al., 2022]
Kelsey Allen, Tatiana Lopez-Guevara, Kimberly L Stachenfeld, Alvaro Sanchez Gonzalez, Peter Battaglia, Jessica B Hamrick, and Tobias Pfaff.
Inverse design for fluid-structure interactions using graph network simulators.
Advances in Neural Information Processing Systems, 35:0 13759–13774, 2022.

[Anderson et al., 2021]
Robert Anderson, Julian Andrej, Andrew Barker, Jamie Bramwell, Jean-Sylvain Camier, Jakub Cerveny, Veselin Dobrev, Yohann Dudouit, Aaron Fisher, Tzanio Kolev, et al.
Mfem: A modular finite element methods library.
Computers & Mathematics with Applications, 81:0 42–74, 2021.

[Andrychowicz et al., 2021]
Marcin Andrychowicz, Anton Raichuk, Piotr Stańczyk, Manu Orsini, Sertan Girgin, Raphaël Marinier, Leonard Hussenot, Matthieu Geist, Olivier Pietquin, Marcin Michalski, Sylvain Gelly, and Olivier Bachem.
What matters for on-policy deep actor-critic methods? a large-scale study.
In International Conference on Learning Representations, 2021.
URL <https://openreview.net/forum?id=nIAxjsniDzg>.

[Arnold et al., 2000]
Douglas N Arnold, Arup Mukherjee, and Luc Pouly.
Locally adapted tetrahedral meshes using bisection.
SIAM Journal on Scientific Computing, 220 (2):0 431–448, 2000.

[Ba et al., 2016]
Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton.
Layer normalization.
stat, 1050:0 21, 2016.

[Baker, 1997]
Timothy J Baker.
Mesh adaptation strategies for problems in fluid dynamics.
Finite Elements in Analysis and Design, 250 (3-4):0 243–273, 1997.

[Bangerth and Rannacher, 2013]
Wolfgang Bangerth and Rolf Rannacher.
Adaptive Finite Element Methods for Differential Equations.
Birkhäuser, 2013.

[Bangerth et al., 2012]
Wolfgang Bangerth, Carsten Burstedde, Timo Heister, and Martin Kronbichler.
Algorithms and data structures for massively parallel generic adaptive finite element codes.
ACM Transactions on Mathematical Software (TOMS), 380 (2):0 1–28, 2012.

[Battaglia et al., 2018]
Peter W Battaglia, Jessica B Hamrick, Victor Bapst, Alvaro Sanchez-Gonzalez, Vinicius Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, et al.
Relational inductive biases, deep learning, and graph networks.
arXiv preprint arXiv:1806.01261, 2018.

[Becker and Rannacher, 1996]
Roland Becker and Rolf Rannacher.
Weighted a posteriori error control in FE methods.
IWR, 1996.

[Becker and Rannacher, 2001]
Roland Becker and Rolf Rannacher.
An optimal control approach to a posteriori error estimation in finite element methods.
Acta numerica, 10:0 1–102, 2001.

[Bentley, 1975]
Jon Louis Bentley.
Multidimensional binary search trees used for associative searching.
Communications of the ACM, 180 (9):0 509–517, 1975.

[Berger and Colella, 1989]
Marsha J Berger and Phillip Colella.
Local adaptive mesh refinement for shock hydrodynamics.
Journal of computational Physics, 820 (1):0 64–84, 1989.

[Binev et al., 2004]
Peter Binev, Wolfgang Dahmen, and Ron DeVore.
Adaptive finite element methods with convergence rates.
Numerische Mathematik, 97:0 219–268, 2004.

[Bohn and Feischl, 2021]
Jan Bohn and Michael Feischl.
Recurrent neural networks as optimal mesh refinement strategies.
Computers & Mathematics with Applications, 97:0 61–76, 2021.

[Brandstetter et al., 2022]
Johannes Brandstetter, Daniel E Worrall, and Max Welling.
Message passing neural pde solvers.
In International Conference on Learning Representations, 2022.

[Brenner and Scott, 2008]
Susanne C Brenner and L Ridgway Scott.
The mathematical theory of finite element methods, volume 3.
Springer, 2008.

[Brockman et al., 2016]
Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba.
Openai gym.
arXiv preprint arXiv:1606.01540, 2016.

[Bronstein et al., 2021]
Michael M Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković.
Geometric deep learning: Grids, groups, graphs, geodesics, and gauges.
arXiv preprint arXiv:2104.13478, 2021.

[Bryan et al., 2014]
Greg L Bryan, Michael L Norman, Brian W O'Shea, Tom Abel, John H Wise, Matthew J Turk, Daniel R Reynolds, David C Collins, Peng Wang, Samuel W Skillman, et al.
Enzo: An adaptive mesh refinement code for astrophysics.
The Astrophysical Journal Supplement Series, 2110 (2):0 19, 2014.

[Carstensen, 2004]
Carsten Carstensen.
An adaptive mesh-refining algorithm allowing for an h 1 stable l 2 projection onto courant finite element spaces.
Constructive Approximation, 20:0 549–564, 2004.

[Cerveny et al., 2019]
Jakub Cerveny, Veselin Dobrev, and Tzanio Kolev.
Nonconforming mesh refinement for high-order finite elements.
SIAM Journal on Scientific Computing, 410 (4):0 C367–C392, 2019.

[Chen and Fidkowski, 2021]
Guodong Chen and Krzysztof J Fidkowski.
Output-based adaptive aerodynamic simulations using convolutional neural networks.
Computers & Fluids, 223:0 104947, 2021.

[Cohen et al., 2021]
Andrew Cohen, Ervin Teng, Vincent-Pierre Berges, Ruo-Ping Dong, Hunter Henry, Marwan Mattar, Alexander Zook, and Sujoy Ganguly.
On the use and misuse of absorbing states in multi-agent reinforcement learning.
arXiv preprint arXiv:2111.05992, 2021.

[Cunningham et al., 2009]
Andrew J Cunningham, Adam Frank, Peggy Varnière, Sorin Mitran, and Thomas W Jones.
Simulating magnetohydrodynamical flow with constrained transport and adaptive mesh refinement: algorithms and tests of the astrobear code.
The Astrophysical Journal Supplement Series, 1820 (2):0 519, 2009.

[Dzanic et al., 2024]
Tarik Dzanic, Ketan Mittal, Dohyun Kim, Jiachen Yang, Socratis Petrides, Brendan Keith, and Robert Anderson.
Dynamo: Multi-agent reinforcement learning for dynamic anticipatory mesh optimization with applications to hyperbolic conservation laws.
Journal of Computational Physics, page 112924, 2024.

[Fidkowski and Chen, 2021]
Krzysztof J Fidkowski and Guodong Chen.
Metric-based, goal-oriented mesh adaptation using machine learning.
Journal of Computational Physics, 426:0 109957, 2021.

[Fidkowski and Darmofal, 2011]
Krzysztof J Fidkowski and David L Darmofal.
Review of output-based error estimation and mesh adaptation in computational fluid dynamics.
AIAA journal, 490 (4):0 673–694, 2011.

[Fortunato et al., 2022]
Meire Fortunato, Tobias Pfaff, Peter Wirnsberger, Alexander Pritzel, and Peter Battaglia.
Multiscale meshgraphnets.
In ICML 2022 2nd AI for Science Workshop, 2022.
URL <https://openreview.net/forum?id=G3TRIsmMhhf>.

[Foucart et al., 2023]
Corbin Foucart, Aaron Charous, and Pierre FJ Lermusiaux.
Deep reinforcement learning for adaptive mesh refinement.
Journal of Computational Physics, 491:0 112381, 2023.

[Freymuth et al., 2023]
Niklas Freymuth, Philipp Dahlinger, Tobias Würth, Simon Reisch, Luise Kärger, and Gerhard Neumann.
Swarm reinforcement learning for adaptive mesh refinement.
Advances in Neural Information Processing Systems, 36, 2023.

[Gibert et al., 2019]
Gaël Gibert, Benoit Prabel, Anthony Gravouil, and Clémentine Jacquemoud.
A 3d automatic mesh refinement x-fem approach for fatigue crack propagation.
Finite Elements in Analysis and Design, 157:0 21–37, 2019.

[Gillette et al., 2022]
Andrew Gillette, Brendan Keith, and Socratis Petrides.
Learning robust marking policies for adaptive mesh refinement.
arXiv preprint arXiv:2207.06339, 2022.

[Guillet et al., 2019]
Thomas Guillet, Rüdiger Pakmor, Volker Springel, Praveen Chandrashekar, and Christian Klingenberg.
High-order magnetohydrodynamics for astrophysics with an adaptive mesh refinement discontinuous galerkin scheme.
Monthly Notices of the Royal Astronomical Society, 4850 (3):0 4209–4246, 2019.

[Guo et al., 2016]
Xiaoxiao Guo, Wei Li, and Francesco Iorio.
Convolutional neural networks for steady flow approximation.
In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD '16, page 481–490, New York, NY, USA, 2016. Association for Computing Machinery.
ISBN 9781450342322.
URL <https://doi.org/10.1145/2939672.2939738>.

[Gustafsson and Mcbain, 2020]
Tom Gustafsson and Geordie Drummond Mcbain.
scikit-fem: A python package for finite element assembly.
Journal of Open Source Software, 50 (52):0 2369, 2020.

[Hayes et al., 2022]
Conor F Hayes, Roxana Rădulescu, Eugenio Bargiacchi, Johan Källström, Matthew Macfarlane, Mathieu Reymond, Timothy Verstraeten, Luisa M Zintgraf, Richard Dazeley, Fredrik Heintz, et al.
A practical guide to multi-objective reinforcement learning and planning.
Autonomous Agents and Multi-Agent Systems, 36:0 26, 2022.

[He et al., 2016]
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.
Deep residual learning for image recognition.
In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.

[Hessel et al., 2018]
Matteo Hessel, Joseph Modayil, Hado Van Hasselt, Tom Schaul, Georg Ostrovski, Will Dabney, Dan Horgan, Bilal Piot, Mohammad Azar, and David Silver.
Rainbow: Combining improvements in deep reinforcement learning.
In Proceedings of the AAAI conference on artificial intelligence, volume 32, 2018.

[Ho-Le, 1988]
K Ho-Le.
Finite element mesh generation methods: a review and classification.
Computer-aided design, 200 (1):0 27–38, 1988.

[Huang et al., 2021]
Keefe Huang, Moritz Krügener, Alistair Brown, Friedrich Menhorn, Hans-Joachim Bungartz, and Dirk Hartmann.
Machine learning-based optimal mesh generation in computational fluid dynamics.
arXiv preprint arXiv:2102.12923, 2021.

[Huang and Russell, 2010]
Weizhang Huang and Robert D Russell.
Adaptive moving mesh methods, volume 174.
Springer Science & Business Media, 2010.

[Huergo et al., 2024]
David Huergo, Gonzalo Rubio, and Esteban Ferrer.
A reinforcement learning strategy for p-adaptation in high order solvers.
Results in Engineering, 21:0 101693, 2024.

[Hüttenrauch et al., 2019]
Maximilian Hüttenrauch, Šošić Adrian, and Gerhard Neumann.
Deep reinforcement learning for swarm systems.
Journal of Machine Learning Research, 200 (54):0 1–31, 2019.

[Ioffe and Szegedy, 2015]
Sergey Ioffe and Christian Szegedy.
Batch normalization: Accelerating deep network training by reducing internal covariate shift.
In International conference on machine learning, pages 448–456. pmlr, 2015.

[John et al., 2016]
Volker John et al.
Finite element methods for incompressible flow problems, volume 51.
Springer, 2016.

[Jones and Plassmann, 1997]
Mark T Jones and Paul E Plassmann.
Adaptive refinement of unstructured finite-element meshes.
Finite Elements in Analysis and Design, 250 (1-2):0 41–60, 1997.

[Kingma and Ba, 2014]
Diederik P Kingma and Jimmy Ba.
Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980, 2014.

[Kita and Kamiya, 2001]
Eisuke Kita and Norio Kamiya.
Error estimation and adaptive mesh refinement in boundary element method, an overview.
Engineering Analysis with Boundary Elements, 250 (7):0 479–495, 2001.

[Larson and Bengzon, 2013]
Mats G. Larson and Fredrik Bengzon.
The Finite Element Method: Theory, Implementation, and Applications, volume 10 of Texts in Computational Science and Engineering.
Springer, Berlin, Heidelberg, 2013.
ISBN 978-3-642-33286-9 978-3-642-33287-6.

[Linkerhägner et al., 2023]
Jonas Linkerhägner, Niklas Freymuth, Paul Maria Scheikl, Franziska Mathis-Ullrich, and Gerhard Neumann.
Grounding graph network simulators using physical sensor observations.
In The Eleventh International Conference on Learning Representations (ICLR), 2023.
URL <https://openreview.net/forum?id=jsZsEd8VEY>.

[Mnih et al., 2013]
Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin Riedmiller.
Playing atari with deep reinforcement learning.
arXiv preprint arXiv:1312.5602, 2013.

[Mnih et al., 2015]
Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al.
Human-level control through deep reinforcement learning.
nature, 5180 (7540):0 529–533, 2015.

[Mukherjee, 1996]
Arup Mukherjee.
An adaptive finite element code for elliptic boundary value problems in three dimensions with applications in numerical relativity.
The Pennsylvania State University, 1996.

[Nagarajan and Soghrati, 2018]
Anand Nagarajan and Soheil Soghrati.
Conforming to interface structured adaptive mesh refinement: 3d algorithm and implementation.
Computational Mechanics, 62:0 1213–1238, 2018.

[Ortiz and Quigley Iv, 1991]
M Ortiz and JJ Quigley Iv.
Adaptive mesh refinement in strain localization problems.
Computer Methods in Applied Mechanics and Engineering, 900 (1-3):0 781–804, 1991.

[Pan et al., 2023]
Jie Pan, Jingwei Huang, Gengdong Cheng, and Yong Zeng.
Reinforcement learning for automatic quadrilateral mesh generation: A soft actor–critic approach.
Neural Networks, 157:0 288–304, 2023.

[Paszke et al., 2019]
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.
Pytorch: An imperative style, high-performance deep learning library.
Advances in neural information processing systems, 32, 2019.

[Perera and Agrawal, 2024]
Roberto Perera and Vinamra Agrawal.
Multiscale graph neural networks with adaptive mesh refinement for accelerating mesh-based simulations.
arXiv preprint arXiv:2402.08863, 2024.

[Pfaff et al., 2021]
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter W. Battaglia.
Learning mesh-based simulation with graph networks.
In International Conference on Learning Representations, 2021.
URL <https://arxiv.org/abs/2010.03409>.

[Plewa et al., 2005]
Tomasz Plewa, Timur Linde, V Gregory Weirs, et al.
Adaptive mesh refinement-theory and applications.
Springer, 2005.

[Quarteroni and Quarteroni, 2009]
Alfio Quarteroni and Silvia Quarteroni.
Numerical models for differential problems, volume 2.
Springer, 2009.

[Rachowicz et al., 1989]
Waldemar Rachowicz, J Tinsley Oden, and Leszek Demkowicz.
Toward a universal hp adaptive finite element strategy part 3. design of hp meshes.
Computer Methods in Applied Mechanics and Engineering, 770 (1-2):0 181–212, 1989.

[Rachowicz et al., 2006]
Waldemar Rachowicz, David Pardo, and Leszek Demkowicz.
Fully automatic hp-adaptivity in three dimensions.
Computer methods in applied mechanics and engineering, 1950 (37-40):0 4816–4842, 2006.

[Raissi et al., 2019]
Maziar Raissi, Paris Perdikaris, and George E Karniadakis.
Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.
Journal of Computational physics, 378:0 686–707, 2019.

[Reddy, 2019]
Junuthula Narasimha Reddy.
Introduction to the finite element method.
McGraw-Hill Education, 2019.

[Rivara, 1984]
M Cecilia Rivara.
Algorithms for refining triangular grids suitable for adaptive and multigrid techniques.
International journal for numerical methods in Engineering, 200 (4):0 745–756, 1984.

[Rong et al., 2019]
Yu Rong, Wenbing Huang, Tingyang Xu, and Junzhou Huang.
Dropedge: Towards deep graph convolutional networks on node classification.
In International Conference on Learning Representations, 2019.

[Roth et al., 2022]
Julian Roth, Max Schröder, and Thomas Wick.
Neural network guided adjoint computations in dual weighted residual error estimation.
SN Applied Sciences, 40 (2):0 62, 2022.

[Sanchez-Gonzalez et al., 2020]
Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia.
Learning to simulate complex physics with graph networks.
In Proceedings of the 37th International Conference on Machine Learning, pages 8459–8468. PMLR, 2020.

[Scarselli et al., 2009]
Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini.
The graph neural network model.
IEEE Transactions on Neural Networks, 200 (1):0 61–80, 2009.

[Schaul et al., 2016]
Tom Schaul, John Quan, Ioannis Antonoglou, and David Silver.
Prioritized experience replay.
In Yoshua Bengio and Yann LeCun, editors, 4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings, 2016.

[Schulman et al., 2015]
John Schulman, Philipp Moritz, Sergey Levine, Michael Jordan, and Pieter Abbeel.
High-dimensional continuous control using generalized advantage estimation.
arXiv preprint arXiv:1506.02438, 2015.

[Schulman et al., 2017]
John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov.
Proximal policy optimization algorithms.
arXiv preprint arXiv:1707.06347, 2017.

[Shi et al., 2023]
Haochen Shi, Huazhe Xu, Samuel Clarke, Yunzhu Li, and Jiajun Wu.
Robocook: Long-horizon elasto-plastic object manipulation with diverse tools.
In 7th Annual Conference on Robot Learning, 2023.

[Służalec et al., 2023]
Tomasz Służalec, Rafał Grzeszczuk, Sergio Rojas, Witold Dzwinel, and Maciej Paszyński.
Quasi-optimal hp-finite element refinements towards singularities via deep neural network prediction.
Computers & Mathematics with Applications, 142:0 157–174, 2023.

[Šošić et al., 2017]
Adrian Šošić, Wasiur R KhudaBukhsh, Abdelhak M Zoubir, and Heinz Koeppl.
Inverse reinforcement learning in swarm systems.
In Proceedings of the 16th Conference on Autonomous Agents and MultiAgent Systems, pages 1413–1421, 2017.

[Stein, 2007]
Erwin Stein.
Adaptive finite elements in linear and nonlinear solid and structural mechanics, volume 416.
Springer Science & Business Media, 2007.

[Stevenson, 2008]
Rob Stevenson.
The completion of locally refined simplicial partitions created by bisection.
Mathematics of computation, 770 (261):0 227–241, 2008.

[Strauss, 2007]
Walter A Strauss.
Partial differential equations: An introduction.
John Wiley & Sons, 2007.

[Suárez et al., 2005]
José P Suárez, Pilar Abad, Angel Plaza, and Miguel A Padron.
Computational aspects of the refinement of 3d tetrahedral meshes.
Journal of Computational Methods in Sciences and Engineering, 50 (4):0 215–224, 2005.

[Sunehag et al., 2017]
Peter Sunehag, Guy Lever, Audrunas Gruslys, Wojciech Marian Czarnecki, Vinicius Zambaldi, Max Jaderberg, Marc Lanctot, Nicolas Sonnerat, Joel Z Leibo, Karl Tuyls, et al.
Value-decomposition networks for cooperative multi-agent learning.
arXiv preprint arXiv:1706.05296, 2017.

[Sutton and Barto, 2018]
Richard S Sutton and Andrew G Barto.
Reinforcement learning: An introduction.
MIT press, 2018.

[Um et al., 2018]
Kiwon Um, Xiangyu Hu, and Nils Thuerey.
Liquid splash modeling with neural networks.
Computer Graphics Forum, 370 (8):0 171–182, 2018.
URL <https://onlinelibrary.wiley.com/doi/abs/10.1111/cgf.13522>.

[Van Hasselt et al., 2016]
Hado Van Hasselt, Arthur Guez, and David Silver.
Deep reinforcement learning with double q-learning.
In Proceedings of the AAAI conference on artificial intelligence, volume 30, 2016.

[Van Moffaert and Nowé, 2014]
Kristof Van Moffaert and Ann Nowé.
Multi-objective reinforcement learning using sets of pareto dominating policies.
The Journal of Machine Learning Research, 150 (1):0 3483–3512, 2014.

[Velickovic et al., 2018]
Petar Velickovic, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.
Graph attention networks.
stat, 1050:0 4, 2018.

[Wallwork, 2021]
Joseph Gregory Wallwork.
Mesh adaptation and adjoint methods for finite element coastal ocean modelling.
PhD thesis, Imperial College London, 2021.

[Wallwork et al., 2022]
Joseph Gregory Wallwork, Jingyi Lu, Mingrui Zhang, and Matthew D Piggott.
E2n: Error estimation networks for goal-oriented mesh adaptation.
arXiv preprint arXiv:2207.11233, 2022.

[Wang et al., 2016]
Ziyu Wang, Tom Schaul, Matteo Hessel, Hado Hasselt, Marc Lanctot, and Nando Freitas.
Dueling network architectures for deep reinforcement learning.
In International conference on machine learning, pages 1995–2003. PMLR, 2016.

[Wanner and Hairer, 1996]
Gerhard Wanner and Ernst Hairer.
Solving ordinary differential equations II, volume 375.
Springer Berlin Heidelberg New York, 1996.

[Wu et al., 2023]
Tailin Wu, Takashi Maruyama, Qingqing Zhao, Gordon Wetzstein, and Jure Leskovec.
Learning controllable adaptive simulation for multi-resolution physics.
In The Eleventh International Conference on Learning Representations, 2023.
URL <https://openreview.net/forum?id=PbfgkZ2HdbE>.

[Würth et al., 2023]
Tobias Würth, Constantin Krauß, Clemens Zimmerling, and Luise Kärger.
Physics-informed neural networks for data-free surrogate modelling and engineering optimization–an example from composite manufacturing.
Materials & Design, 231:0 112034, 2023.

[Würth et al., 2024]
Tobias Würth, Niklas Freymuth, Clemens Zimmerling, Gerhard Neumann, and Luise Kärger.
Physics-informed meshgraphnets (pi-mgn): Neural finite element solvers for non-stationary and nonlinear simulations on arbitrary meshes.
arXiv:2402.10681, 2024a.

[Würth et al., 2024]
Tobias Würth, Niklas Freymuth, Clemens Zimmerling, Gerhard Neumann, and Luise Kärger.
Physics-informed meshgraphnets (pi-mgns): Neural finite element solvers for non-stationary and nonlinear simulations on arbitrary meshes.
arXiv preprint arXiv:2402.10681, 2024b.

[Yang et al., 2023]
Jiachen Yang, Tarik Dzanic, Brenden K Petersen, Jun Kudo, Ketan Mittal, Vladimir Tomov, Jean-Sylvain Camier, Tuo Zhao, Hongyuan Zha, Tzanio Kolev, Robert Anderson, and Daniel Faissol.
Reinforcement learning for adaptive mesh refinement.
26th International Conference on Artificial Intelligence and Statistics (AISTATS), 2023a.

[Yang et al., 2023]
Jiachen Yang, Ketan Mittal, Tarik Dzanic, Socratis Petrides, Brendan Keith, Brenden Petersen, Daniel Faissol, and Robert Anderson.
Multi-agent reinforcement learning for adaptive mesh refinement.
22nd International Conference on Autonomous Agents and Multiagent Systems (AAMAS), 2023b.

[Yano and Darmofal, 2012]
Masayuki Yano and David L Darmofal.
An optimization-based framework for anisotropic simplex mesh adaptation.
Journal of Computational Physics, 2310 (22):0 7626–7649, 2012.

[Zhang et al., 2020]
Zheyan Zhang, Yongxing Wang, Peter K Jimack, and He Wang.
Meshingnet: A new mesh generation method based on deep learning.
In Computational Science–ICCS 2020: 20th International Conference, Amsterdam, The Netherlands, June 3–5, 2020, Proceedings, Part III 20, pages 186–198. Springer, 2020.

[Zienkiewicz et al., 2005]
Olek C Zienkiewicz, Robert Leroy Taylor, and Jian Z Zhu.
The finite element method: its basis and fundamentals.
Elsevier, 2005.

[Zienkiewicz and Zhu, 1992]
Olgierd Cecil Zienkiewicz and Jian Zhong Zhu.
The superconvergent patch recovery and a posteriori error estimates. part 1: The recovery technique.
International Journal for Numerical Methods in Engineering, 330 (7):0 1331–1364, 1992.

[Zimmerling et al., 2019]
Clemens Zimmerling, Dominik Dörr, Frank Henning, and Luise Kärger.
A machine learning assisted approach for textile formability assessment and design improvement of composite components.
Compos Part A, 124:0 105459, September 2019a.
ISSN 1359-835X.

[Zimmerling et al., 2019]
Clemens Zimmerling, Daniel Trippe, Benedikt Fengler, and Luise Kärger.
An approach for rapid prediction of textile draping results for variable composite component geometries using deep neural networks.
AIP Conference Proceedings, 21130 (1), 07 2019b.
ISSN 0094-243X.
URL <https://doi.org/10.1063/1.5112512>.

[Zimmerling et al., 2022]
Clemens Zimmerling, Christian Poppe, Oliver Stein, and Luise Kärger.
Optimisation of manufacturing process parameters for variable component geometries using reinforcement learning.
Materials & Design, 214:0 110423, February 2022.
ISSN 02641275.