# MTAC: Hierarchical Reinforcement Learning-based Multi-gait Terrain-adaptive
Quadruped Controller

Nishaant Shah, Kshitij Tiwari, and Aniket Bera The authors are with the
Department of Computer Science, Purdue University, USA,
<EMAIL_ADDRESS>
###### Abstract

Urban search and rescue missions require rapid first response to minimize loss
of life and damage. Often, such efforts are assisted by humanitarian robots
which need to handle dynamic operational conditions such as uneven and rough
terrains, especially during mass casualty incidents like an earthquake.
Quadruped robots, owing to their versatile design, have the potential to
assist in such scenarios. However, control of quadruped robots in dynamic and
rough terrain environments is a challenging problem due to the many degrees of
freedom of these robots. Current locomotion controllers for quadrupeds are
limited in their ability to produce multiple adaptive gaits, solve tasks in a
time and resource-efficient manner, and require tedious training and manual
tuning procedures. To address these challenges, we propose MTAC: a multi-gait
terrain-adaptive controller, which utilizes a Hierarchical reinforcement
learning (HRL) approach while being time and memory-efficient. We show that
our proposed method scales well to a diverse range of environments with
similar compute times as state-of-the-art methods. Our method showed greater
than 75% on most tasks, outperforming previous work on the majority of test
cases.

##  I Introduction

The past few decades have witnessed a high number of natural and man-made
disasters. In the wake of these natural disasters and mass casualty incidents
such as earthquakes or building collapses, hazardous conditions make these
environments extremely difficult to access for first responders. Moreover, the
potential presence of asbestos, dust, extreme temperatures, and unstable
structures endanger the lives of first responders. In the search for advanced
technologies to alleviate these problems, mobile robots have been increasingly
discussed as potential solutions. Robots have been used to assist first
responders in recent devastation such as the collapse of the World Trade
Center, the 2004 Mid-Niigata earthquake in Japan, the 2005 Hurricanes Katrina,
Rita, and Wilma in the United States, as well as the 2011 Tohoku earthquake
and tsunami in Japan [1]. Nevertheless, the challenge of robot-assisted urban
search and rescue is far from being solved. In order to fully explore post-MCI
environments, robots must be able to navigate dynamic obstacles, as well as
cluttered and unstructured terrains. In recent years, tracked robots have been
used to drive across uneven terrain. For, in response to the Great Eastern
Japan Earthquake, tracked robots such as the KOHGA3 were used for inspection
of disaster sites [2]. Although this teleoperated robot was able to drive over
uneven ground and use its sensors to provide the operators with information
about the site, it was limited in its ability to climb stairs or steep slopes.

Figure 1: The proposed controller solves rough terrain navigation using
hierarchical learning. A high-level policy is trained to act over a family of
pre-trained low-level experts trained to execute unique gaits. Together, these
policies are able to execute adaptive locomotion on unstructured terrains

Quadruped robots, thanks to the high degree of versatility that comes with
their legged design, have great potential to navigate complex environments
that their counterparts, such as wheeled or tracked robots, cannot [3, 4, 5].
Because quadrupeds use legs to move around, they have the freedom to choose
the points at which they make contact with their surroundings and quickly
change their footholds based on their environment. These features give
quadruped robots the potential to perform agile locomotion that previously
only animals were capable of [6]. Despite their highly versatile design,
however, the control of these legged systems has proven to be a challenging
task. Although their physical design is bio-inspired, existing control
algorithms for quadrupeds struggle to replicate animals’ natural ability for
dynamic locomotion over diverse terrains. Locomotion control for quadrupeds in
rough terrain is a challenging problem, because it requires the robot to adapt
its gait to diverse conditions, and the resources on these mobile robots are
constrained [7, 8]. Currently, both conventional model-based controllers, as
well as learning-based model-free controllers, have shown some success in the
field of quadrupedal locomotion for rough terrain. Nevertheless, they have
limitations when it comes to addressing a diverse set of environmental
conditions and resource efficiency. Furthermore, existing controllers
typically converge to a standard trotting gait, which is not always optimal
for certain velocity ranges or all terrains [6, 9]. In fact, specific gaits
are more optimal for certain velocity ranges than others [10]. Multiple gait
modes allow the robot to switch between gaits, thereby performing more complex
and dynamic motions, such as changing from walking to pronking [11, 12]. It
has been shown that behaviors learned from scratch are more energy efficient
for locomotion tasks and capable of more complex locomotion compared to
primitive gaits like trot [13]. Together, these works demonstrate the need for
a learning-based controller that can illicit multiple emergent gaits from the
robot, while accomplishing tasks in a time and energy-efficient manner. In
this paper, we propose the use of hierarchical reinforcement learning (HRL) to
train a multi-gait terrain-adaptive locomotion controller for quadruped robots
(see Fig. 1). To this end, we make the following novel contributions:

  * •

The policies are trained end-to-end, without using human-engineered trajectory
generators (TG) or motion primitives. We train our locomotion controller from
scratch, using the design of the terrain in the simulated training curriculum
to elicit different gaits rather than developing gaits from motion primitives.
The resulting gaits are more natural and diverse because they do not draw
inspiration from existing motions. These gaits are referred to as emergent
behaviors, and they are specialized for the conditions in which they are
trained. This yields more efficient performance than motion primitives that
are generalized to multiple terrain types. In this work, we focus on training
environments that contain stairs, gaps, and uneven ground. These three
distinct environments elicit unique behaviors from their respective policies.
Training the high-level policy to modulate over these pre-trained policies
allows for these expert skills to be synthesized within a single navigation
task/mission.

  * •

The model reduces navigation time by generating multiple emergent gaits and
performing at high velocities. The methodology allows for a reduction of
navigation time, by training a family of complimentary and unique gaits.
Previous works have shown that certain gaits are more efficient and
appropriate for certain velocity ranges, and terrain types. By straying from
traditional pre-determined gaits and eliciting emergent behaviors from our
low-level expert policies, our proposed controller reduces the time taken to
complete a navigation task/mission. In doing so, we extend the potential range
of the robot for a single mission, without the need for external batteries,
which would be a cumbersome addition to the quadruped’s payload. This is
necessary for applications in urban search and rescue, such as earthquakes,
when missions must be completed quickly- especially due to the high chance of
after-shocks which could further compromise the integrity of buildings and
human lives.

Figure 2: Shown above is MTAC’s pre-trained stair expert successfully
navigating a difficult grade stair pyramid, as a model trained through the
previous work’s generalized curriculum falls on a lower difficulty level.

##  II Related Works

This section discusses the varied control methods adopted for quadrupedal
locomotion, focusing on model-based controls, impact-aware locomotion, and
learning-based methodologies.

###  II-A Quadrupedal Locomotion

Quadrupeds, with their myriad degrees of freedom, pose significant control
challenges, especially when aiming for agile motion. The literature reveals a
range of strategies that address this intricacy.

####  II-A1 Model-based Control

Traditional approaches in this category leverage motion planning algorithms,
deeply rooted in the robot’s dynamics [14, 12, 15]. These approaches usually
employ a blend of state machines and meticulously crafted foot trajectories.
They consider various parameters, including potential slippage and ensuring
appropriate ground contact [16]. While these methods demonstrate
effectiveness, they sometimes falter in dynamic scenarios where real-time
adaptability is paramount. Notably, Fankhauser et al. introduced a more
adaptive controller. It continuously maps the robot’s environment, adjusting
to changes such as new obstacles. However, it does have limitations in
recognizing deformable surfaces or accounting for terrain alterations as the
robot moves [17].

####  II-A2 Impact-aware control

The evolution of robotic hardware has sparked interest in controllers that
integrate direct feedback from the robot’s appendages. Guo et al.’s work,
which introduced a sensorized foot pad capable of terrain classification,
exemplifies this shift [18]. By merging such sensor feedback with existing
control algorithms, there’s a marked improvement in how robots adapt to
terrains [19, 20]. A noteworthy development in this domain is Zhang et al.’s
risk assessment network, which employs a combination of elevation maps and
learning policies to navigate rough terrains. While promising, this model
acknowledges some limitations in teaching agile behaviors [21].

####  II-A3 Learning-based control

The increasing dominance of reinforcement learning in robotic control is
evident in quadrupedal locomotion as well. Policies Modulating Trajectory
Generators (PMTG) is a pioneering architecture in this space, where neural
networks determine trajectory parameters based on robot states [22]. The
Evolved Environmental Trajectory Generators (EETG) framework builds upon PMTG.
It introduces a level of diversity in trajectory generation, even if its
primary results are limited to simulations [11, 23]. Real-world
implementations of PMTG-like architectures have seen success, especially Tan
et al.’s work that exhibited robustness across varied terrains [9]. In
contrast, Lee et al. emphasized training under diverse simulated conditions,
though their results showed a preference for a single gait [24]. Finally,
Rudin et al. took a different track, advocating for a holistic end-to-end
approach. This method demonstrated unique gaits and better energy efficiency,
but it had its set of challenges, especially in reverse locomotion [13, 25].

##  III Multi-gait Terrain-adaptive Controller (MTAC)

In this section, we first present an overview of the proposed multi-gait
controller MTAC (Sec. III-A) followed by detailed descriptions of the
corresponding high and low-level controllers in Sec. III-B and Sec. III-C,
respectively.

Figure 3: The methodology for MTAC is pictured above. At the highest level, is
the input of a high-level action into the controller. Based on this command,
the high-level policy then selects a low-level action to be executed by one of
the 3 expert policies. The training process for these experts is also
pictured. Then the robot executes the low-level action and interacts with the
environment. These contacts then start the next iteration of the feedback
loop.

while controller running do
Get command velocity
Get observation $o^{high}_{t}$
$a^{high}_{t}$, $d$ = $\pi^{high}(a_{t}|o_{t})$
while $t_{elapsed}<d$
Get observation $o^{low}_{t}$
Run $\pi^{low}$

Algorithm 1 Locomotion Control Algorithm

###  III-A Overview

The proposed controller splits the task of locomotion into multiple sub-tasks
using a hierarchical learning architecture. An overview of the method is shown
in Fig. 3. The Multi-gait Terrain-Adaptive Controller
(MTAC)111https://ideas.cs.purdue.edu/research/QuadrupedLocomotion consists of
a high-level control policy, which regulates a family of low-level policies
trained to exhibit diverse emergent skills. These skills are optimized for a
range of sub-tasks. Each low-level policy becomes an expert in certain skills
that are optimized for its respective sub-task. The high-level policy selects
which expert to apply to unseen environments with the goal of successfully
navigating to a desired goal. This method is tested in cluttered and
unstructured environments in simulation. The execution of the controller is
described in Algorithm 1. First, the controller receives a commanded velocity
vector for the robot’s base. Then an observation of the robot’s current state
is made. Based on this observation, the high-level policy outputs a high-level
action and the estimated duration of this action. The high-level action
activates a low-level gait mode, which controls the motion of the legs of the
robot as it traverses the terrain.

Load environment $env$
Send command $c$ or $t_{elasped}$ $<$ 25 do
Get observations $o^{high}_{t}$, $o^{low}_{t}$
$a^{high}_{t}$= $\pi^{high}(a_{t}|o_{t})$
$\pi_{expert}$ = $a^{high}_{t}$
$a^{low}_{t}$ = $\pi_{expert}(a_{t}|o_{t})$
$env$.increment-step($a^{low}_{t}$)
Calculate reward $R(t)$
Pass reward to $\pi_{high}$

Algorithm 2 High Level Controller Training

Load environment $env$
while $t_{elapsed}$ $<$ 20 or $n_{steps}$ $<$ 25 do
Get observation $o^{low}_{t}$
$a^{low}_{t}$= $\pi^{low}(a_{t}|o_{t})$
$env$.increment-step($a^{low}_{t}$)
Calculate reward $R(t)$
Pass reward to $\pi_{low}$

Algorithm 3 Low Level Controller Training

###  III-B High Level Controller

One of the novel components of our approach, MTAC, is the high-level
controller which acts as the policy over policies. Instead of a single policy
controlling locomotion, the high-level controller is trained using DRL to
allocate certain periods of locomotion to the more specialized expert
policies. We represent the problem of locomotion control in discrete time
steps. At each time step $t$, the robot obtains an observation $o_{t}\in O$,
performs an action $a_{t}\in A$, and receives a reward $r_{t}\in R$. The robot
chooses $a_{t}$ based on a policy $\pi(a_{t}|O_{t})$. We use deep
reinforcement learning (DRL) to solve this control problem. Using DRL, the
goal is to obtain a policy $\pi^{*}$ that maximizes the sum of the discounted
rewards over an infinite time horizon. In order to implement the DRL
framework, we first formulate the high-level task of locomotion as a Markov
decision process (MDP). MDP is a mathematical tool used in reinforcement
learning for modeling discrete-time decision-making problems that are
partially controllable and partially stochastic. An MDP is defined by the
tuple $(A,S,p,r)$, where $A\in\mathbb{R}^{3}$ represents the continuous action
space, $S$ where $S\in\mathbb{R}^{15}$ represents the continuous state space,
$p$ is the unknown state-transition probability, and the reward is $r$. This
tuple is stored in memory in a replay buffer so that the policy is able to
keep track of its history. For the high-level controller, we define a
continuous action space consisting of a confidence rating for the performance
of each pre-trained expert policy. The expert with the highest confidence is
then selected to control the robot’s gait. For this controller, the state
space or observation consists of the quadruped’s 12 joint positions and
velocities, as well as the base twist, orientation, and direction command. The
training procedure for the high-level controller is presented in Algorithm 2.

###  III-C Low Level Control

For the low-level controllers, we train a family of expert policies. Each
expert policy consists of a neural network trained using the proximal policy
optimization (PPO) algorithm [25]. By using PPO, we take advantage of entropy
regularization to encourage the robot to explore a wider range of behaviors or
take actions with more entropy. As done with the high-level controller, we
also formulate the problem of low-level control as an MDP. We define the
continuous action space$A\in\mathbb{R}^{D}$ as a set of joint reference angles
for each of the robot’s actuated joints, D being the number of degrees of
freedom of the particular quadruped. These desired joint states are then fed
into the robot’s joint PD controller. The controller’s observation space is
defined by the robot’s joint positions and velocities as well as base twist,
orientation, and direction command. For training the experts, we expand on the
massively parallel training procedure presented in [26] as this work is open
source and allows for rapid training of legged locomotion policies, by
simulating thousands of robots in a single environment and learning a policy
in a parallel manner. During training, domain randomization is applied in the
form of friction randomization, mass randomization, and randomized
disturbances to the robot’s base.

Figure 4: Expert policies are developed by training in specialized
environments that consist of focused terrain types. Here we show three
environments: bumpy terrain (top left), stair pits (top right), stair pyramids
(bottom left) and steps (bottom right). All of these terrain types are
parameterized so that difficulty can be gradually varied across the map.

We train three expert policies using a specialized terrain curriculum with
increasing difficulty at each level within the map. By doing so, we allow each
policy to learn a subset of skills that is more specialized towards its
respective environment. Each simulation environment is shown in Fig. 4, and
described in detail later in this paper. Each policy is trained in its
respective environment with the same process. First, A 10 by 10 grid is
generated consisting of the environment’s respective terrain. The grid is
arranged so that the terrains of varying difficulty are presented. In each
row, the parameters of the terrain are incremented in order to make the
terrain more difficult to traverse. Then the robots are placed in the
simulation in large clusters in the middle of each grid square. Short episodes
of 20 seconds are implemented, and an always positive rewards scale is
implemented in order to avoid early termination. At the end of each episode,
agents are moved to other points on the grid based on performance, so that
over the course of the training period, all of the agents will first learn
from the easy examples of the terrain, before learning to address more severe
and difficult terrains. One training iteration contains 25 steps, and 1500
iterations of this process are performed. After 1500 iterations, a final
policy is generated. In training, the maximum episode length is 20 seconds,
and episodes last over multiple iterations and, therefore, multiple policy
updates. The low-level training procedure is detailed in Algorithm 3. As
opposed to the previous work, which aims to improve performance by approaching
the upper limit of robots simulated, we reduce the number of robots per
environment from 4096 to 1024, allowing for training to be more lightweight,
while still converging on more specialized behaviors. Furthermore, by training
each expert in its own environment, we increase the ceiling for the maximum
difficulty of terrain faced by the agents in training. We also focus more
robots on a single locomotion task at a time, rather than splitting their
focus across learning multiple types of terrains.

####  III-C1 Bumpy Terrain

In this environment (shown in Fig. 4 (top left)), the quadruped is trained to
follow a goal velocity while traversing a terrain with many random bumps and
elevation changes. In this environment, the policy is learned in which the
robot can quickly and accurately follow its commanded velocity while
maintaining balance and safe footholds. The minimum and maximum height of the
irregularities parameterize this terrain. As the difficulty scalar increases,
the range of the min and max heights does as well. Domain randomization is
applied, consisting of random force applied to the base as well as variations
in the mass of the payload of the robot. Here, the quadruped learns a unique
gait similar to the trot; however, it includes dynamic movements in the
shoulder joints, which allow for more versatile motion in multiple directions.
The robot also learns how to walk both forward and backward, as well as a
sideways shuffling motion. Unlike the policy trained in stairs, this policy
exhibits more economic steps, in which the robot lifts its feet just high
enough to step over the irregularities in the ground, without lifting the foot
too high and spending unnecessary time and energy. It can also be observed
that this policy executes similar efficient motions when moving laterally as
well, by not over-actuating its shoulder joints. The agents also show smooth
motion along non-uniform trajectories, such as a smooth transition from
lateral motion to forward motion as the direction of the goal velocity
changes.

####  III-C2 Stair pits and pyramids

In this environment (shown in Fig. 4 (top right and bottom left,
respectively)), the quadruped is trained to follow a goal velocity while
climbing up, down, and over pyramids of stairs with increasing difficulty as
they move around the map. Both the stair pyramids and stair pits are
parameterized in terms of step width, step height, and the size of the
platform at the apex. As the robot advances in the training curriculum, it
faces pyramids and pits with larger step heights and smaller platform sizes.
This gradual introduction of more difficult conditions allows the policy to
gradually learn how to locomote over the terrain, allowing for more
exploration of the solution space rather than an early introduction of
failure. In this environment, the robot learns a foot-trapping motion to lift
its legs in order to climb over the steps. Additionally, the agents exhibit
more independent leg movements on the way down the staircases, as well as
sideways movements down the stairs, to maintain better balance. Additionally,
the robot replicates the foot-trapping motion in the horizontal axis as well,
using its shoulder joint to lift its legs while climbing up the stairs
sideways. By being able to move and traverse the terrain in multiple
directions, this policy saves navigation time and efficiency, unlike hand-
engineered trajectory generators, which do not include backward and sideways
motions. These motions are learned with respect to the terrain in which the
agent is trained, so unique behaviors such as foot trapping, and horizontal
foot trapping are specialized for the task of navigating stairs. When the same
model is run in a more general terrain curriculum, these behaviors are lost as
the agent tries to converge on a gait that is more applicable to the
generalized terrain. Fig.2 shows how these unique behaviors increase
performance.

####  III-C3 Steps

In this environment (shown in Fig. 4 (bottom right)), the quadruped is trained
to follow a goal velocity while traversing a terrain with randomly placed
rectangular steps on the ground. This terrain is parameterized by the height,
length, and width of the steps. As the difficulty scalar increases, the
dimensions of the steps increase. In this environment, the robot learns to
negotiate the terrain in a similar manner to the stair-climbing expert;
however, it adapts more to the fluctuating ups and downs. The trained policy
in this scenario exhibited a pronking-like behavior, however, with more
caution and without the airborne motion typically associated with pronking.
Instead, the model learns to synchronize its front and back legs to overcome
the height barrier of the presented steps.

##  IV Empirical Evaluation

###  IV-A Quadruped robot platform

The ANYmal quadruped is a 12-degree-of-freedom robot, with a torque-controlled
actuator on each joint [4]. When training and evaluating our policies, an
accurate physical model of the ANYmal C robot is used in simulation, together
with pre-trained actuator networks, to realistically mirror the real-world
actuation of ANYmal’s joints.

###  IV-B Scenario

We evaluated each of the pre-trained expert policies in environments
consisting of their respective terrains and compared the results with those of
a policy trained using the same method as [26]. Each environment consisted of
a single row of 5 sections of terrains. For the stair expert, these 5 sections
consisted of stairs of varying difficulty. For the bumpy ground expert, the
environment consisted of hills of various heights with varying difficulty in
bumps. For the step expert, step dimensions of various heights were included.
The expert test cases can be seen in Fig. 4. In each test, the robot is given
a forward velocity command, and the completion rate is measured over 10-15
trials. For each test scenario, the robot’s performance is evaluated on the
terrain of low and high difficulty, as well as with low and high commanded
velocity. For low difficulty, the difficulty scalar was reduced by 50%, this
scalar affects the height and severity of the changes in the terrain. For high
difficulty, the difficulty scalar was increased to 100%, while the exact
values vary for each task. For high velocity, the robots were commanded to
move forward at 1.75 meters per second, and for low velocity, they were
commanded to move forward at .75 meters per second.

###  IV-C Results

The results of these experiments are shown for stairs in Table I, for bumpy
terrain in Table II, and for stepped terrain in Table III. In each table the
completion rate is represented by the column C.R. After conducting multiple
trials, the number of trials in which the quadruped successfully makes it
across the terrain is recorded and represented as the Task Completion Rate.
Furthermore, the accuracy at which the expert policies are able to track their
speed and direction commands as well as the error between the desired joint
positions robot’s legs and the actual measured values. These results are shown
for the bumpy terrain expert in Fig. 6. Here we can see how this policy
converges to closely match the velocity commands as well as closely resembles
the targeted joint positions. Similar results are shown for the stair-climbing
expert in Fig. 7. However, for this policy, the measured velocities show more
oscillation and above and below the target curve due to the more drastic
changes in slope and elevation for the robot base. The success rate of the
policy shows that this is not an issue in terms of performance and rather
caused by the nature of the task. As shown in I, the MTAC stair expert showed
the best results with a high success rate. Similarly, the MTAC step expert
also had a high percentage of completion. For the bumpy terrain test, both
models performed similarly, performing well on the easier difficulties, and
with lower completion rates as the speed and difficulty increased.
Nevertheless, in the most difficult test case, with full difficulty for the
bumpy terrain and high velocity, MTAC did perform better, with a completion
rate more than double that of the previous work. These results show how gaits
generated using MTAC’s methodology result in greater adaptability to
environments with difficulties that were not faced in training, and with
navigating difficult terrain at velocity ranges higher than those faced in
training. In the case of an MCI, these traits would be more desirable as fast
response, and navigation of unpredictable terrain are both necessary in MCI
scenarios.

Figure 5: Here is shown the performance of the MTAC stepping expert during training. It shows how the body velocities converge closely to the commanded velocities. Similarly, the DOF positions also converge to match the desired positions. TABLE I: Stair pyramid performance metrics. Each model was tested under two types of difficulty scales and at high and low velocity, with many trials for each scenario. Policy | Terr. Type | Terr. Difficulty | Vel. (m/s) | C.R.
---|---|---|---|---
MTAC | Stairs | 50% | 0.75 | 100%
Generalized PPO | Stairs | 50% | 0.75 | 50%
MTAC | Stairs | 100% | 0.75 | 100%
Generalized PPO | Stairs | 100% | 0.75 | 30%
MTAC | Stairs | 50% | 1.75 | 100%
Generalized PPO | Stairs | 50% | 1.75 | 25%
MTAC | Stairs | 100% | 1.75 | 82%
Generalized PPO | Stairs | 100% | 1.75 | 13%
TABLE II: Bumpy terrain performance metrics. Each model was tested under two types of difficulty scales and at high and low velocity, with many trials for each scenario. Policy | Terr. Type | Terr. Difficulty | Vel. (m/s) | C.R.
---|---|---|---|---
MTAC | Bumpy | 50% | 0.75 | 100%
Generalized PPO | Bumpy | 50% | 0.75 | 100%
MTAC | Bumpy | 100% | 0.75 | 92%
Generalized PPO | Bumpy | 100% | 0.75 | 91%
MTAC | Bumpy | 50% | 1.75 m/s | 75%
Generalized PPO | Bumpy | 50% | 1.75 m/s | 79%
MTAC | Bumpy | 100% | 1.75 m/s | 34%
Generalized PPO | Bumpy | 100% | 1.75 m/s | 13%
TABLE III: Stepped terrain performance metrics. Each model was tested under two types of difficulty scales and at high and low velocity, with many trials for each scenario. Policy | Terr. Type | Terr. Difficulty | Vel. (m/s) | C.R.
---|---|---|---|---
MTAC | Stepped | 50% | 0.75 | 83%
Generalized PPO | Stepped | 50% | 0.75 | 41%
MTAC | Stepped | 100% | 0.75 | 61%
Generalized PPO | Stepped | 100% | 0.75 | 33%
MTAC | Stepped | 50% | 1.75 m/s | 75%
Generalized PPO | Stepped | 50% | 1.75 m/s | 42%
MTAC | Stepped | 100% | 1.75 m/s | 58%
Generalized PPO | Stepped | 100% | 1.75 m/s | 34%
Figure 6: Metrics from training in the bumpy terrain environment are shown. We
can see that the x and y base velocities, as well as the yaw of the robot,
converge towards the commanded values as training time progresses. Similarly,
as time progresses, the positions of the robot’s degrees of freedom also more
closely resemble the commanded values. Figure 7: Metrics from training in the
stair pyramid terrain environment are shown. We can see that the x and y base
velocities, as well as the yaw of the robot, converge towards the commanded
values and oscillate around the endpoint as training time progresses.
Similarly, as time progresses, the reference positions of the robot’s degrees
of freedom also closely resemble the commanded values.

##  V Conclusion

We presented a hierarchical reinforcement learning-based controller, MTAC, for
rough terrain locomotion in urban search and rescue scenarios for quadruped
robots. Without additional sensors or previous knowledge of the simulated
terrain, the proposed controller can learn unique and terrain-specialized
motions and behaviors to increase the speed and efficiency at which it
traverses over different types of terrains and irregularities. The controller
is evaluated in multiple test scenes. First, each pre-trained expert is
compared against the generalized controller in its native environment type,
however, with randomized terrain parameters. Next, the integrated MTAC
controller is tested against the previous work in a test map consisting of
randomly generated terrains with a goal destination at the end. Each model is
evaluated on the time taken to complete its given task, the success rate at
which it completes tasks, and the average error between its actual and command
base velocity over the time it completes the given task. In this methodology,
basic domain randomization is used to allow for sim-to-real transfer. However,
in this work, we do not focus on transferring the policies learned in
simulation to a real robot. In the future, more focus should be applied to
improving the existing sim-to-real strategies and applying them to this model.
Furthermore, future work should include the use of a gating neural network
(GNN) to synthesize the weights of the low-level policies and derive a new
policy at each time step. This method would increase the adaptive nature of
the generated gaits, and create more natural and complex behaviors. Lastly,
future work should explore the use of non-markovian reward scales during
training in order to encourage or elicit specific natural behaviors and gaits.

## References

  * [1] Y. Liu, G. Nejat, and B. Doroodgar, “Learning based semi-autonomous control for robots in urban search and rescue,” in _2012 IEEE International Symposium on Safety, Security, and Rescue Robotics (SSRR)_ , 2012, pp. 1–6.
  * [2] F. Matsuno, N. Sato, K. Kon, H. Igarashi, T. Kimura, and R. Murphy, _Utilization of Robot Systems in Disaster Sites of the Great Eastern Japan Earthquake_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2014, pp. 1–17.
  * [3] G. Bledt, M. J. Powell, B. Katz, J. Di Carlo, P. M. Wensing, and S. Kim, “Mit cheetah 3: Design and control of a robust, dynamic quadruped robot,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 2245–2252.
  * [4] M. Hutter, C. Gehring, D. Jud, A. Lauber, C. D. Bellicoso, V. Tsounis, J. Hwangbo, K. Bodie, P. Fankhauser, M. Bloesch, R. Diethelm, S. Bachmann, A. Melzer, and M. Hoepflinger, “Anymal - a highly mobile and dynamic quadrupedal robot,” in _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2016, pp. 38–44.
  * [5] B. Katz, J. D. Carlo, and S. Kim, “Mini cheetah: A platform for pushing the limits of dynamic quadruped control,” in _2019 International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 6295–6301.
  * [6] J. Hwangbo, J. Lee, A. Dosovitskiy, D. Bellicoso, V. Tsounis, V. Koltun, and M. Hutter, “Learning agile and dynamic motor skills for legged robots,” _Science Robotics_ , vol. 4, no. 26, p. eaau5872, 2019. [Online]. Available: https://www.science.org/doi/abs/10.1126/scirobotics.aau5872
  * [7] K. Tiwari, X. Xiao, and N. Y. Chong, “Estimating achievable range of ground robots operating on single battery discharge for operational efficacy amelioration,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 3991–3998.
  * [8] K. Tiwari, X. Xiao, A. Malik, and N. Y. Chong, “A unified framework for operational range estimation of mobile robots operating on a single discharge to avoid complete immobilization,” _Mechatronics_ , vol. 57, pp. 173–187, 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0957415818301946
  * [9] W. Tan, X. Fang, W. Zhang, R. Song, T. Chen, Y. Zheng, and Y. Li, “A hierarchical framework for quadruped locomotion based on reinforcement learning,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2021, pp. 8462–8468.
  * [10] Y. Kim, B. Son, and D. Lee, “Learning multiple gaits of quadruped robot using hierarchical reinforcement learning,” _CoRR_ , vol. abs/2112.04741, 2021\. [Online]. Available: https://arxiv.org/abs/2112.04741
  * [11] S. Surana, B. Lim, and A. Cully, “Efficient learning of locomotion skills through the discovery of diverse environmental trajectory generator priors,” in _2023 IEEE International Conference on Robotics and Automation (ICRA)_ , 2023, pp. 12 134–12 141.
  * [12] C. D. Bellicoso, F. Jenelten, C. Gehring, and M. Hutter, “Dynamic locomotion through online nonlinear motion optimization for quadrupedal robots,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 3, pp. 2261–2268, 2018\.
  * [13] N. Rudin, D. Hoeller, M. Bjelonic, and M. Hutter, “Advanced skills by learning locomotion and local navigation end-to-end,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022, pp. 2497–2503.
  * [14] M. Kalakrishnan, J. Buchli, P. Pastor, M. Mistry, and S. Schaal, “Fast, robust quadruped locomotion over challenging terrain,” in _2010 IEEE International Conference on Robotics and Automation_ , 2010, pp. 2665–2670.
  * [15] Y. H. Lee, Y. H. Lee, H. Lee, H. Kang, Y. B. Kim, J. H. Lee, L. T. Phan, S. Jin, H. Moon, J. C. Koo, and H. R. Choi, “Whole-body motion and landing force control for quadrupedal stair climbing,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2019, pp. 4746–4751.
  * [16] F. Jenelten, J. Hwangbo, F. Tresoldi, C. D. Bellicoso, and M. Hutter, “Dynamic locomotion on slippery ground,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 4170–4176, 2019.
  * [17] P. Fankhauser, M. Bjelonic, C. Dario Bellicoso, T. Miki, and M. Hutter, “Robust rough-terrain locomotion with a quadrupedal robot,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 5761–5768.
  * [18] X. Guo, B. Blaise, J. Molnar, J. Coholich, S. Padte, Y. Zhao, and F. L. Hammond, “Soft foot sensor design and terrain classification for dynamic legged locomotion,” in _2020 3rd IEEE International Conference on Soft Robotics (RoboSoft)_ , 2020, pp. 550–557.
  * [19] G. Valsecchi, R. Grandia, and M. Hutter, “Quadrupedal locomotion on uneven terrain with sensorized feet,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 1548–1555, 2020.
  * [20] C. Dario Bellicoso, F. Jenelten, P. Fankhauser, C. Gehring, J. Hwangbo, and M. Hutter, “Dynamic locomotion and whole-body control for quadrupedal robots,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2017, pp. 3359–3365.
  * [21] H. Zhang, J. Wang, Z. Wu, Y. Wang, and D. Wang, “Terrain-aware risk-assessment-network-aided deep reinforcement learning for quadrupedal locomotion in tough terrain,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2021, pp. 4538–4545.
  * [22] A. Iscen, K. Caluwaerts, J. Tan, T. Zhang, E. Coumans, V. Sindhwani, and V. Vanhoucke, “Policies modulating trajectory generators,” in _Proceedings of The 2nd Conference on Robot Learning_ , ser. Proceedings of Machine Learning Research, A. Billard, A. Dragan, J. Peters, and J. Morimoto, Eds., vol. 87. PMLR, 29–31 Oct 2018, pp. 916–926. [Online]. Available: https://proceedings.mlr.press/v87/iscen18a.html
  * [23] J. Mouret and J. Clune, “Illuminating search spaces by mapping elites,” _CoRR_ , vol. abs/1504.04909, 2015. [Online]. Available: http://arxiv.org/abs/1504.04909
  * [24] J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning quadrupedal locomotion over challenging terrain,” _Science Robotics_ , vol. 5, no. 47, p. eabc5986, 2020. [Online]. Available: https://www.science.org/doi/abs/10.1126/scirobotics.abc5986
  * [25] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _CoRR_ , vol. abs/1707.06347, 2017. [Online]. Available: http://arxiv.org/abs/1707.06347
  * [26] N. Rudin, D. Hoeller, P. Reist, and M. Hutter, “Learning to walk in minutes using massively parallel deep reinforcement learning,” _CoRR_ , vol. abs/2109.11978, 2021. [Online]. Available: https://arxiv.org/abs/2109.11978


# PARALLEL METRIC-BASED MESH ADAPTATION IN PETSC USING PARMMG

Joseph G. Wallwork1,∗ Matthew G. Knepley2 Nicolas Barral3,4 Matthew D.
Piggott1

( 1Department of Earth Science and Engineering, Imperial College London,
London, U.K.
2Computer Science and Engineering, State University of New York, Buffalo, NY,
U.S.A.
3University of Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, 33400, Talence,
France
4INRIA, IMB, UMR 5251, 33400, Talence, France
∗Corresponding author. E-mail<EMAIL_ADDRESS>
)

###### Abstract

This research note documents the integration of the MPI-parallel metric-based
mesh adaptation toolkit _ParMmg_ into the solver library _PETSc_. This
coupling brings robust, scalable anisotropic mesh adaptation to a wide
community of PETSc users, as well as users of downstream packages. We
demonstrate the new functionality via the solution of Poisson problems in
three dimensions, with both uniform and spatially-varying right-hand sides.

##  1 Introduction

Metric-based mesh adaptation is now a mature technology, which has the
capacity to produce highly anisotropic, multi-scale meshes of scientifically
and industrially relevant domains in two or three dimensions. Its
effectiveness as an advanced discretisation method has been demonstrated on
numerous occasions, in application areas ranging from geoscience to aerospace
engineering [1, 2]. The core idea is to use a Riemannian metric space within
the mesher, to guide the concept of mesh optimality. Whilst an element of an
adapted mesh might look distorted in Euclidean space, it is ‘unit’ when viewed
in the Riemannian space – its edges having (close to) unit length. A great
advantage of the metric-based approach is that it allows control of element
shape and orientation, as well as size. This can be particularly beneficial
for problems with strong direction-dependence and/or anisotropy.

PETSc (Portable, Extensible Toolkit for Scientific Computing) [3, 4] is a
widely used scientific library, written in C, which provides the linear and
nonlinear solvers and associated data structures required to solve partial
differential equations in parallel, among many other things. Metric-based mesh
adaptation is already available to PETSc’s unstructured mesh manager, DMPlex
[5, 6], using _Pragmatic_ , as documented in a previous IMR research note [7].
The software development efforts documented in this research note build upon
the previous work, integrating two new toolkits – Mmg (serial) [8] and ParMmg
(MPI-parallel) [9] – as well as providing an interface to DMPlex, with
routines for constructing, modifying and combining Riemannian metrics. Python
bindings are provided through _petsc4py_ [10], so that downstream packages,
such as the _Firedrake_ finite element library [11], may take advantage of
these enhancements.

The remainder of this research note is as follows. Section 2 describes the
implementation of the PETSc-ParMmg coupling, including information on how it
can be integrated into an existing workflow. Section 3 demonstrates its
application to Poisson problems with two different choices of RHS. Finally, an
outlook on future work is made in Section 4. An appendix is provided, which
includes the commands required to reproduce the experimental results.

##  2 Implementation

###  2.1 Mmg and ParMmg

_Mmg_ comprises a collection of tools for anisotropic metric-based adaptation
of simplicial meshes [8]. It includes a 2D local remesher (_Mmg2d_), a 3D
local remesher (_Mmg3d_) and a surface remesher (_Mmgs_). The 3D local
remesher has been parallelized using MPI as _ParMmg_ [9]. Parallel
implementations do not currently exist in 2D or for surface remeshing. Given
an input mesh and a metric defined upon it, Mmg applies a sequence of
operations to the mesh in order to optimize the quality of its elements, as
measured by a quality functional, for which optimality means a perfectly
isotropic element with edges of unit length. In the 3D case, the operations
used to optimize the mesh are ($h$-adaptive) node insertion, node deletion and
face swapping, as well as ($r$-adaptive) node movement.

ParMmg performs parallel mesh adaptation iteratively, first applying Mmg3d to
the mesh partition owned by each process, then repartitioning and doing the
same again. Three iterations are typically sufficient to remove dependence on
the initial partition.

###  2.2 PETSc DMPlex

In PETSc, unstructured meshes are managed by the DMPlex data structure. The
representation uses a Hasse diagram, or directed acyclic graph [5], whereby
all entities are treated equally as vertices of the graph. This enables
algorithms to be written independently of mesh properties such as dimension
and cell type [12]. DMPlex supports all of the operations required to handle
meshing, such as generation, partitioning and distribution, creating of
missing edges and faces, regular refinement, extrusion, traversal, selection,
manipulation and I/O [13]. This work adds Riemannian metric utilities, to
support metric-based adaptation.

###  2.3 Coupling

#### Overview

Given a linear simplicial mesh described using a DMPlex, a Riemannian metric
field may be defined in the finite element context as a piece-wise linear and
continuous ($\mathbb{P}1$) field, with its degrees of freedom (DoFs) at mesh
vertices. The values at these DoFs can be encapsulated in a PETSc Vec. The new
functionality introduced by this paper provides routines for the construction
(e.g. DMPlexMetricCreate) and modification (e.g. DMPlexMetricEnforceSPD) of
such objects, as well as for combining different instances (e.g.
DMPlexMetricIntersection). We make use of the $L^{p}$ normalization techniques
introduced by [14], which both ensure that resulting adapted meshes have the
desired complexity and enable them to be appropriately multi-scale
(DMPlexMetricNormalize).

Given a metric constructed using these routines, DMAdaptMetric provides the
interface so that ParMmg can be used to drive the mesh adaptation process. The
Pragmatic interface described in [7] is also maintained. When invoked, the Vec
containing the metric data is converted to an array of PetscScalars (i.e.
doubles in real mode) and passed to the appropriate remesher. In addition to
the metric operations above, Mmg and ParMmg have in-built metric gradation
functionality, which ensures that two adjacent edges in the mesh do not differ
by too large a ratio.

#### Hessian Metric

An appropriately scaled Hessian is a common basis for the construction of an
error metric, since it incorporates anisotropic information related to the
curvature of the solution and is related to interpolation error [14]. Its
symmetry implies an orthogonal eigendecomposition at each DoF, which can be
computed straightforwardly. A metric can then be constructed by simply taking
the eigenvalues in modulus. This is justified because error magnitudes are
typically more important than their signs.

#### Integration into an Existing Workflow

PETSc’s nonlinear solver module, SNES, has inbuilt capacity for mesh
adaptation. It allows the user to take a fixed point iteration approach,
whereby the PDE is solved, a Hessian is recovered from its solution and the
mesh is adapted w.r.t. the Hessian at each iteration. For a sufficiently well-
behaved steady-state PDE and appropriate metric parameters, the process
typically converges in a small number of iterations. The timestepping object,
TS, has a similar facility.

The Hessian recovery step is implemented using two applications of the Clément
interpolant [15] – one to recover the gradient and another for the Hessian.
The parallel computation of Clément interpolation requires that the _star_ of
every vertex is available, i.e. the set of all cells containing the vertex.
DMPlex is able to generate an arbitrary mesh overlap in parallel using
DMPlexDistributeOverlap. Then we project the average gradient over the star
into $\mathbb{P}1$. A metric is constructed from the Hessian using
$L^{p}$-normalization, with the default $p=1$ allowing multi-scale meshing.

If the executable for a particular run involving nonlinear solves is called ex
then the application of metric-based mesh adaptation with n fixed point
iterations can be as simple as using the command line options

    ./ex -dm_adaptor mmg -snes_adapt_sequence <n>.

To run in parallel with np MPI processes, simply call

    mpiexec -np <np> ./ex -dm_adaptor parmmg \
                          -snes_adapt_sequence <n>.

##  3 Demonstration

To demonstrate the coupling between PETSc and ParMmg, we solve the Poisson
equation in parallel on adaptive meshes of the unit cube, $\Omega=[0,1]^{d}$,
$d=3$:

$\Delta u=f\quad\text{in}\>\Omega,\qquad u=g\quad\text{on}\>\partial\Omega,$
(1)

with solution $u\in H^{2}(\Omega)$, forcing $f\in L^{2}(\Omega)$ and Dirichlet
condition defined by $g\in L^{2}(\partial\Omega)$. By choosing different $f$
and $g$, we may construct manufactured solutions with different features. We
discretize by seeking weak solutions in order $k\in\mathbb{N}$ Lagrange space,
$\mathbb{P}k\subset H^{1}(\Omega)$. The resulting linear systems are solved
using conjugate gradients, preconditioned with SOR.

###  3.1 Hessian of a Quadratic Function

Suppose we seek the quadratic manufactured solution,

$u_{1}(\mathbf{x}):=\tfrac{2}{3}\mathbf{x}\cdot\mathbf{x},\quad\mathbf{x}:=(x,y,z).$
(2)

This can be achieved with $f\equiv 4$ and appropriate Dirichlet conditions.
Three fixed point iterations are performed, to reduce dependence on the
initial mesh.

The mesh plot in Figure 1 verifies that the uniform Hessian can be recovered
and that the resulting adapted mesh is indeed isotropic. Some elements are
moderately anisotropic (with aspect ratios greater than two), although this is
only true for 0.6% of the elements. The presence of some slightly anisotropic
elements is unavoidable in tetrahedral mesh adaptation since unit tetrahedra
cannot tile space in general.

Figure 1:  Left: mesh generated by adapting w.r.t. the Hessian of the weak
solution of (1) with RHS $f\equiv 4$. Right: $L^{2}$ error convergence
analysis. Mesh statistics: 88,362 elements, 18,534 vertices, max. aspect ratio
5.2981, mean aspect ratio 1.4051, std. dev. 0.1898.

By increasing target metric complexity, we may produce increasingly refined
meshes and thereby perform convergence analysis. The right-hand plot in Figure
1 demonstrates the $(k+1)/d$ convergence rate in the $L^{2}$ norm for degree
$k=1$. This corresponds to the expected $\mathcal{O}(h^{k+1})$ convergence on
isotropic meshes, which is matched by standard uniform refinement.

###  3.2 Interface Capturing

Next, we demonstrate the capability to capture interfaces. We seek the
manufactured solution,

$u_{2}(\mathbf{x}):=\tanh(\alpha(r^{2}-(\mathbf{x}-\mathbf{x}_{0})\cdot(\mathbf{x}-\mathbf{x}_{0}))),$
(3)

which is a smooth approximation of an indicator function for a sphere of
radius $r>0$, centered on $\mathbf{x}_{0}:=(0.5,0.5,0.5)$. The value of
$\alpha>0$ determines the sharpness of the interface is between the exterior
and the interior of the sphere. Take $\alpha=500$ and $r=0.15$.
Differentiation implies the forcing

$f(\mathbf{x}):=2\alpha(u_{2}^{2}(\mathbf{x})-1)(4\alpha(\mathbf{x}-\mathbf{x}_{0})\cdot(\mathbf{x}-\mathbf{x}_{0})\,u_{2}(\mathbf{x})+d).$
(4)

The left-hand plot in Figure 2 shows the adapted mesh for target complexity
32,000. High resolution is used to capture the interface, with some fairly
anisotropic elements. Notably coarse resolution is used near the domain
boundaries. The caption indicates that we have a multi-scale mesh, whose
element volumes span eight orders of magnitude. Elements are generally more
anisotropic for this test case, with 14.8% having aspect ratio greater than
two.

Figure 2:  Left: mesh generated by adapting w.r.t. the Hessian of the weak
solution of (1) with RHS (4). Right: $L^{2}$ error convergence analysis. Mesh
statistics: 595,594 elements, 105,960 vertices, max. aspect ratio 21.0797,
mean aspect ratio 1.6564, std. dev. 0.3907, min./max. volumes: $3.7504\times
10^{-10}$/$1.3294\times 10^{-3}$.

The Hessian-based curves in Figure 2 exhibit the same $(k+1)/d$ convergence
rate as in Subsection 3.1. However, uniform refinement does not attain the
optimal rate for this particular problem. As such, we demonstrate that,
together, PETSc and ParMmg are able to achieve optimal $L^{2}$ convergence for
a Poisson problem with a spatially varying RHS, for which standard uniform
refinement is unable to do this.

##  4 Outlook

The developments described herein comprise the first phase of a project that
aims to bring metric-based mesh adaptation to the user communities of both
PETSc and the downstream Firedrake finite element library [11]. The next phase
is to facilitate mesh adaptation in Firedrake, by implementing consistent
mesh-to-mesh interpolation schemes in parallel, perhaps based on those
currently in PETSc. The project focuses on _goal-oriented_ adaptation, which
enables the accurate approximation of diagnostic quantities using relatively
few DoFs overall. Such technologies are currently being implemented as part of
this project, based on previous proof-of-concept work including [16]. The
final phase will include performing benchmarking experiments and writing
extensive documentation.

## Appendix

The numerical experiments conducted in this paper use one of PETSc’s standard
test cases, which can be found at $PETSC_DIR/src/snes/tutorials/ex12.c in the
source code [4]. Having compiled this example, adaptation to a uniform metric
(as in Subsection 3.1) can be done by running

    mpiexec -np <np> ./ex12 -run_type full
        -dm_plex_dim 3 -dm_distribute
        -dm_plex_box_faces 10,10,10
        -bc_type dirichlet -petscspace_degree <k>
        -variable_coefficient none
        -ksp_type cg -pc_type sor
        -snes_adapt_sequence <n>
        -adaptor_target_num <target>
        -dm_plex_metric_h_max 0.5
        -dm_adaptor parmmg

Here np and n again encode the number of processors and the number of fixed
point iterations in the adaptation loop, whilst k and target encode the
polynomial degree for the finite element space and the target metric
complexity. The convergence analysis experiments may be performed by varying
these values. For the interface capturing example in Subsection 3.2, simply
replace with

        -variable_coefficient ball

in the appropriate line above.

Acknowledgments

This work was funded under the embedded CSE program of the ARCHER2 UK National
Supercomputing Service (http://www.archer2.ac.uk). Support is also
acknowledged from EPSRC under grant EP/R029423/1.

## References

  * [1] Piggott M., Farrell P., Wilson C., Gorman G., Pain C.  “Anisotropic mesh adaptivity for multi-scale ocean modelling.”  _Philos T R Soc A_ , vol. 367, no. 1907, 4591–4611, 2009
  * [2] Alauzet F., Loseille A.  “A decade of progress on anisotropic mesh adaptation for computational fluid dynamics.”  _Comput Aided Des_ , vol. 72, 13–39, 2016
  * [3] Balay S., Abhyankar S., Adams M.F., Benson S., Brown J., Brune P., Buschelman K., Constantinescu E., Dalcin L., Dener A., Eijkhout V., Gropp W.D., Hapla V., Isaac T., Jolivet P., Karpeev D., Kaushik D., Knepley M.G., Kong F., Kruger S., May D.A., McInnes L.C., Mills R.T., Mitchell L., Munson T., Roman J.E., Rupp K., Sanan P., Sarich J., Smith B.F., Zampini S., Zhang H., Zhang H., Zhang J.  “PETSc/TAO Users Manual.”  Tech. Rep. ANL-21/39 - Revision 3.16, Argonne National Laboratory, 2021
  * [4] Balay S., Abhyankar S., Adams M.F., Brown J., Brune P., Buschelman K., Dalcin L., Dener A., Eijkhout V., Gropp W.D., Karpeyev D., Kaushik D., Knepley M.G., May D.A., McInnes L.C., Mills R.T., Munson T., Rupp K., Sanan P., Smith B.F., Zampini S., Zhang H., Zhang H.  “PETSc Web page.”, 2021.  URL www.mcs.anl.gov/petsc
  * [5] Knepley M.G., Karpeev D.A.  “Mesh algorithms for PDE with Sieve I: Mesh distribution.”  _Sci Program_ , vol. 17, no. 3, 215–230, 2009
  * [6] Lange M., Mitchell L., Knepley M.G., Gorman G.J.  “Efficient mesh management in Firedrake using PETSc-DMPlex.”  _SIAM J Sci Comp_ , vol. 38, no. 5, S143–S155, 2016
  * [7] Barral N., Knepley M.G., Lange M., Piggott M.D., Gorman G.J.  “Anisotropic mesh adaptation in Firedrake with PETSc DMPlex.”  _25th Intl Meshing Roundtable_. 2016.  URL https://arxiv.org/pdf/1610.09874.pdf
  * [8] Dobrzynski C., Frey P.  “Anisotropic Delaunay mesh adaptation for unsteady simulations.”  _17th Intl Meshing Roundtable_ , pp. 177–194. Springer, 2008
  * [9] Cirrottola L., Froehly A.  “Parallel unstructured mesh adaptation using iterative remeshing and repartitioning.”  Research Report RR-9307, INRIA Bordeaux, équipe CARDAMOM, Nov. 2019\.  URL https://hal.inria.fr/hal-02386837
  * [10] Dalcin L.D., Paz R.R., Kler P.A., Cosimo A.  “Parallel distributed computing using Python.”  _Adv Water Resour_ , vol. 34, no. 9, 1124–1139, 2011
  * [11] Rathgeber F., Ham D.A., Mitchell L., Lange M., Luporini F., McRae A.T., Bercea G.T., Markall G.R., Kelly P.H.  “Firedrake: automating the finite element method by composing abstractions.”  _ACM Trans Math Softw_ , vol. 43, no. 3, 24, 2016
  * [12] Knepley M.G., Lange M., Gorman G.J.  “Unstructured Overlapping Mesh Distribution in Parallel.”, 2015\.  URL https://arxiv.org/abs/1506.06194
  * [13] Hapla V., Knepley M.G., Afanasiev M., Boehm C., van Driel M., Krischer L., Fichtner A.  “Fully Parallel Mesh I/O using PETSc DMPlex with an Application to Waveform Modeling.”  _SIAM J Sci Comp_ , vol. 43, no. 2, C127–C153, 2021
  * [14] Loseille A., Alauzet F.  “Continuous mesh framework part II: validations and applications.”  _SIAM J Numer Anal_ , vol. 49, no. 1, 61–86, 2011
  * [15] Clément P.  “Approximation by finite element functions using local regularization.”  _Rev fr autom inform rech opér, Anal numér_ , vol. 9, no. R2, 77–84, 1975
  * [16] Wallwork J.G., Barral N., Ham D.A., Piggott M.D.  “Anisotropic Goal-Oriented Mesh Adaptation in Firedrake.”  _28th Intl Meshing Roundtable_ , pp. 83–100. Zenodo, 2020.  URL https://doi.org/10.5281/zenodo.3653101