# Swarm Reinforcement Learning for Adaptive Mesh Refinement

Niklas Freymuth1 Philipp Dahlinger1 Tobias Würth2 Luise Kärger2 Gerhard
Neumann1
1Autonomous Learning Robots, Karlsruhe Institute of Technology, Karlsruhe
2Institute of Vehicle Systems Technology, Karlsruhe Institute of Technology,
Karlsruhe correspondence to<EMAIL_ADDRESS>
###### Abstract

Adaptive Mesh Refinement (AMR) is crucial for mesh-based simulations, as it
allows for dynamically adjusting the resolution of a mesh to trade off
computational cost with the simulation accuracy. Yet, existing methods for AMR
either use task-dependent heuristics, expensive error estimators, or do not
scale well to larger meshes or more complex problems. In this paper, we
formalize AMR as a Swarm Reinforcement Learning problem, viewing each element
of a mesh as part of a collaborative system of simple and homogeneous agents.
We combine this problem formulation with a novel agent-wise reward function
and Graph Neural Networks, allowing us to learn reliable and scalable
refinement strategies on arbitrary systems of equations. We experimentally
demonstrate the effectiveness of our approach in improving the accuracy and
efficiency of complex simulations. Our results show that we outperform learned
baselines and achieve a refinement quality that is on par with a traditional
error-based AMR refinement strategy without requiring error indicators during
inference.

##  1 Introduction

Adaptive Mesh Refinement (AMR) is a powerful technique for complex Finite
Element Method (FEM)-based simulations (Reddy, 2019; Sabat & Kundu, 2021) that
allows to dynamically allocate resources to relevant parts of the simulated
domain to trade off computational speed and simulation accuracy (Plewa et al.,
2005; Anderson et al., 2021). Applications of AMR include fluid dynamics
(Berger & Colella, 1989; Baker, 1997; Borker et al., 2019), structural
mechanics (Ortiz & Quigley Iv, 1991; Provatas et al., 1998; Stein, 2007;
Gibert et al., 2019) and astrophysics (Cunningham et al., 2009; Bryan et al.,
2014; Guillet et al., 2019). Yet, classical approaches for AMR usually rely on
problem-dependent heuristics or error indicators (Mukherjee, 1996; Arnold et
al., 2000; Kita & Kamiya, 2001; Bangerth & Rannacher, 2003; Cerveny et al.,
2019). Here, we instead formalize AMR as a Reinforcement Learning (RL) (Sutton
& Barto, 2018) problem. Following previous work (Yang et al., 2021; Foucart et
al., 2022; Yang et al., 2022) in this direction, we encode the state of the
current simulation as local observations that we feed to RL agents. The chosen
actions then determine which elements of a mesh to refine. Opposed to e.g.,
Physics-Informed Neural Networks (Raissi et al., 2019; Cai et al., 2021),
which are neural networks that are trained to satisfy the governing equations
of a physical system, learned AMR strategies combine Machine Learning with
classical mesh-based methods. As such, learned AMR methods present a more
robust approach to solving complex physics problems that require high accuracy
and precision, as they identify regions of interest which are subsequently
used for a classical numerical solution.

While existing RL-based AMR approaches are trained on refinements of
individual elements (Yang et al., 2021; Foucart et al., 2022) or use a learned
value decomposition of parallel refinements (Yang et al., 2022), we instead
employ Swarm Reinforcement Learning (Šošić et al., 2017; Hüttenrauch et al.,
2019) to efficiently train reliable and scalable refinement strategies. For
this, our approach uses a novel spatial reward formulation that assigns each
element of the mesh an agent and each agent its own reward throughout an
episode. We combine this reward with Message Passing Networks (MPNs) (Sanchez-
Gonzalez et al., 2020), a class of Graph Neural Networks (GNNs) (Scarselli et
al., 2009; Battaglia et al., 2018; Wu et al., 2020) that has proven highly
effective for physical simulations (Pfaff et al., 2021; Linkerhägner et al.,
2023). The resulting method, Adaptive Swarm Mesh Refinement (ASMR), produces
highly refined meshes with thousands of elements while being applicable to
arbitrary Partial Differential Equations (PDEs). Experimentally, we show the
effectiveness of ASMR on challenging elliptic PDEs, outperforming strong
learned baselines (Yang et al., 2021; 2022). Further, our results are on par
with a traditional error-based AMR heuristic without requiring expensive error
indicators during inference. We additionally conduct a series of ablations to
show which parts of the approach make it uniquely effective.

Figure 1:  Overview of Adaptive Swarm Mesh Refinement (ASMR). Given a mesh
$\Omega_{t}$, the observation graph ${\mathcal{G}}_{\Omega^{t}}$ encodes the
elements as graph nodes, their neighborhood as edges, and general information
as global features. The graph is given to a learned policy, which marks mesh
elements for refinement. A remesher refines these elements, resulting in a
refined mesh.

##  2 Adaptive Swarm Mesh Refinement

Message Passing Networks. Let
${\mathcal{G}}=(\mathbf{V},\mathbf{E},\mathbf{X}_{\mathbf{V}},\mathbf{X}_{\mathbf{E}},\mathbf{g})$
be a directed graph with nodes $\mathbf{V}$ and edges
$\mathbf{E}\subseteq\mathbf{V}\times\mathbf{V}$. Then
$\mathbf{X}_{\mathbf{V}}:\mathbf{V}\rightarrow\mathbb{R}^{d_{\mathbf{V}}}$ and
$\mathbf{X}_{\mathbf{E}}:\mathbf{E}\rightarrow\mathbb{R}^{d_{\mathbf{E}}}$ are
node and edge features of dimensions $d_{\mathbf{V}}$ and $d_{\mathbf{E}}$
respectively, and $\mathbf{g}\in\mathbb{R}^{d_{\mathbf{g}}}$ are global
features. A MPN (Sanchez-Gonzalez et al., 2020; Pfaff et al., 2021;
Linkerhägner et al., 2023) is a GNN consisting of $L$ Message Passing Steps.
Each step $l$ receives the output of the previous step and updates the
features for all nodes $v\in\mathbf{V}$, edges $e\in\mathbf{E}$, and globals
$\mathbf{g}$. Using linear embeddings $\mathbf{x}_{v}^{0}$,
$\mathbf{x}_{e}^{0}$, and $\mathbf{g}^{0}$ of the initial node, edge, and
global features, the step is given as

$\mathbf{x}^{l+1}_{e}=f^{l}_{\mathbf{E}}(\mathbf{x}^{l}_{v},\mathbf{x}^{l}_{u},\mathbf{x}^{l}_{e},\mathbf{g}^{l}),\textrm{
with }e=(u,v)\text{,}$
$\mathbf{x}^{l+1}_{v}=f^{l}_{\mathbf{V}}(\mathbf{x}^{l}_{v},\bigoplus_{e=(v,u)\in\mathbf{E}}\mathbf{x}^{l+1}_{e},\mathbf{g}^{l})\text{,}\quad\text{~{}and~{}}\quad\mathbf{g}^{l+1}=f^{l}_{\mathbf{g}}(\bigoplus_{v\in\mathbf{V}}\mathbf{x}^{l+1}_{v},\bigoplus_{e\in\mathbf{E}}\mathbf{x}^{l+1}_{e},\mathbf{g}^{l})\text{.}$

The operator $\oplus$ is a permutation-invariant aggregation such as a sum,
max, or mean operator. Each $f^{l}_{\cdot}$ is a learned function that we
generally parameterize as a simple Multilayer Perceptron (MLP). The network’s
final output is a learned representation $\mathbf{x}^{L}_{v}$ for each node
$v\in\mathbf{V}$. Appendix A.1 provides a schematic overview of the MPN
architecture.

Swarm Markov Decision Process We frame AMR as an instance of a Swarm Markov
Decision Process (SwarmDP) (Šošić et al., 2017; Hüttenrauch et al., 2019),
which can be seen as a special case of a decentralized partially observable
MDP for swarm systems. For this, we adapt the process to per-agent rewards, a
shared observation structure, and varying numbers of agents. Details can be
found in Appendix A.2.

Agents and Actions. Given a mesh $\Omega^{t}$ for a domain $\Omega$, we view
each element $\Omega^{t}_{i}\subseteq\Omega^{t}$ as an agent. Each element’s
action space comprises a binary decision to mark it for refinement or not. A
remesher refines all marked elements $i$ as
$\Omega_{i}^{t}=\bigcup_{j}\Omega_{i,j}^{t+1}$, resulting in a finer mesh
$\Omega^{t+1}=\bigcup_{j}\Omega_{j}^{t+1}$. The remesher may also refine
unmarked elements to assert a conforming solution (Arnold et al., 2000).

Observations. We encode an observation graph ${\mathcal{G}}_{\Omega^{t}}$ by
creating a node for each mesh element and (directed) edges between all pairs
of neighboring elements. The node, edge, and global features of the resulting
graph depend on the considered system of equations. However, we encode
positions on the mesh only as relative Euclidean distances in the edge
features $\mathbf{X}_{\mathbf{E}}$ to ensure that our observations are
equivariant under the Euclidean group (Bronstein et al., 2021; Pfaff et al.,
2021).

Figure 2: Qualitative results for ASMR on Poisson’s Equation on an L-shaped
domain. Ranging from high (left) to low (right) element penalties $\alpha$,
our approach produces highly accurate refinements for different numbers of
mesh elements. In all cases, the policy focuses on and refines the relevant
parts of the underlying system of equations. We list more examples in Appendix
E.

Reward. A good refinement strategy trades off the accuracy of the mesh
$\Omega^{t}$ with its total number of elements
$\Omega^{t}_{i}\subseteq\Omega^{t}$. Let $\tilde{\text{err}}(\Omega^{t}_{i})$
be an error estimate of element $\Omega^{t}_{i}$,
$\text{Area}(\Omega^{t}_{i})$ its total area and
$\delta^{t}(\Omega^{t}_{i})=\\{\Omega^{t+1}_{j}\subseteq\Omega^{t}_{i}\\}$ the
mapping that assigns the element to all elements $\Omega^{t+1}_{j}$ that it
refines into. We then define a reward per element as

$R(\Omega^{t}_{i})=\frac{1}{\text{Area}(\Omega^{t}_{i})}\left(\tilde{\text{err}}(\Omega^{t}_{i})-\sum_{\Omega_{j}^{t+1}\in\delta^{t}(\Omega^{t}_{i})}\tilde{\text{err}}(\Omega^{t+1}_{j})\right)-\alpha\left(|\delta^{t}(\Omega_{i}^{t})|-1\right)\text{,}$
(1)

where $\alpha$ is a hyperparameter that acts as a penalty for adding new
elements. Note that the reference mesh $\Omega^{*}$ is only needed during
training and that the trained policy does not depend on it when refining a new
mesh during inference. To track the influence of the agents even though they
may appear and disappear throughout an episode, we can iteratively assign the
mapping $\delta^{t}(\Omega^{t}_{i})$ for $k$ time steps as
$\delta^{t}_{k}(\Omega^{t}_{i})=\\{\Omega^{t+k}_{j}\subseteq\Omega^{t}_{i}\\}$.
This mapping allows us to compute a return per element

$J^{t}_{i}=\mathbb{E}_{\pi}\left[\sum_{k=0}^{\infty}\gamma^{k}\sum_{\Omega^{t+k}_{j}\in\delta^{t}_{k}(\Omega^{t}_{i})}R(\Omega^{t+k}_{j})\right]\text{.}$

We provide additional details on the construction of the error estimate
$\tilde{\text{err}}(\Omega^{t}_{i})$, as well as an intuition on the reward in
Appendix A.3.

Adaptive Swarm Mesh Refinement. Our approach trains a policy
$\pi(\mathbf{a}|{\mathcal{G}}_{\Omega^{t}})$ that computes a joint action
vector $\mathbf{a}\in\mathcal{A}^{|\\{\Omega^{t}_{i}\in\Omega^{t}\\}|}$ for
all mesh elements $\Omega^{t}_{i}$ by passing messages along the observation
graph ${\mathcal{G}}$. Figure 1 gives a schematic overview.

##  3 Experiments

Systems of Equations and Graph Features. We experiment on two families of
elliptical PDEs, namely Laplace’s Equation on domains with an inner and an
outer boundary with different boundary conditions and Poisson’s Equation on
L-shaped domains for Gaussian Mixture Model load functions. Both equations are
ubiquitous in physics (Strauss, 2007), with applications including heat flow
(Zill, 2016), electrodynamics (Griffiths, 2005) and astrophysics (Misner et
al., 1973). We refer to Appendix B for further details. The nodes features
$\mathbf{X}_{v}$ for each node $v\in\mathbf{V}$ are given as the element area,
the distance to the closest boundary, and the mean and standard deviation of
the solution on the element’s vertices. Edge features are defined as euclidean
distances between element midpoints, and global features $\mathbf{g}$ consist
of the number of mesh elements and mesh vertices, as well as the current
environment step. For Laplace’s Equation, the closest distance to the inner
boundary is included as an additional node feature, while for Poisson’s
Equation, the load function $f$ is evaluated at the respective face midpoint.

Setup. We use Proximal Policy Optimization (PPO) (Schulman et al., 2017) for
all experiments. We experiment on different $2$d domains with triangular
meshes and linear elements. The PDEs are implemented using scikit-fem
(Gustafsson & Mcbain, 2020). All methods are evaluated for different $\alpha$
to produce a wide range of solutions, and we repeat each experiment
configuration for $10$ different random seeds. We report the average
performance on $100$ randomly sampled but fixed evaluation PDEs for each seed.
Appendix D.1 lists all PPO and network training hyperparameters.

All environment episodes start with a coarse initial mesh $\Omega^{0}$ and
iteratively refine this mesh $T=4$ times unless mentioned otherwise. During
training, the reference mesh $\Omega^{*}$ is computed by uniformly refining
the initial mesh $\Omega^{0}$ a total of $4$ times. We evaluate a given
refinement via its normalized numerically integrated squared error compared to
a $5$ times refined reference mesh to ensure a small numerical error in the
evaluation. Details on this metric are given in Appendix C.2. We also provide
results for a normalized linear error in Appendix C.3.

(a) Laplace’s Equation

(b) Poisson’s Equation

Figure 3: Normalized squared error plotted against the number of elements for
(3(a)) the Laplace Equation and (3(b)) the Poisson Equation. Each point
represents the average performance over $100$ evaluation environments for one
random seed. The lines are a Log-Log linear regression over the aggregated
results for each method. For both experiments, all methods can provide
refinements for different numbers of elements. However, ASMR is significantly
more stable than the other learned methods, producing more consistent and
higher-quality refinements, especially for larger meshes.

Baselines. We adapt an Argmax (Yang et al., 2021) and a VDGN-inspired (Yang et
al., 2022) baseline, called PPO-VDQN, to our setup. Detailed information on
the baselines is given in Appendix C.1. We also compare to a traditional
error-based AMR baseline (Binev et al., 2004; Bangerth et al., 2012; Foucart
et al., 2022). Given a parameter $\theta$, this Heuristic iteratively refines
all elements $\Omega^{t}_{i}$ for which
$\tilde{\text{err}}(\Omega^{t}_{i})>\max_{j}\tilde{\text{err}}(\Omega^{t}_{j})$.
This baseline uses the reference mesh $\Omega^{*}$ for error calculation as
oracle information, which is usually unavailable during inference.
Additionally, the heuristic is purely local and thus may be sub-optimal for
globally propagating errors, which is a well-known issue for elliptic PDEs
(Strauss, 2007; Foucart et al., 2022). Opposed to this, RL methods learn to
maximize an expected return, allowing them to find more strategic marking
strategies.

We additionally conduct ablations to determine which parts of our approach
make it effective. For this, we consider an Unscaled variant for our method
that drops the area scaling in Equation 1, as well as a Shared Reward version
that considers the sum of all agent’s rewards as a single global reward. Note
that the combination of these two changes results in PPO-VDGN.

Results. Figure 2 shows qualitative results of ASMR on an exemplary Poisson
problem. For different values of $\alpha$, our approach provides refinements
of different granularity. Quantitatively, Figure 3 shows that our approach
outperforms other learned methods, especially for larger meshes.
Interestingly, ASMR matches or surpasses the error-based heuristic baseline,
indicating its ability to learn sophisticated refinement strategies that
minimize error throughout the entire environment episode, even without access
to error estimates. We provide additional evaluations for a linear instead of
a squared error in Appendix C.3. Appendix C.4 provides further ablation
experiments on different observation graphs. Finally, Appendix E lists
additional visualizations for all methods.

##  4 Conclusion

We propose Adaptive Swarm Mesh Refinement (ASMR), a novel method for Adaptive
Mesh Refinement that utilizes Swarm Reinforcement Learning to iteratively
refine meshes to efficiently solve Partial Differential Equations. By viewing
the mesh elements as agents in a swarm system, our method can efficiently
train a shared policy for all agents. Combining this approach with Graph
Neural Networks and a novel per-agent reward formulation, our method yields
stable and efficient mesh refinements for meshes with thousands of elements
without requiring access to an error estimate during inference.
Experimentally, we show that our method outperforms existing Reinforcement
Learning mesh refinement techniques and compares well to an error-based
heuristic. In future work, we want to extend our approach to more complex
systems of equations, including time-dependent Partial Differential Equations.
Other promising directions include extending our approach to quadrilateral
meshes and combining the current per-agent reward with a global measure of
mesh quality.

#### Acknowledgments

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.

## References

  * 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: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_ , 22(2):431–448, 2000.
  * Ba et al. (2016) Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton.  Layer normalization.  _stat_ , 1050:21, 2016.
  * Baker (1997) Timothy J Baker.  Mesh adaptation strategies for problems in fluid dynamics.  _Finite Elements in Analysis and Design_ , 25(3-4):243–273, 1997.
  * Bangerth & Rannacher (2003) Wolfgang Bangerth and Rolf Rannacher.  _Adaptive finite element methods for differential equations_.  Springer Science & Business Media, 2003.
  * 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)_ , 38(2):1–28, 2012.
  * Battaglia et al. (2018) Peter W. Battaglia, Jessica B. Hamrick, Victor Bapst, Alvaro Sanchez-Gonzalez, Vinícius Flores Zambaldi, Mateusz Malinowski, Andrea Tacchetti, David Raposo, Adam Santoro, Ryan Faulkner, Çaglar Gülçehre, H. Francis Song, Andrew J. Ballard, Justin Gilmer, George E. Dahl, Ashish Vaswani, Kelsey R. Allen, Charles Nash, Victoria Langston, Chris Dyer, Nicolas Heess, Daan Wierstra, Pushmeet Kohli, Matthew Botvinick, Oriol Vinyals, Yujia Li, and Razvan Pascanu.  Relational inductive biases, deep learning, and graph networks.  _CoRR_ , abs/1806.01261, 2018.  URL http://arxiv.org/abs/1806.01261.
  * Berger & Colella (1989) Marsha J Berger and Phillip Colella.  Local adaptive mesh refinement for shock hydrodynamics.  _Journal of computational Physics_ , 82(1):64–84, 1989.
  * Binev et al. (2004) Peter Binev, Wolfgang Dahmen, and Ron DeVore.  Adaptive finite element methods with convergence rates.  _Numerische Mathematik_ , 97:219–268, 2004.
  * Borker et al. (2019) Raunak Borker, Daniel Huang, Sebastian Grimberg, Charbel Farhat, Philip Avery, and Jason Rabinovitch.  Mesh adaptation framework for embedded boundary methods for computational fluid dynamics and fluid-structure interaction.  _International Journal for Numerical Methods in Fluids_ , 90(8):389–424, 2019.
  * 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_ , 211(2):19, 2014.
  * Cai et al. (2021) Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis.  Physics-informed neural networks (pinns) for fluid mechanics: A review.  _Acta Mechanica Sinica_ , 37(12):1727–1738, 2021\.
  * Cerveny et al. (2019) Jakub Cerveny, Veselin Dobrev, and Tzanio Kolev.  Nonconforming mesh refinement for high-order finite elements.  _SIAM Journal on Scientific Computing_ , 41(4):C367–C392, 2019.
  * 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_ , 182(2):519, 2009.
  * Foucart et al. (2022) Corbin Foucart, Aaron Charous, and Pierre FJ Lermusiaux.  Deep reinforcement learning for adaptive mesh refinement.  _arXiv preprint arXiv:2209.12351_ , 2022.
  * 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:21–37, 2019\.
  * Griffiths (2005) David J Griffiths.  Introduction to electrodynamics, 2005.
  * 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_ , 485(3):4209–4246, 2019.
  * Gustafsson & Mcbain (2020) Tom Gustafsson and Geordie Drummond Mcbain.  scikit-fem: A python package for finite element assembly.  _Journal of Open Source Software_ , 5(52):2369, 2020.
  * 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_ , pp. 770–778, 2016.
  * Hüttenrauch et al. (2019) Maximilian Hüttenrauch, Šošić Adrian, Gerhard Neumann, et al.  Deep reinforcement learning for swarm systems.  _Journal of Machine Learning Research_ , 20(54):1–31, 2019.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kita & Kamiya (2001) Eisuke Kita and Norio Kamiya.  Error estimation and adaptive mesh refinement in boundary element method, an overview.  _Engineering Analysis with Boundary Elements_ , 25(7):479–495, 2001.
  * 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.
  * Misner et al. (1973) Charles W Misner, Kip S Thorne, and John Archibald Wheeler.  _Gravitation_.  Macmillan, 1973.
  * 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.
  * Ortiz & Quigley Iv (1991) M Ortiz and JJ Quigley Iv.  Adaptive mesh refinement in strain localization problems.  _Computer Methods in Applied Mechanics and Engineering_ , 90(1-3):781–804, 1991.
  * 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.
  * 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.  2005\.
  * Provatas et al. (1998) Nikolas Provatas, Nigel Goldenfeld, and Jonathan Dantzig.  Efficient computation of dendritic microstructures using adaptive mesh refinement.  _Physical Review Letters_ , 80(15):3308, 1998\.
  * 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:686–707, 2019\.
  * Reddy (2019) Junuthula Narasimha Reddy.  _Introduction to the finite element method_.  McGraw-Hill Education, 2019.
  * Sabat & Kundu (2021) Lovely Sabat and Chinmay Kumar Kundu.  History of finite element method: a review.  _Recent Developments in Sustainable Infrastructure_ , pp. 395–404, 2021.
  * 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_ , pp. 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_ , 20(1):61–80, 2009.  doi: 10.1109/TNN.2008.2005605.
  * 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.
  * Š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_ , pp. 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.
  * Strauss (2007) Walter A Strauss.  _Partial differential equations: An introduction_.  John Wiley & Sons, 2007.
  * 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 & Barto (2018) Richard S Sutton and Andrew G Barto.  _Reinforcement learning: An introduction_.  MIT press, 2018.
  * Wu et al. (2020) Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip.  A comprehensive survey on graph neural networks.  _IEEE transactions on neural networks and learning systems_ , 32(1):4–24, 2020.
  * Yang et al. (2021) Jiachen Yang, Tarik Dzanic, Brenden Petersen, Jun Kudo, Ketan Mittal, Vladimir Tomov, Jean-Sylvain Camier, Tuo Zhao, Hongyuan Zha, Tzanio Kolev, et al.  Reinforcement learning for adaptive mesh refinement.  _arXiv preprint arXiv:2103.01342_ , 2021.
  * Yang et al. (2022) 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.  _arXiv preprint arXiv:2211.00801_ , 2022.
  * Zill (2016) Dennis G Zill.  _Differential equations with boundary-value problems_.  Cengage Learning, 2016.

##  Appendix A Method Details

###  A.1 Message Passing Network Architecture

Figure 4 shows the Message Passing Network employed by ASMR.

Figure 4: A detailed view of the MPN architecture employed by ASMR. Given a
graph ${\mathcal{G}}$, the global features $\mathbf{g}$, node features
$\mathbf{X}_{\mathbf{V}}$ and edge features $\mathbf{X}_{\mathbf{E}^{\prime}}$
are linearly embedded into a latent space. From this latent space, $L$ message
passing steps are performed. The resulting latent features per node are
interpreted as a local observation encoding and can be given to an RL policy
or value function.

###  A.2 Swarm Markov Decision Process

Formally, we define an abstract agent as a tuple
$\mathbb{A}:=\langle\mathcal{S},\mathcal{O},\mathcal{A},\pi\rangle$. Here,
$\mathcal{S}$, $\mathcal{O}$, $\mathcal{A}$ are sets of local states,
observations and actions, respectively and
$\pi:\mathcal{O}\times\mathcal{A}\rightarrow[0,1]$ is the local policy of the
agent. From this definition, the Swarm MDP is characterized by a tuple
$\langle N,\mathbb{A},T,R,\xi\rangle$. More precisely, given a number of
agents $N$, the system follows a global transition model
$T:\mathcal{S}^{N}\times\mathcal{A}^{N}\rightarrow\mathcal{S}^{M}$, where
potentially $M\neq N$ and uses a per-agent reward function
$R:\mathcal{S}^{N}\times\mathcal{A}^{N}\rightarrow\mathbb{R}^{N}$. The
observations of the agents are calculated from their states via the
observation model $\xi:\mathcal{S}^{N}\rightarrow\mathcal{O}^{N}$.

In our case, the local states and observations share a lot of structure. Given
a mesh $\Omega^{t}$ with $N$ elements $\Omega^{t}_{i}\subseteq\Omega^{t}$,
each element can be considered an agent. The agent’s state then is its local
view on the full mesh, and its observations is a local encoding of the mesh.
We make use of this shared structure by encoding the observations as a graph,
where each node represents a mesh element and each edge the neighborhood of
two elements. This allows us to compute joint actions for all agents by using
the whole observation graph as the input for an MPN-based policy.

###  A.3 Reward Design

Adaptive Mesh Refinement fundamentally trades off the accuracy of the refined
mesh with its usage of resources. To estimate the accuracy, we define an error
per element as the numerically integrated difference in solution of this
element when compared to some fine-grained reference mesh $\Omega^{*}$. For
this, we consider the midpoints $p_{\Omega_{r}^{*}}\in\Omega^{*}_{r}$ of all
elements $\Omega^{*}_{r}$ of the reference mesh. For each element
$\Omega^{t}_{i}$ we then integrate over all midpoints that fall into it,
scaling each by the area of its respective element in $\Omega^{*}$. This
results in an error estimate

$\hat{\text{err}}(\Omega_{i}^{t})\approx\sum_{\Omega_{r}^{*}\subseteq\Omega_{i}^{t}}\text{Area}(\Omega^{*}_{r})\Big{(}|u_{\Omega^{*}}(p_{\Omega_{r}^{*}})-u_{\Omega^{t}}(p_{\Omega_{r}^{*}})|\Big{)}$

for the solution $u_{\Omega^{*}}$ on the fine mesh and the solution
$u_{\Omega^{t}}$ on the current mesh. To get an error estimate that is
consistent across different geometries, we normalize the error with the total
error of the initial mesh $\Omega^{0}$, i.e.,

$\tilde{\text{err}}(\Omega_{i}^{t})=\frac{\hat{\text{err}}(\Omega_{i}^{t})}{\sum_{\Omega_{e}^{0}\in\Omega^{0}}\hat{\text{err}}(\Omega_{e}^{0})}\text{.}$

From here, we formulate a reward per element as given in Equation 1. This
reward intuitively evaluates whether a given refinement decreases the overall
error enough to warrant the extra use of resources, and becomes $0$ for
unrefined faces. The added area scaling term leads to more accurate
refinements, as it discourages the refinements of large elements with a
relatively low average error. Mathematically, optimizing this reward
corresponds to minimizing the error in elements with a high error density
while keeping the total usage of elements low. We expect a high reward of the
refinement of elements with a high error density, because the error is locally
high compared to the fine solution. A refinement thus leads to a more accurate
solution in this area. We note that this reward formulation is purely local,
allowing the approach to scale to a large number of agents, at the cost of
providing less accurate rewards for PDEs with global dependencies.

##  Appendix B Systems of Equations

In its most general form, the FEM is used to approximate the solution $u(x)$
for the set of test functions $\psi$, which satisfies the weak form $\forall
x:~{}\forall\psi(x):~{}a(u(x),\psi(x))=l(\psi(x))$ of the underlying system of
equations.

###  B.1 Laplace’s Equation

Let $\Omega$ be a domain with an inner boundary $\partial\Omega_{0}$ and an
outer boundary $\partial\Omega_{1}$. Laplace’s Equation is given as

$\nabla^{2}u(x)=0\quad\text{in }\Omega\text{.}$

The weak solution to this equation consists of finding $u(x)$ such that

$\int_{\Omega}\nabla u(x)\cdot\nabla v(x)~{}\text{d}x=0$

for all test functions $\psi(x)$. Additionally, the solution has to satisfy
the Dirichlet boundary conditions

$u(x)=0,x\in\partial\Omega_{0}\text{,}\qquad
u(x)=1,x\in\partial\Omega_{1}\text{.}$

We use a unit square $(0,1)^{2}$ for the outer boundary $\partial\Omega_{0}$
of the domain and add a randomly sampled square hole, whose borders are
considered to be the inner boundary $\partial\Omega_{1}$. The size of the hole
is sampled from the uniform distribution $U(0.05,0.25)^{2}$, and its mean
position is sampled from $U(0.2,0.8)^{2}$.

###  B.2 Poisson’s Equation

We additionally consider the Poisson problem

$\nabla^{2}u(x)=f(x)\quad\text{in }\Omega\text{,}\qquad u(x)=0\quad\text{on
}\partial\Omega$

for a load function $f(x):\Omega\rightarrow\mathbb{R}$. Here, the weak
solution consists of finding $u(x)$ such that

$\displaystyle\int_{\Omega}\nabla
u(x)\cdot\nabla\psi(x)~{}\text{d}x=\int_{\Omega}f(x)\psi(x)~{}\text{d}x$

for all test functions $\psi(x)$ while being $0$ on the boundary
$\partial\Omega$.

We model Poisson’s Equation on L-shaped domains, using a rectangular cutoff
whose lower left corner is sampled from $p_{0}\sim U(0.2,0.95)^{2}$, resulting
in a domain $\Omega=(0,1)^{2}\backslash(p_{0}\times(1,1))$. On this domain, we
sample a Gaussian Mixture Model with $o=3$ components. The mean of each
component is sampled from $U(0.1,0.9)^{2}$, using rejection sampling to ensure
that all means lie within the domain. The components’ covariances are
determined by first drawing diagonal covariances, where each dimension is
drawn independently from a log-uniform distribution
$\exp(U(\log(0.0003,0.003)))$. The diagonal covariances are then rotated by a
random angle in $U(0,180)$ to produce Gaussians with a full covariance matrix.
The component weights are drawn from the distribution $\exp(N(0,1))+1$ and
subsequently normalized, where the $1$ in the end is used to ensure that all
components have relevant weight.

##  Appendix C Experiments and Evaluations

###  C.1 Reinforcement Learning Baselines

Both the Argmax and the PPO-VDGN baseline utilize a global reward, which is
the sum over all local rewards without area scaling, and use our MPN for the
policy and value functions. For the Argmax baseline, each element predicts a
continuous action, and in each environment step the element with the highest
action is marked. We remove the face penalty and train for up to $100$
environment steps since the method refines only one element at a time. For
PPO-VDQN, we compute a value function for each agent and then sum them up to
calculate a value decomposition (Sunehag et al., 2017).

###  C.2 Evaluation Metric

Given a reference mesh $\Omega^{*}$ and a final refined mesh $\Omega^{T}$, we
compute the normalized squared error by evaluating the squared difference in
the reference solution $u_{\Omega^{*}}$ and the solution $u_{\Omega^{T}}$ as

$\frac{\sum_{m}\text{Area}(\Omega^{*}_{m})\left(u_{\Omega^{*}}(p_{\Omega_{m}^{*}})-u_{\Omega^{T}}(p_{\Omega_{m}^{*}})\right)^{2}}{\sum_{m}\text{Area}(\Omega^{*}_{m})\left(u_{\Omega^{*}}(p_{\Omega_{m}^{*}})-u_{\Omega^{0}}(p_{\Omega_{m}^{*}})\right)^{2}}\text{.}$
(2)

Here, we use the squared error during evaluation as it is more sensitive to
outliers and thus provides a cleaner measure of the quality of the mesh. We
additionally provide results for the linear error, i.e., for Equation 2
without the squares, in Appendix C.3.

###  C.3 Linear Error Metric

Figure 5 shows results for the linear error. The general trends for the linear
errors are consistent with the quadratic errors in Figure 3. However, we find
that the linear metrics does not punish meshes with very high local errors as
much, resulting in slightly better results for methods that produce less
adaptive meshes.

(a) Laplace’s Equation

(b) Poisson’s Equation

Figure 5: Normalized linear error plotted against the number of elements for
(5(a)) the Laplace Equation and (5(b)) the Poisson Equation. For the linear
error, areas with a high error are weighted less. As a result, methods that
produce more regular meshes perform slightly better, because their main source
of error is not taken into account as much.

###  C.4 Additional Ablations

(a) Laplace’s Equation

(b) Poisson’s Equation

Figure 6: Normalized squared error plotted against the number of elements for
ablations of ASMR on (6(a)) the Laplace Equation and (6(b)) the Poisson
Equation. Finding a good refinement for Laplace’s Equation seems to be
possible without either global information or information about the solution
$u(x)$. However, for Poisson’s Equation, the performance degrades when not
providing global features. It further decreases for policies without any
access to global messages and when not providing solution information per mesh
element.

We experiment with a version of ASMR that does not use global features
$\mathbf{g}$ (ASMR (No Global Features)) but still relies on a global update
to pass messages across the full mesh, and with a version that removes the
global message passing altogether (ASMR (No Global Messages)). In both cases,
the global information of the mesh is limited, and for the latter, the
policies have to act purely from local information. Additionally, we consider
an ablation of our approach that does not use the solution information $u(x)$
as a node feature. For this ablation, the RL algorithm has to refine relevant
areas of the domain from just an encoding of the mesh and basic information
about the underlying PDE.

Figure 6 shows the results for these ablations. We find that all methods
perform roughly equally well on the simpler Laplace Equation. However, the
solution information and global features and messages significantly improve
performance on the more challenging Poisson problem. Interestingly, the
results for ASMR (No Solution) suggest that the RL algorithm is able to
calculate relevant regions of the PDE from just an encoding of the domain and
the boundary conditions. This is likely achieved by implicitly construction
some form of solution within the policy.

##  Appendix D Hyperparameters

###  D.1 General Hyperparameters

We use the same general hyperparameters across all methods and environments
for the sake of simplicity. All experiments are repeated for $n=10$ random
seeds with randomized PDEs and network parameters. During training, a new PDE
is sampled for every environment reset. All PDEs are normalized to be in
$(0,1)^{2}$. We evaluate all approaches for the final policies on $100$
randomly sampled PDEs that we keep consistent across random seeds for better
comparability.

For PPO, we largely follow the suggestions of Andrychowicz et al. (2021). We
train for a total of $500$ iterations, drawing $256$ samples per iteration and
then training for $5$ epochs with a batch size of $32$. All networks are
implemented in PyTorch (Paszke et al., 2019) and trained using the (Kingma &
Ba, 2014) ADAM optimizer with a learning rate of $3.0e-4$. The loss of the
value function is multiplied with a factor of $0.5$. We clip the norm of the
gradient during backpropagation to $0.5$, and set the clip range for the
policy and the value function to $0.2$. We normalize the observations with a
running mean and standard deviation, use a discount factor of $\gamma=0.99$
and estimate the advantages using Generalized Advantage Estimate (Schulman et
al., 2015) with $\lambda=0.95$. The policy and value function heads MLPs with
$2$ hidden layers with a latent dimension of $32$ and use a tanh activation
function.

We use separate MPNs for the policy and the value function. Each MPN consists
of $2$ Message Passing Steps, where each function $f_{\cdot}$ is represented
as a $2$-layer MLP with a latent dimension of $32$. All message aggregations
$\bigoplus$ are mean aggregations. Additionally, we apply Layer Normalization
(Ba et al., 2016) and Residual Connections (He et al., 2016) independently for
the node, edge and global features after each message passing step.

Table 1: Ranges for the refinement hyperparameters for the Laplace Equation. The Argmax baseline considers a different number of timesteps $T$,while the error-based Heuristic uses different thresholds $\theta$. All other methods make use of an element penalty $\alpha$. Note that values of $\alpha$ need to be different when scaling/not scaling the reward in Equation 1 with the element area term. Method | Parameter |
---|---|---
ASMR | $\alpha$ | $\\{3e-1,1e-1,3e-2,1e-2\\}$
ASMR (Shared Reward) | $\alpha$ | $\\{3e-1,1e-1,3e-2,1e-2\\}$
ASMR (Unscaled) | $\alpha$ | $\\{1e-5,3e-5,1e-4,3e-4,1e-3,3e-3,1e-2\\}$
PPO-VDGN | $\alpha$ | $\\{1e-5,3e-5,1e-4,3e-4,1e-3,3e-3,1e-2\\}$
Argmax | $T$ | $\\{20,60,100\\}$
Heuristic | $\theta$ | $\\{0,0.02,0.04,...,0.98\\}$
Table 2: Ranges for the refinement hyperparameters for the Poisson Equation. The Argmax baseline considers a different number of timesteps $T$, while the error-based Heuristic uses different thresholds $\theta$. All other methods make use of an element penalty $\alpha$. Note that values of $\alpha$ need to be different when scaling/not scaling the reward in Equation 1 with the element area term. Method | Parameter |
---|---|---
ASMR | $\alpha$ | $\\{2e-1,1e-1,5e-2,2e-2,1e-2,5e-3\\}$
ASMR (Shared Reward) | $\alpha$ | $\\{2e-1,1e-1,5e-2,2e-2,1e-2,5e-3\\}$
ASMR (Unscaled) | $\alpha$ | $\\{1e-5,3e-5,1e-4,3e-4,1e-3,3e-3,1e-2\\}$
PPO-VDGN | $\alpha$ | $\\{1e-5,3e-5,1e-4,3e-4,1e-3,3e-3,1e-2\\}$
Argmax | $T$ | $\\{20,60,100\\}$
Heuristic | $\theta$ | $\\{0,0.02,0.04,...,0.98\\}$

###  D.2 Method-Specific Parameters

All methods control the number of mesh elements that they produce with a
parameter. The error-based Heuristic uses the refinement threshold $\theta$,
the Argmax baseline a number of timesteps $T$ to refine for, and all other
methods make use of an element penalty $\alpha$. Table 1 lists the different
ranges for these parameters for our experiments for Laplace’s Equation, and
Table 2 for Poisson’s Equation.

##  Appendix E Additional Visualizations

We provide additional visualizations for both systems of equations on randomly
selected evaluation environments and algorithm random seeds. For Laplace’s
Equation, visualizations are shown in Figure 7(ASMR), Figure 8 (PPO-VDGN),
Figure 9 (Argmax), and Figure 10 (error-based Heuristic). For Poissons’s
Equation, visualizations are provided in Figure 11(ASMR), Figure 12 (PPO-
VDGN), Figure 13 (Argmax), and Figure 14 (error-based Heuristic). The
visualizations show that ASMR is able to provide consistent refinements of
high quality for different mesh resolutions. Opposed to this, VDGN sometimes
over- or underrefines meshes for the same element penalty $\alpha$. The Argmax
baseline generally performs well, but sometimes focuses on uninteresting
regions of the mesh. Finally, the error-based Heuristic greedily refines the
elements with the largest errors. This may lead to issues for globally
propagating errors (Strauss, 2007).

Figure 7: Refinements for Laplace’s Equation for refinements of ASMR on
randomly sampled PDEs for different values of the element penalty $\alpha$.
Figure 8: Refinements for Laplace’s Equation for refinements of PPO-VDGN on
randomly sampled PDEs for different values of the element penalty $\alpha$.
Figure 9: Refinements for Laplace’s Equation for refinements of the Argmax
baseline on randomly sampled PDEs for different environment rollout lengths
$T$. Figure 10: Refinements for Laplace’s Equation for refinements of the
error-based Heuristic on randomly sampled PDEs for different refinement
thresholds $\theta$. Figure 11: Refinements for Poisson’s Equation for
refinements of ASMR on randomly sampled PDEs for different values of the
element penalty $\alpha$ Figure 12: Refinements for Poisson’s Equation for
refinements of PPO-VDGN on randomly sampled PDEs for different values of the
element penalty $\alpha$ Figure 13: Refinements for Poisson’s Equation for
refinements of the Argmax baseline on randomly sampled PDEs for different
environment rollout lengths $T$. Figure 14: Refinements for Poisson’s Equation
for refinements of the error-based Heuristic on randomly sampled PDEs for
different refinement thresholds $\theta$.