# AdjointNet: Constraining machine learning models with physics-based codes

Satish Karra1,∗, Bulbul Ahmmed1, and Maruti K. Mudunuru2

1Computational Earth Science Group (EES-16), Earth and Environmental Sciences
Division, Los Alamos National Laboratory, Los Alamos, NM 87545
2Watershed & Ecosystem Science, Pacific Northwest National Laboratory,
Richland, WA 99352

∗Corresponding author<EMAIL_ADDRESS>Last updated: ; LA-UR-21-28763;
PNNL-SA-158825

## Abstract

Physics-informed Machine Learning has recently become attractive for learning
physical parameters and features from simulation and observation data.
However, most existing methods do not ensure that the physics, such as balance
laws (e.g., mass, momentum, energy conservation), are constrained. Some recent
works (e.g., physics-informed neural networks) softly enforce physics
constraints by including partial differential equation (PDE)-based loss
functions but need re-discretization of the PDEs using auto-differentiation.
Training these neural nets on observational data showed that one could solve
forward and inverse problems in one shot. They evaluate the state variables
and the parameters in a PDE. This re-discretization of PDEs is not necessarily
an attractive option for domain scientists that work with physics-based codes
that have been developed for decades with sophisticated discretization
techniques to solve complex process models and advanced equations of state.
This paper proposes a physics constrained machine learning framework,
AdjointNet, allowing domain scientists to embed their physics code in neural
network training workflows. This embedding ensures that physics is constrained
everywhere in the domain. Additionally, the mathematical properties such as
consistency, stability, and convergence vital to the numerical solution of a
PDE are still satisfied. We show that the proposed AdjointNet framework can be
used for parameter estimation (and uncertainty quantification by extension)
and experimental design using active learning. The applicability of our
framework is demonstrated for four cases – (1) flow in a homogeneous porous
medium, (2) data assimilation for homogeneous porous media flow, (3) flow in a
heterogeneous porous medium, and (4) cavity flow using the Navier-Stokes
equation. Results show that AdjointNet-based inversion can estimate process
model parameters with reasonable accuracy. These examples demonstrate the
applicability of using existing software with no changes in source code to
perform accurate and reliable inversion of model parameters.

Keywords: machine learning, physics constraints, partial differential
equations, adjoint, neural networks, inversion.

##  1\. Introduction

State-of-the-art physics-informed machine learning (PIML) [1, 2, 3] or
knowledge-guided machine learning (KGML) [4, 5, 6] approaches have been
primarily data-driven. PIML model development involves data generation from
physics-based codes. Then, unsupervised (e.g., NMFk, NTFk, PCA, autoencoders)
[7] or supervised (e.g., deep neural networks) [1] learning methods are used
to train on this data to obtain state variables or quantities of interest.
These approaches decompose complex coupled processes into various
features/signals that help domain scientists identify the unique physical
mechanisms [8, 7]. Furthermore, popular PIML methods such as physics-informed
neural networks (PINNs) [9, 10, 11] softly (e.g., in $L^{2}$-norm) ensure
physics constraints by training the PDE variables as neural networks. Auto-
differentiation (e.g., in TensorFlow [12], PyTorch [13, 14]) is used to
discretize the PDE. One then uses a weighted loss function (e.g., based on
mean squared error) that linearly combines the loss due to PDE residual and
discrepancies from data. This loss function ensures that the governing laws
are ‘reasonably satisfied’ although not entirely accurate locally at
mesh/collocation points in the domain [15].

Existing PIML and KGML approaches also require large volumes of data to make
reasonably accurate predictions. As a result, ML model development becomes
computationally intensive and time-consuming as large data needs to be
generated. For example, it takes several hours to a day to run a single at-
scale subsurface reservoir simulation [16] with degrees-of-freedom in the
$\mathcal{O}(10^{7})$ on state-of-the-art high-performance computing (HPC)
machines [17, 18, 19]. For this reason, existing PIML/KGML methods are not
ideal for approximating complex and natural systems, which require at least
hundreds of realizations as training data. Furthermore, for applications where
heterogeneity is present in the systems, for instance, different materials in
different parts of the domain, optimizing the number of parameters can be
significant. For instance, in subsurface flow applications, properties such as
permeability, a measure of how much a porous material can allow flow, can be
different in each grid cell of the model. In such scenarios, traditional ML or
PIML methods require Monte Carlo approaches (e.g., MCMC) to sample from the
parameter spaces and generate data [20]. These requirements make data
generation even more challenging.

Most PIML/KGML methods do not include simulation data and experimental
observations in the loop during the ML training process. Training is usually a
two-step process, where an ML model is first built from simulations. Then
real-time data streams are integrated later to infer system decision
parameters (as shown in Fig. (1)(a)). If new data is assimilated beyond the
training range, one needs to generate more training data from simulations and
then re-train the ML model. Our proposed method overcomes this serious
technical gap by directly incorporating physics-based codes into the ML
training procedure. We can use our approach to perform ML while ensuring that
a given physics constraint (e.g., the balance of mass, the balance of
momentum) is satisfied everywhere, as modeled by a physics-based code. Our
innovation is the use of adjoint sensitivities, which can be obtained directly
from physics calculations, to calculate the gradients of the loss function in
the training process of a neural network. Since one needs to solve the forward
problem using the physics code to calculate these adjoint sensitivities, the
underlying physics will be ‘automatically’ satisfied and constrained. Due to
the ‘correct’ constraining of physics in the entire domain, one can still get
highly predictive ML models with limited, sparsely-sampled, and noisy data.

Our computational approach (see Fig. (1)(b)), AdjointNet, allows domain
scientists to incorporate simulators into scalable ML workflows without any
loss of physics (up to the accuracy of a simulator) and simultaneously
integrate observational data. Moreover, our method can directly assimilate
real-time data streams in the training step. The trained ML models can then be
used to forecast system behavior or inform system optimization. Furthermore,
this approach has the advantage that one can invert the unknown parameters in
training the ML model to measured data. We can use these parameters to
forecast the system behavior towards further decision-making. In this paper,
we demonstrate the utility of this AdjointNet framework by applying it to
inversion and learn model parameters in porous media flow and Navier-Stokes
flow problems, where the balance laws – the balance of mass, the balance of
momentum – are solved. The following specific examples are solved: (1) flow in
a homogeneous porous medium, (2) data assimilation for homogeneous porous
media flow, (3) flow in a heterogeneous porous medium, and (4) cavity flow
using the Navier-Stokes equation.

Below, we summarize the novelty of the proposed AdjointNet framework compared
to current state-of-the-art:

  1. (1)

Incorporating existing codes and extendability: Our framework can incorporate
any existing physics-based codes such as PFLOTRAN[21, 22], E3SM [23], OpenFOAM
[24], SWAT [25], PRMS[26], HOSS [27], and WRF-Hydro [28]. For instance, these
codes can simulate complex coupled processes such as flow, chemical and
thermal transport, and mechanics for various applications such as climate,
subsurface energy, and carbon and nuclear waste storage, etc. Most of the
above codes have undergone rigorous verification and validation over the
years. As a result, we can simulate, estimate, and invert for material or
conceptual properties with no changes in existing codes reliably. This
provides a significant advantage over existing methods such as PINNs [9],
where process models and equations of state need to be re-written to perform
the inversion.

  2. (2)

Balance laws: Even though the trained PIML models seemingly give reasonable
predictions on quantities of interest such as pressures, displacements, and
temperatures, there are no guarantees that the governing laws of physics
(e.g., the balance of mass, energy, or momentum) will truly be satisfied by
the machine learning (ML) model. That is, the model can violate these laws
locally or globally, or both, as reported by the developers of PINNs [15]. For
example, a significant limitation of PINNs [29, 30, 15] is the accuracy of the
solution and associated local balance. In PINNs, the absolute error does not
go below $\mathcal{O}(10^{-5})$ [15]. This is because of the reduced accuracy
in solving the high-dimensional non-convex optimization problem, resulting in
local minima. Since the AdjointNet framework directly calls the existing
codes, which ensure that the balance laws are satisfied up to machine
precision, the final trained model will always be locally conservative [31,
32].

  3. (3)

Scalability: Over the decades, some of these codes mentioned above have been
parallelized to scale on leadership-class HPC machines [21, 22, 23]. Our
AdjointNet framework can be scalable as it uses the existing HPC codes under
the hood. Moreover, the proposed algorithm uses discrete adjoint-based
methods, which are also scalable [33, 34]. Hence, the overall approach can be
scaled on the next generation HPC machines geared towards exascale computing
(e.g., Perlmutter, Aurora) [35, 36, 37]. This provides a significant advantage
over other PIML methods (e.g., PINNs).

  4. (4)

Agnostic to underlying numerical discretization: AdjointNet works with any
numerical method implemented in a physics code. Here, we demonstrate the
utility with two physics codes implemented with different numerical
methods–finite volume and finite difference–but can be extended to any
numerical method of choice.

  5. (5)

Differentiable programming: Recently, programming languages like Julia have
made differentiable programming more accessible [38]. However, the governing
equations in physics codes have to be re-written in Julia to take advantage of
this capability. As mentioned previously, this re-writing of codes may not be
attractive for domain scientists. AdjointNet leverages the infrastructure in
Tensorflow or PyTorch to perform auto-differentiation on the neural networks
but does not require differentiable programming of the underlying physics
codes.

  6. (6)

Computational cost: Another major limitation of PINNs in addition to
inaccuracy is the significant training cost associated with deep neural
networks and long-time integration of the PDEs [15]. This is because the
addition of a PDE-based loss function results in a stiff and non-convex
optimization problem [39, 40]. Our AdjointNet framework does not do this but
poses ML training as a PDE-constrained optimization problem. Furthermore,
since we can directly integrate existing HPC physics-based codes, we can
obtain the final trained model faster, thus allowing us to work with more
realistic and large-scale problems.

  7. (7)

Desirable mathematical properties: Most of the existing physics codes satisfy
the desirable properties such as structure-preserving, coercivity, error
estimates, consistency, convergence, accuracy, and stability [41, 31].
Whichever of these properties are satisfied by the underlying PDE solver used
in the AdjointNet framework, the final trained model will also satisfy that.
On the other hand, existing PIML methods do not ensure this.

  8. (8)

eXplainable Artificial Intelligence (XAI): Typically, numerical methods like
finite element, finite volume, etc., [41, 31] are employed to obtain the
solution of a PDE. These numerical methods are explainable, which results in
the trustworthiness of the obtained solution. On the other hand, existing PIML
methods such as PINNs use neural networks to compute the solution of the
variables in a PDE. These neural networks are difficult to interpret, and XAI
methods [42, 43] (e.g., SHAPley values [44], LRP [45], deep Taylor
decomposition [46], LIME [47]) are needed to solve the solution. When neural
networks are employed to obtain the solution, this interpretation process
needs to be performed for each scenario of interest (e.g., changes in boundary
conditions, initial conditions, or domain). However, note that most of the
PIML related works have rarely used XAI methods to explain why such a method
works.

  9. (9)

Reproducability: Typically, in any PIML method, the solution to the PDE
depends on the trainable weights that need to be initialized randomly [48, 49,
50]. As a result, the procedure of weights initialization plays a predominant
role in the obtained solution, leading to inconsistency and hence may not be
reproducible. However, our AdjointNet is not affected by such initialization
as we rely on the existing physics codes to always provide reproducible
solutions of the PDE.

The outline of our paper is as follows: Sec. (2) presents the mathematical
formulation and the linkage between physics-based codes, which is the
AdjointNet framework. Section (3) presents examples of our proposed approach.
Conclusions are drawn in Sec. (4).

##  2\. AdjointNet Framework

Machine learning methods (e.g., neural networks) typically involve minimizing
a loss function $\mathcal{L}$. We add a constraint that the PDE of interest
needs to be satisfied in our underlying formulation. This turns into a PDE-
constrained optimization problem, that is,

$\displaystyle\textrm{min}~{}\mathcal{L}$ (2.1)
$\displaystyle\textrm{s.t.}~{}\mathcal{F}\left(\boldsymbol{u}(\boldsymbol{x},t);\boldsymbol{p}\right)=\boldsymbol{0}$
(2.2)

where $\boldsymbol{x}$ is the position, $\boldsymbol{u}$ is the variable of
interest (e.g., pressure, temperature, concentration), $\boldsymbol{p}$ is the
vector of parameters (e.g., permeability, thermal conductivity, diffusivity),
and $t$ is time. Here, Eq. (2.2) is the discrete form of the PDE that is
solved by a physics-based code. For instance, in the case of porous media, the
parallel code PFLOTRAN solves the discrete form of the non-linear diffusion
PDE [22]:

$\displaystyle\dfrac{\partial{\left[\phi\rho(u)u\right]}}{\partial{t}}-\mathrm{div}\left[\frac{\rho(u)k}{\mu}\nabla{u}\right]=0$
(2.3)

where $u$ is the fluid pressure, $k$ is the permeability parameter, $\phi$ is
porosity, $\mu$ is fluid viscosity, and $\rho$ is the fluid density. PDE
discretization in PFLOTRAN is performed using a two-point flux finite volume
method with backward Euler for time-stepping. The resulting nonlinear
algebraic equations are solved using a Newton-Krylov solver.

If we assume that the parameters $\boldsymbol{p}$ are to be learnt from a set
of observational data $\boldsymbol{u}_{\mathrm{obs}}$, then one can use a loss
function of $\mathcal{L}=\|u-u_{obs}\|$. Also, suppose that $\boldsymbol{p}$
is a ML model such as a neural network, i.e.,
$\boldsymbol{p}=NN(\boldsymbol{x})$. If $W^{(k)}$ and $b^{(k)}$ are the
weights and biases for this neural network $NN(\boldsymbol{x})$, then the
gradients $\dfrac{\partial\mathcal{L}}{\partial W^{(k)}}$,
$\dfrac{\partial\mathcal{L}}{\partial b^{(k)}}$ are needed. These gradients
must be supplied to the optimization algorithm that minimizes the loss
function; for instance, stochastic gradient descent is a commonly used
algorithm. Now, using the chain rule, one can write these gradients as

$\displaystyle\dfrac{\partial\mathcal{L}}{\partial
W^{(k)}}=\dfrac{\partial\mathcal{L}}{\partial\boldsymbol{u}}\dfrac{\partial\boldsymbol{u}}{\partial\boldsymbol{p}}\dfrac{\partial\boldsymbol{p}}{\partial
W^{(k)}},~{}\dfrac{\partial\mathcal{L}}{\partial
b^{(k)}}=\dfrac{\partial\mathcal{L}}{\partial\boldsymbol{u}}\dfrac{\partial\boldsymbol{u}}{\partial\boldsymbol{p}}\dfrac{\partial\boldsymbol{p}}{\partial
b^{(k)}}$ (2.4)

Notice that in both sets of gradients, we need the sensitivity
$\dfrac{\partial\boldsymbol{u}}{\partial\boldsymbol{p}}$, calculated through
the adjoint form of Eq. (2.2), which is obtained by taking the derivative of
Eq. (2.2) with respect to $\boldsymbol{p}$:

$\displaystyle\dfrac{\partial\mathcal{F}}{\partial\boldsymbol{p}}+\dfrac{\partial\mathcal{F}}{\partial\boldsymbol{u}}\dfrac{\partial\boldsymbol{u}}{\partial\boldsymbol{p}}=\boldsymbol{0}$
(2.5)

We have dropped time $t$ in Eq. (2.5) assuming steady-state, for the sake of
illustration. However, in the case of transient problems, if one discretizes
the PDE with time first, we get a similar form for the adjoint equation Eq.
(2.5). In each epoch of the optimization loop, one needs to solve Eq. (2.5) to
get the sensitivity. The other terms in the chain rule, the derivatives
$\dfrac{\partial\mathcal{L}}{\partial\boldsymbol{u}}$,
$\dfrac{\partial\boldsymbol{p}}{\partial b^{(k)}}$,
$\dfrac{\partial\boldsymbol{p}}{\partial W^{(k)}}$ can be obtained directly
through auto-differentiation; ML packages like TensorFlow or PyTorch provide
them. One can use the discrete adjoint method [51] to solve Eq. (2.5).

Discrete adjoint methods involve solving the forward problem and checkpointing
the solution; then the adjoint (backward) solve is performed by integrating
Eq. (2.5) using the saved forward solution at checkpoints. Hence, the cost of
an adjoint solve is approximately the cost of a forward solve, which is
independent of the number of parameters. For this reason, discrete adjoint
methods are popular in the area of design, for example, in shape optimization
of airfoils in the aerospace industry [52, 53]. Since in the discrete adjoint
method to solve Eq. (2.5), one must solve the forward problem Eq. (2.2), we
ensure that the physics is constrained within each epoch. Alternatively, one
can obtain the sensitivities
$\dfrac{\partial\boldsymbol{u}}{\partial\boldsymbol{p}}$, by perturbing each
parameter in $\boldsymbol{p}$, re-solving the PDE Eq. (2.2), and then
calculating the gradients numerically. Even for this approach, since one needs
to solve the PDE, the physics is thus constrained in each epoch. However, the
cost of evaluating gradients with the proposed approach is $N_{p}+1$ forward
solves, where $N_{p}$ is the number of parameters. Thus as the unknown
parameters increase, this approach becomes relatively expensive. An approach
to overcome this problem is to integrate prior information (e.g., domain
expertise, local and global sensitivity analysis, mutual information theory)
with unsupervised learning (e.g., $k$-mean clustering) to reduce the
dimensions of unknown parameters before applying our AdjointNet framework.
This provides an efficient way to parameterize the high-dimensional space of
unknown parameters before performing inversion efficiently.

(a) State-of-the-art ML workflow (b) Proposed workflow

Figure 1. AdjointNet framework: Comparison between the state-of-the-art ML
workflow with the proposed workflow. A major advantage is that the proposed
workflow requires minimal simulations, as it calls the physics-based code on
the fly, to perform data assimilation and ML training.

##  3\. Numerical Examples

In this section, we present representative numerical examples to demonstrate
the applicability of our AdjointNet framework for parameter inversion and data
assimilation.

###  3.1. Porous media flow: Inversion of homogeneous permeability under
sparse data

Through a synthetic case, we perform porous media simulations using PFLOTRAN
to invert for homogeneous permeability. PFLOTRAN is an open-source, state-of-
the-art massively parallel subsurface flow and reactive transport code [21,
22]. We use the pressure data sampled at different locations in the domain as
ground truth to perform the inversion. The model domain is $100\times 1\times
1$ ${\mathrm{m}}^{3}$ in dimension and consists of one layer. A uniform mesh
of 100 grid cells is used to generate data. The permeability value $10^{-14}$
m2 is assigned for all grids, which is assumed to be the ground truth. The
initial condition for the PFLOTRAN simulation includes the pressure of
1.5$\times 10^{5}$ Pa throughout the model domain. Flow is driven by pressure
boundary conditions, 1.0$\times 10^{6}$ Pa and 1.5$\times 10^{5}$ Pa at the
left and right faces, respectively.

(a) Permeability estimation vs. epochs (b) Pressure profile prediction and
comparison

Figure 2. AdjointNet for inversion of homogeneous permeability: Training of
permeability using AdjointNet. As the model trains to the pressure data (with
noise added) as shown in (b), the model learns the permeability as shown in
(a). (b) also compares the prediction using the trained permeability against
the original pressure data.

To check the validity of the proposed framework, first, we generate ground
truth pressure data using the homogeneous numerical simulation. Next, we
sparsely sampled ten out of 100 data points for training the permeability
neural network. Random noise (with a magnitude of $\approx 0.1$ MPa) is added
to the training data to perform inversion under noisy data. A permeability
neural network is created that has two layers containing 256 neurons in each
layer. The loss function is formulated based on the mean-squared error between
the true and predicted pressure values at these ten points. The learning rate
is set to $10^{-3}$, and the neural network weights and biases are initialized
based on Glorot uniform initialization. The homogeneous permeability neural
network is trained for 100 epochs, and convergence is reached within 30 epochs
as shown in Fig. (2)(a). During each epoch, the neural network is updated as
the flow process is simulated by PFLOTRAN. Neural network provides gradient of
permeability with respect to weights $\dfrac{\partial\boldsymbol{p}}{\partial
W}$ while the numerical model provides
$\dfrac{\partial\mathcal{L}}{\partial\boldsymbol{u}}$ and
$\dfrac{\partial\boldsymbol{u}}{\partial\boldsymbol{p}}$. The loss value
plateaued to 0.0015. Figure (2)(a) shows that the AdjointNet estimated
permeability is very close to the ground truth. Figure (2)(b) compares the
predicted pressures to ground truth and shows that pressure data predictions
are robust to noise compared to actual values in the entire domain. This case
study instills confidence that AdjointNet methodology can provide reasonably
accurate permeability estimates even under sparsely sampled data and its model
predictions are robust to observational noise.

(a) Phase 1 (b) Phase 2 (c) Phase 3 (d) Permeability training

Figure 3. AdjointNet for data assimilation with homogeneous permeability:
Training of permeability using a neural network as data is assimilated in
three sequential phases shown in (a), (b), and (c). The trained permeability
is shown in (d).

###  3.2. Porous media flow: Data assimilation

This case study shows a proof-of-concept towards data assimilation
applications. The simulation setup is the same as the previous use case.
However, the homogeneous permeability neural network is trained in three
sequential phases (Phase 1, Phase 2, and Phase 3). The neural network is first
trained against two randomly selected pressure data points with added random
noise (with a magnitude of $0.1$ MPa) in Phase 1 (Fig. (3)(a)). This is
followed by further training by assimilating an additional data point in each
of Phase 2 (Fig. (3)(b)) and Phase 3 (Fig. (3)(c)). Training in each phase is
run for 200 epochs, which totals 600 epochs. Figure (3)(d) shows the
permeability predictions for three phases and compares them with the ground
truth. It is evident that as we progressively add more pressure data, we can
better constrain permeability. As a result, Phase 3 permeability prediction is
closer to the ground truth than Phase 2 and Phase 1 (see Fig. (3)(d)).

###  3.3. Porous media flow: Heterogeneous permeabilities

For this case study, the model domain consists of two zones with different
permeabilities. The domain is meshed in the same manner as in Sec. (3.1). The
total number of grid cells is equal to 100. The initial condition for the
PFLOTRAN simulation includes the pressure of 1.5$\times 10^{5}$ Pa throughout
the model domain. Pressure boundary conditions are prescribed on the left and
right sides of the model domain, which are 1.0 and 0.5 MPa. For ground truth,
we assume the first and last $50$ grid cells have permeability values of
$10^{-13}$ and $10^{-15}$ m2, respectively. Based on these permeability
values, we generate pressure data for AdjointNet-enabled inversion. The solid
red line in Fig. (4) shows the ground truth pressure distribution. Two neural
networks with two layers and 256 neurons in each layer are trained to estimate
the permeabilities until a loss value of $10^{-8}$ is reached. The training
process converged within 100 epochs, as shown in Fig. (5). The AdjointNet
learned permeabilities were able to predict the pressure distribution
accurately, as seen in Fig. (4).

Figure 4. AdjointNet for inversion of heterogeneous permeability: This figure
shows the actual and AdjointNet predicted pressure distribution over the model
domain. AdjointNet can capture the pressure profiles with reasonably good
accuracy even under discontinuities.

(a) $k_{1}$ training (b) $k_{2}$ training

Figure 5. AdjointNet for inversion of heterogeneous permeability: This figure
shows the training progress of the two permeability neural networks. We can
see that the permeability estimates plateau within 100 epochs.

###  3.4. Lid-driven cavity flow using the Navier-Stokes equation

In this section, we show the applicability of AdjointNet methodology to
simulate 2D lid-driven cavity flow. The domain is $4\times 4$
${\mathrm{m}}^{2}$ discretized by $41\times 41$ uniform grid points. We used a
pressure Poisson finite difference solver written in Python [54]. The initial
condition is $u,v,p=0$ everywhere. The boundary conditions are: $u=1$ at
$y=4$; $u,v=0$ at $y=0$; $\frac{\partial p}{\partial y}=0$ at $y=0$, $p=0$ at
$y=4$, $\frac{\partial p}{\partial x}=0$ at $x=0,4$. The ground truth is
simulated using $\rho$ and $\nu$ to be equal to 1
$\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{{}^{3}}\mathrm{m}$ and 0.1
$\mathrm{{}^{2}}\mathrm{m}\mathrm{/}\mathrm{s}$, respectively.

The model is run for 100 iterations, and for each iteration total time step is
300. Simulations are performed to generate pressure data for AdjointNet-
enabled inversion. A pressure snapshot (Figure (6)) with 1,681 pressure points
is used to learn the kinematic viscosity. A neural network for the kinematic
viscosity is created that has two layers containing 64 and 32 neurons in each
layer. Convergence is reached in around 200 epochs with the loss value of less
than $10^{-8}$ (Fig. (7)). AdjointNet estimated kinematic viscosity is 0.1021
compared to the ground truth of kinematic viscosity 0.1, which shows the
applicability of AdjointNet for inversion of flow parameters with the Navier-
Stokes equation.

Figure 6. AdjointNet for inversion of kinematic viscosity: Lid-driven cavity
Navier-Stokes flow. Colors represent pressure distribution for the ground
truth $\rho$ and $\nu$.

(a) Loss versus epochs (b) True and estimated kinematic viscosity

Figure 7. AdjointNet for inversion of kinematic viscosity: Kinematic viscosity
prediction by AdjointNet (a) shows the progression of loss over epochs and (b)
compares the true and kinematic viscosity as training progresses.

##  4\. Conclusions

Currently, neural network models are black-box models, making it challenging
to interpret and trust the results. Recently, XAI techniques (e.g., deep
Taylor decomposition) were developed to make neural networks interpretable,
but they are not universally applicable (i.e., trustworthiness is always
problem-specific). This paper proposed a new ML methodology called AdjointNet
and demonstrated its utility to solve inverse problems. The proposed
AdjointNet framework overcomes these pitfalls, as the existing codes are
interpretable and physics is constrained everywhere. Through representative
numerical examples, we showed that one could invert material parameters by
embedding a physics-code (e.g., PFLOTRAN, Navier-Stokes solver) in AdjointNet
without modifying the underlying physics-code, or without the need to re-write
a code for a PDE from scratch, which is a significant bottleneck for existing
PIML methods such as PINNs. We demonstrated AdjointNet’s capability on four
examples: (1) homogeneous permeability inversion; (2) data assimilation; (3)
heterogeneous permeability inversion; and (4) estimating material properties
in a Navier-Stokes lid-driven cavity flow.

In addition to satisfying the physics constraint, the results also showed that
we could successfully interpret and trust them as the existing codes are
verified and validated by domain experts over decades. Since AdjointNet
requires sensitivities of the underlying variable (e.g., pressure,
temperature) with respect to the input parameters, this can be a limitation in
cases where there are many parameters, such as spatially distributed
permeability, in the case of porous flow. In such cases, numerically
evaluating the sensitivities by perturbing the input parameters may be
computationally intractable. One will have to resort to the discrete adjoint
method, which is agnostic to the number of model parameters but may need some
intervention to the underlying physics code. Alternatively, one can resort to
a priori dimensionality reduction, where one can integrate local/global
sensitivity analysis with unsupervised learning (e.g., $k$-mean clustering) to
identify the most critical parameters and estimate them using AdjointNet.
Finally, since AdjointNet can learn the underlying parameters while
constraining the physics, it can perform well beyond the training range. As a
result, this framework is attractive for domain scientists who use complex
codes to simulate physical processes.

## ABBREVIATIONS

  * •

HPC: High Performance Computing

  * •

LIME: Local Interpretable Model-Agnostic Explanations

  * •

LRP: Layer-Wise Relevance Propagation

  * •

KGML: Knowledge-Guided Machine Learning

  * •

MCMC: Markov Chain Monte Carlo

  * •

ML: Machine Learning

  * •

NMFk: Non-negative Matrix Factorization with custom $k$-means Clustering

  * •

NTFk: Non-negative Tensor Factorization with custom $k$-means Clustering

  * •

PDE: Partial Differential Equation

  * •

PFLOTRAN: An open source, state-of-the-art massively parallel subsurface flow
and reactive transport code

  * •

PIML: Physics-informed Machine Learning

  * •

PINN: Physics-informed Neural Network

  * •

XAI: eXplainable Artificial Intelligence

## Acknowledgments

The authors thank the U.S. Department of Energy’s Biological and Environmental
Research Program for support through the SciDAC4 program. SK and BA also thank
the Center for Space and Earth Science and the Information Science &
Technology Institute at Los Alamos National Laboratory. MKM thanks the support
from the U.S. Department of Energy Office of Science’s River Corridor Science
Focus Area at PNNL. The views and opinions of authors expressed herein do not
necessarily state or reflect those of the United States Government or any
agency thereof.

## References

  * [1] D. R. Harp, D. O’Malley, B. Yan, and R. Pawar, “On the feasibility of using physics-informed machine learning for underground reservoir pressure management,” _Expert Systems with Applications_ , vol. 178, p. 115006, 2021\.
  * [2] K. Kashinath, M. Mustafa, A. Albert, J. L. Wu, C. Jiang, S. Esmaeilzadeh, K. Azizzadenesheli, R. Wang, A. Chattopadhyay, and A. Singh, “Physics-informed machine learning: Case studies for weather and climate modelling,” _Philosophical Transactions of the Royal Society A_ , vol. 379, p. 20200093, 2021.
  * [3] B. Yan, D. R. Harp, B. Chen, and R. Pawar, “A physics-constrained deep learning model for simulating multiphase flow in 3D heterogeneous porous media,” _arXiv preprint arXiv:2105.09467_ , 2021.
  * [4] A. Karpatne, G. Atluri, J. H. Faghmous, M. Steinbach, A. Banerjee, A. Ganguly, S. Shekhar, N. Samatova, and V. Kumar, “Theory-guided data science: A new paradigm for scientific discovery from data,” _IEEE Transactions on knowledge and data engineering_ , vol. 29, pp. 2318–2331, 2017.
  * [5] P. C. Hanson, A. B. Stillman, X. Jia, A. Karpatne, H. A. Dugan, C. C. Carey, J. Stachelek, N. K. Ward, Y. Zhang, and J. S. Read, “Predicting lake surface water phosphorus dynamics using process-guided machine learning,” _Ecological Modelling_ , vol. 430, p. 109136, 2020.
  * [6] A. Khandelwal, S. Xu, X. Li, X. Jia, M. Stienbach, C. Duffy, J. Nieber, and V. Kumar, “Physics guided machine learning methods for hydrology,” _arXiv preprint arXiv:2012.02854_ , 2020.
  * [7] V. Vesselinov, M. Mudunuru, S. Karra, D. O’Malley, and B. Alexandrov, “Unsupervised machine learning based on non-negative tensor factorization for analyzing reactive-mixing,” _Journal of Computational Physics_ , vol. 395, pp. 85–104, 2019.
  * [8] A. Cichocki, R. Zdunek, A. H. Phan, and S.-I. Amari, _Nonnegative matrix and tensor factorizations: Applications to exploratory multi-way data analysis and blind source separation_. John Wiley & Sons, 2009.
  * [9] M. Raissi, P. Perdikaris, and G. 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_ , vol. 378, pp. 686–707, 2019.
  * [10] A. Tartakovsky, C. Marrero, P. Perdikaris, G. Tartakovsky, and D. Barajas-Solano, “Physics-informed deep neural networks for learning parameters and constitutive relationships in subsurface flow problems,” _Water Resources Research_ , vol. 56, no. 5, p. e2019WR026731, 2020.
  * [11] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis, “Deepxde: A deep learning library for solving differential equations,” _SIAM Review_ , vol. 63, no. 1, pp. 208–228, 2021.
  * [12] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, and M. Isard, “Tensorflow: A system for large-scale machine learning,” in _12th $\\{$USENIX$\\}$ symposium on operating systems design and implementation ($\\{$OSDI$\\}$ 16)_, 2016, pp. 265–283.
  * [13] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer, “Automatic differentiation in pytorch,” 2017.
  * [14] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, and L. Antiga, “Pytorch: An imperative style, high-performance deep learning library,” _Advances in neural information processing systems_ , vol. 32, pp. 8026–8037, 2019.
  * [15] A. D. Jagtap, E. Kharazmi, and G. E. Karniadakis, “Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 365, p. 113028, 2020.
  * [16] M. K. Mudunuru, D. O’Malley, S. Srinivasan, J. D. Hyman, M. R. Sweeney, L. Frash, B. Carey, M. R. Gross, N. J. Welch, S. Karra, V. V. Vesselinov, Q. Kang, H. Xu, R. J. Pawar, T. Carr, L. Li, G. D. Guthrie, and H. S. Viswanathan, “Physics-informed machine learning for real-time reservoir management.” in _AAAI Spring Symposium: MLPS_ , 2020.
  * [17] (2021) NERSC – National Energy Research Scientific Computing Center. Accessed on: 2021-09-03. [Online]. Available: https://www.nersc.gov/
  * [18] (2021) OLCF – Oak Ridge Leadership Computing Facility. Accessed on: 2021-09-03. [Online]. Available: https://www.olcf.ornl.gov/
  * [19] (2021) ALCF – Argonne Leadership Computing Facility. Accessed on: 2021-09-03. [Online]. Available: https://www.alcf.anl.gov/
  * [20] J. E. Santos, M. Mehana, H. Wu, M. Prodanovic, Q. Kang, N. Lubbers, H. Viswanathan, and M. J. Pyrcz, “Modeling nanoconfinement effects using active learning,” _The Journal of Physical Chemistry C_ , vol. 124, pp. 22 200–22 211, 2020.
  * [21] G. Hammond, P. Lichtner, and R. Mills, “Evaluating the performance of parallel subsurface simulators: An illustrative example with PFLOTRAN,” _Water Resources Research_ , vol. 50, no. 1, pp. 208–228, 2014.
  * [22] P. C. Lichtner, G. E. Hammond, C. Lu, S. Karra, G. Bisht, B. Andre, R. Mills, and J. Kumar, “Pflotran user manual: A massively parallel reactive flow and transport model for describing surface and subsurface processes,” Tech. Rep., 2015.
  * [23] J.-C. Golaz, P. M. Caldwell, L. P. Van R., M. R. Petersen, Q. Tang, J. D. Wolfe, G. Abeshu, V. Anantharaj, X. S. Asay-D., and D. C. Bader, “The DOE E3SM coupled model version 1: Overview and evaluation at standard resolution,” _Journal of Advances in Modeling Earth Systems_ , vol. 11, pp. 2089–2129, 2019.
  * [24] H. Jasak, A. Jemcov, and Z. Tukovic, “OpenFOAM: A C++ library for complex physics simulations,” in _International workshop on coupled methods in numerical dynamics_ , vol. 1000, 2007, pp. 1–20.
  * [25] J. G. Arnold, D. N. Moriasi, P. W. Gassman, K. C. Abbaspour, M. J. White, R. Srinivasan, C. Santhi, R. D. Harmel, A. V. Griensven, and M. W. Van Liew, “SWAT: Model use, calibration, and validation,” _Transactions of the ASABE_ , vol. 55, pp. 1491–1508, 2012.
  * [26] G. H. Leavesley and L. G. Stannard, “The precipitation-runoff modeling system-PRMS,” _Computer models of watershed hydrology_ , pp. 281–310, 1995\.
  * [27] E. E. Knight, E. Rougier, Z. Lei, B. Euser, V. Chau, S. H. Boyce, K. Gao, K. Okubo, and M. Froment, “HOSS: An implementation of the combined finite-discrete element method,” _Computational Particle Mechanics_ , vol. 7, pp. 765–787, 2020.
  * [28] K. Sampson and D. Gochis, “RF Hydro GIS Pre-Processing Tools, Version 5.0, Documentation,” _Boulder, CO: National Center for Atmospheric Research, Research Applications Laboratory_ , 2018.
  * [29] O. Fuks and H. A. Tchelepi, “Limitations of physics informed machine learning for nonlinear two-phase transport in porous media,” _Journal of Machine Learning for Modeling and Computing_ , vol. 1, 2020.
  * [30] S. Wang, Y. Teng, and P. Perdikaris, “Understanding and mitigating gradient pathologies in physics-informed neural networks,” _arXiv preprint arXiv:2001.04536_ , 2020.
  * [31] R. Eymard, T. Gallouët, and R. Herbin, “Finite Volume Methods,” _Handbook of Numerical Analysis_ , vol. 7, pp. 713–1018, 2000.
  * [32] T. J. R. Hughes, G. Engel, L. Mazzei, and M. G. Larson, “The continuous Galerkin method is locally conservative,” _Journal of Computational Physics_ , vol. 163, pp. 467–488, 2000.
  * [33] M. Towara, M. Schanen, and U. Naumann, “MPI-parallel discrete adjoint OpenFOAM,” _Procedia Computer Science_ , vol. 51, pp. 19–28, 2015.
  * [34] H. Zhang, E. M. Constantinescu, and B. F. Smith, “PETSc TSAdjoint: A discrete adjoint ODE solver for first-order and second-order sensitivity analysis,” _arXiv preprint arXiv:1912.07696_ , 2019.
  * [35] P. Messina, “The exascale computing project,” _Computing in Science & Engineering_, vol. 19, pp. 63–67, 2017.
  * [36] R. Stevens, J. Ramprakash, P. Messina, M. Papka, and K. Riley, “Aurora: Argonne’s next-generation exascale supercomputer,” ANL (Argonne National Laboratory (ANL), Argonne, IL (United States)), Tech. Rep., 2019.
  * [37] C. Yang and J. Deslippe, “Accelerate Science on Perlmutter with NERSC,” _Bulletin of the American Physical Society_ , vol. 65, 2020.
  * [38] M. Innes, A. Edelman, K. Fischer, C. Rackauckas, E. Saba, V. B. Shah, and W. Tebbutt, “A differentiable programming system to bridge machine learning and scientific computing,” _arXiv preprint arXiv:1907.07587_ , 2019.
  * [39] K. Xu, A. Tartakovsky, J. Burghardt, and E. Darve, “Inverse modeling of viscoelasticity materials using physics constrained learning,” _arXiv preprint arXiv:2005.04384_ , 2020.
  * [40] K. Xu and E. Darve, “Physics constrained learning for data-driven inverse modeling from sparse observations,” _arXiv preprint arXiv:2002.10521_ , 2020\.
  * [41] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu, _The finite element method: Its basis and fundamentals_ , 2005.
  * [42] F. K. Došilović, M. Brčić, and N. Hlupić, “Explainable artificial intelligence: A survey,” in _2018 41st International convention on information and communication technology, electronics and microelectronics (MIPRO)_. IEEE, 2018, pp. 0210–0215.
  * [43] A. Adadi and M. Berrada, “Peeking inside the black-box: A survey on explainable artificial intelligence (XAI),” _IEEE access_ , vol. 6, pp. 52 138–52 160, 2018.
  * [44] M. Sundararajan and A. Najmi, “The many Shapley values for model explanation,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 9269–9278.
  * [45] G. Montavon, A. Binder, S. Lapuschkin, W. Samek, and K.-R. Müller, “Layer-wise relevance propagation:An overview,” _Explainable AI: interpreting, explaining and visualizing deep learning_ , pp. 193–209, 2019.
  * [46] G. Montavon, S. Lapuschkin, A. Binder, W. Samek, and K.-R. Müller, “Explaining nonlinear classification decisions with deep Taylor decomposition,” _Pattern Recognition_ , vol. 65, pp. 211–222, 2017.
  * [47] M. T. Ribeiro, S. Singh, and C. Guestrin, “Model-agnostic interpretability of machine learning,” _arXiv preprint arXiv:1606.05386_ , 2016.
  * [48] S. K. Kumar, “On weight initialization in deep neural networks,” _arXiv preprint arXiv:1704.08863_ , 2017.
  * [49] V. Nagarajan and J. Z. Kolter, “Generalization in deep networks: The role of distance from initialization,” _arXiv preprint arXiv:1901.01672_ , 2019\.
  * [50] Y. Zhang, Z.-Q. J. Xu, T. Luo, and Z. Ma, “A type of generalization error induced by initialization in deep neural networks,” in _Mathematical and Scientific Machine Learning_. PMLR, 2020, pp. 144–164.
  * [51] M. B. Giles, M. C. Duta, J.-D. Muller, and N. A. Pierce, “Algorithm developments for discrete adjoint methods,” _AIAA journal_ , vol. 41, pp. 198–205, 2003.
  * [52] S. K. Nadarajah and A. Jameson, “Optimum shape design for unsteady flows with time-accurate continuous and discrete adjoint method,” _AIAA journal_ , vol. 45, pp. 1478–1491, 2007.
  * [53] M. P. Rumpfkeil and D. W. Zingg, “The optimal control of unsteady flows with a discrete adjoint method,” _Optimization and Engineering_ , vol. 11, pp. 5–22, 2010.
  * [54] (2021) 12 steps to Navier–Stokes. Accessed on: 2021-09-03. [Online]. Available: https://nbviewer.jupyter.org/github/barbagroup/CFDPython/blob/master/lessons/14_Step_11.ipynb


[a]Tomás E. Gonzalo

# BSM Developments and Tools

###### Abstract

In this conference paper I introduce a selection of BSM tools and describe
their most recent developments. I choose to focus on tools for the
reinterpretation of LHC searches, tools that compute dark matter constraints,
beyond-the-Standard-Model (BSM) inference tools, and tools that calculate
amplitudes and cross sections directly from the Lagrangian.

TTK-21-40

##  1 Introduction

We are spoiled for choices when it comes to BSM tools. For many years the
community has developed tools for every aspect of BSM physics, from various
computations of observables, to abstract level tools able to handle the full
Lagrangian of a theory. Due to the large amount of software available, in this
conference article I will only describe a selection of tools and their recent
developments. I will focus on tools for the recasting of LHC searches, dark
matter (DM) tools, statistical inference tools and Lagrangian-level tools.

##  2 Reinterpretation of LHC searches

The reinterpretation of BSM searches has been a profilic field of study in
particle physics. Significant advances have been made in the last few years in
the interaction between the experimental and theoretical communities, and the
exchange of information between them is at its highest, especially after the
release of simplified and full likelihoods by CMS and ATLAS, respectively [1].
This has had profound positive effects on the whole community as well as the
physics outcomes, as it has helped with the long-term preservation of the
results and provides motivation for the design of new analyses. Consequently,
many recasting tools have been developed, following two different approaches:
the Simplified Models approach, or the Fast Simulation approach.

The Simplified Models approach tackles the issue of recasting by decomposing a
given model into simplified model topologies and comparing the predicted cross
section times branching ratio with that reported by the experiments [2]. As it
does not require any simulation of collisions, this approach is extremely fast
and perfectly suitable for large parameter fits. The prominent tool that
implements this approach is SModelS, which is in fact responsible for
pioneering the idea [2]. SModelS currently contains the largest database of
implemented analyses in the community, with over $100$ run I and II available
analyses. With its latest version, SModelS can be interfaced with the full
likelihoods provided by the ATLAS experiment [3]. Another notable tool that
implements the Simplified Model approach is DarkCAST, a recasting tool for
dark photon searches [4].

On the other hand is the Fast Simulation approach, which requires to perform a
Monte-Carlo simulation of the hard scattering process first, then drag the
generated events through some form of detector modelling and lastly pass them
through an event analysis framework. There are many available event analysis
frameworks, some of them with native detector modelling. A non-exhaustive list
of these include MadAnalysis 5 [7, 8, 9, 10, 11, 12], which can be used with a
detector simulation or transfer functions and contains around 40 run I and II
analyses; CheckMATE [13, 14], which contains over 50 run I and II analyses,
and was recently extended to support long-lived particle searches in addition
to prompt particle searches; ColliderBit [15], which includes its own detector
modelling in the form fast 4-vector smearing (Buckfast), and has a database of
around 40 analyses, from runs I and II111ColliderBit is part of the GAMBIT
framework, but it will soon be released as a standalone package; and Rivet
[17, 18], which has currently “only” 30 BSM analyses, but boasts a huge
library of SM analyses, upwards of 800. In fact, it has been shown that SM
inclusive measurements also set strong constraints on new physics, and the
tool Contur has been developed for this purpose [19, 20, 21, 22, 23]. A
comparison of the performance of these tools for a specific CMS search for
supersymmetric (SUSY) particles [24] was performed [25] and strong agreement
was seen across all included tools.

##  3 Dark Matter

Searches for DM are a strong component of the physics agenda of the LHC
experiments, but the strongest constraints on DM models arise often from
direct and indirect detection, as well as from the relic density (RD) of DM,
and thus many tools have been developed over the years for that purpose222A
thorough review of many DM tools can be found in [26]. In recent times, there
is a tendency towards the development of model-independent tools, able to
constrain a large variety of BSM models. In addition, many existing and new
tools attempt a comprehensive approach to constraints of DM, intending to
simultaneously provide constraints from RD, direct and indirect detection. The
most employed comprehensive tools developed by the community are DarkSUSY 6
[27, 28], recently expanded to compute the RD beyond the standard kinetic
equilibrium; MicrOMEGAs [29, 30, 31, 32, 33, 34, 35], whose most recent
upgrade included computations of freeze-in RD; MadDM [36, 37, 38, 39], part of
the MadGraph toolchain; and SuperISO Relic [40, 41, 42], which is widely used
for SUSY models.

Even though comprehensive tools can cover all aspects of DM constraints, some
more dedicated tools are better suited to compute a single constraint with
higher accuracy. In the case of constraints from direct detection (DD), the
front runner is DDCalc [43, 44], which contains the largest collection of DD
experiments to date, and can be used for complete as well as effective DM
models, with RG evolution provided by DirectDM [45, 46]. In the indirect
direction frontier the amount of available tools is too large to recount here.
A short list of selected tools include Capt’n General [47], which imposes
constraints on DM capture in the Sun, nulike [48, 49], which calculates
constraints from the DM annihilation to neutrinos, and gamLike [43], which
deals with DM annihilation to $\gamma$-rays. Lastly one can impose constraints
on DM models from cosmological observations. In particular, the energy
injected into the Cosmic Microwave Background from DM annihilations in the
early Universe can be computed with DarkAges [50], part of the ExoCLASS branch
of CLASS [51].

##  4 Inference and Sampling

It has long been the goal of many particle theorists to assess the validity of
BSM models in light of data. However, the vast quantities of data produced by
the experiments, and the increasingly large parameter spaces of the models,
makes this a difficult task. The common practices of overlaying exclusing
regions and sampling the parameter space using random or grid samplers are
insufficient and much information is lost along the way. Composite likelihoods
that keep all the information about experimental data, smart sampling
strategies to explore large and intricate parameter spaces [52], and a
rigorous treatment of statistics are critical features that significantly
improve and simplify this task [53]. And for this purpose many global fitting
tools have been developed that include those features and that are targeted,
though not exclusively, to BSM models.

Bayesian inference has been a staple in astrophysics for many years, and there
are many public tools developed for that purpose, but in the field of BSM
phenomenology the tendency has been to develop private fitters for specific
purposes, e.g. flavour fits, neutrino oscillations, SUSY, etc. Only in recent
years various public inference tools have been released. The first “wave” of
publicly available global fitters, such as Fittino [54], SuperBayes [55] or
EasyScanHEP [56], were developed to fit various SUSY models to LHC data, and
some are still in use today.

The revolution in the field of BSM inference came with the advent of model-
independent tools. The front-runner for SM fits is HEPFit [57], built upon a
Bayesian Markov chain Monte Carlo, which includes a vast library of
electroweak precision, Higgs and flavour observables, and can also work with
BSM models [58, 59]. Nevertheless, the indisputable leader of BSM inference
studies is the GAMBIT tool [16, 60]. GAMBIT contains a huge library of
likelihoods and observables for collider physics [15], dark matter [43],
flavour [61], precision observables [62], neutrino physics [63] and cosmology
[64], as well as an extensive database of models and sampling strategies, both
frequentist and Bayesian [65]. GAMBIT has been employed in studies of a large
variety of BSM modes, including SUSY models [66, 67, 68], Higgs-portal DM
models [69, 70, 44], EFT DM models [71], axion and axion-like particle models
[72, 73], models with right-handed neutrinos [63], fits of neutral flavour
anomalies [74] and cosmological fits of neutrino masses [75].

##  5 Lagrangians and Model Building

In addition to the software developed for the computation of observables,
there is a small collection of tools that go one step further, and attempt to
bridge the gap between the full Lagrangian of the BSM theory and the required
amplitudes and cross sections. Some of these tools, such as FeynRules [76,
77], SARAH [78, 79, 80] and LanHEP [81, 82, 83], compute the Feynman rules
directly from the Lagrangian, while others, such as FeynArts/FormCalc [84, 85,
86] or CalcHEP [87], use those to compute the required amplitudes. Although
widely used, the main disadvantage of most of these tools is that they are
written in Mathematica, which is a proprietary software framework and it is
notably slower than other interpreted or compiled languages. Very recently a
new tool has been released, MARTY [88], a fully open-source framework that
uses symbolic computations in C++ to obtain amplitudes and cross sections from
the Lagrangian for any BSM theory. Finally, the last step in the BSM toolchain
is to automatically connect the amplitudes, cross sections, etc, obtained by
the Lagrangian-level tools, with the calculation of likelihoods and
observables. This step is most often done by hand, by manually embedding the
code or expressions outputted by the above tools where needed. Fortunately, an
automatic integration of this step can now be achieved with GUM [89]. By
interfacing with SARAH and FeynRules, GUM extracts all information about
particles, parameters and expressions, auto-generates code in GAMBIT, and uses
it to impose DM, collider and vacuum stability constraints. Using both GUM and
GAMBIT one can start with a Lagrangian and end with a complete likelihood
analysis without writing a single line of code.

##  6 Summary

In this conference article I have introduced a selection of tools for BSM
physics and described their most recent developments. The state of the
communication between the theory and experimental communities is currently at
its best, and that improves the quality and relevance of the results from
recasting tools. There is a community trend towards the development of model-
independent tools, as it is the case for some of the DM comprehensive tools,
as well as the leading global fitting frameworks, HEPFit and GAMBIT. And
lastly, with MARTY and GUM we have managed to complete the BSM toolchain, from
Lagrangian to likelihoods, using only open-source software.

## References

  * [1] LHC Reinterpretation Forum collaboration, W. Abdallah et al., _Reinterpretation of LHC Results for New Physics: Status and Recommendations after Run 2_ , _SciPost Phys._ 9 (2020) 022, [2003.07868].
  * [2] S. Kraml, S. Kulkarni, U. Laa, A. Lessa, W. Magerl, D. Proschofsky-Spindler et al., _SModelS: a tool for interpreting simplified-model results from the LHC and its application to supersymmetry_ , _Eur. Phys. J. C_ 74 (2014) 2868, [1312.4175].
  * [3] G. Alguero, S. Kraml and W. Waltenberger, _A SModelS interface for pyhf likelihoods_ , _Comput. Phys. Commun._ 264 (2021) 107909, [2009.01809].
  * [4] P. Ilten, Y. Soreq, M. Williams and W. Xue, _Serendipity in dark photon searches_ , _JHEP_ 06 (2018) 004, [1801.04847].
  * [5] T. Sjostrand, S. Mrenna and P. Z. Skands, _A Brief Introduction to PYTHIA 8.1_ , _Comput. Phys. Commun._ 178 (2008) 852–867, [0710.3820].
  * [6] DELPHES 3 collaboration, J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens et al., _DELPHES 3, A modular framework for fast simulation of a generic collider experiment_ , _JHEP_ 02 (2014) 057, [1307.6346].
  * [7] E. Conte, B. Fuks and G. Serret, _MadAnalysis 5, A User-Friendly Framework for Collider Phenomenology_ , _Comput. Phys. Commun._ 184 (2013) 222–256, [1206.1599].
  * [8] E. Conte, B. Dumont, B. Fuks and C. Wymant, _Designing and recasting LHC analyses with MadAnalysis 5_ , _Eur. Phys. J. C_ 74 (2014) 3103, [1405.3982].
  * [9] B. Dumont, B. Fuks, S. Kraml, S. Bein, G. Chalons, E. Conte et al., _Toward a public analysis database for LHC new physics searches using MADANALYSIS 5_ , _Eur. Phys. J. C_ 75 (2015) 56, [1407.3278].
  * [10] E. Conte and B. Fuks, _Confronting new physics theories to LHC data with MADANALYSIS 5_ , _Int. J. Mod. Phys. A_ 33 (2018) 1830027, [1808.00480].
  * [11] J. Y. Araz, M. Frank and B. Fuks, _Reinterpreting the results of the LHC with MadAnalysis 5: uncertainties and higher-luminosity estimates_ , _Eur. Phys. J. C_ 80 (2020) 531, [1910.11418].
  * [12] J. Y. Araz, B. Fuks and G. Polykratis, _Simplified fast detector simulation in MADANALYSIS 5_ , _Eur. Phys. J. C_ 81 (2021) 329, [2006.09387].
  * [13] M. Drees, H. Dreiner, D. Schmeier, J. Tattersall and J. S. Kim, _CheckMATE: Confronting your Favourite New Physics Model with LHC Data_ , _Comput. Phys. Commun._ 187 (2015) 227–265, [1312.2591].
  * [14] D. Dercks, N. Desai, J. S. Kim, K. Rolbiecki, J. Tattersall and T. Weber, _CheckMATE 2: From the model to the limit_ , _Comput. Phys. Commun._ 221 (2017) 383–418, [1611.09856].
  * [15] GAMBIT collaboration, C. Balázs et al., _ColliderBit: a GAMBIT module for the calculation of high-energy collider observables and likelihoods_ , _Eur. Phys. J. C_ 77 (2017) 795, [1705.07919].
  * [16] GAMBIT collaboration, P. Athron et al., _GAMBIT: The Global and Modular Beyond-the-Standard-Model Inference Tool_ , _Eur. Phys. J. C_ 77 (2017) 784, [1705.07908].
  * [17] A. Buckley, J. Butterworth, D. Grellscheid, H. Hoeth, L. Lonnblad, J. Monk et al., _Rivet user manual_ , _Comput. Phys. Commun._ 184 (2013) 2803–2819, [1003.0694].
  * [18] C. Bierlich et al., _Robust Independent Validation of Experiment and Theory: Rivet version 3_ , _SciPost Phys._ 8 (2020) 026, [1912.05451].
  * [19] J. M. Butterworth, D. Grellscheid, M. Krämer, B. Sarrazin and D. Yallup, _Constraining new physics with collider measurements of Standard Model signatures_ , _JHEP_ 03 (2017) 078, [1606.05296].
  * [20] S. Amrith, J. M. Butterworth, F. F. Deppisch, W. Liu, A. Varma and D. Yallup, _LHC Constraints on a $B-L$ Gauge Model using Contur_, _JHEP_ 05 (2019) 154, [1811.11452].
  * [21] A. Buckley, J. M. Butterworth, L. Corpe, D. Huang and P. Sun, _New sensitivity of current LHC measurements to vector-like quarks_ , _SciPost Phys._ 9 (2020) 069, [2006.07172].
  * [22] J. M. Butterworth, M. Habedank, P. Pani and A. Vaitkus, _A study of collider signatures for two Higgs doublet models with a Pseudoscalar mediator to Dark Matter_ , _SciPost Phys. Core_ 4 (2021) 003, [2009.02220].
  * [23] A. Buckley et al., _Testing new physics models with global comparisons to collider measurements: the Contur toolkit_ , _SciPost Phys. Core_ 4 (2021) 013, [2102.04377].
  * [24] CMS collaboration, A. M. Sirunyan et al., _Search for new physics in events with two soft oppositely charged leptons and missing transverse momentum in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _Phys. Lett. B_ 782 (2018) 440–467, [1801.01846].
  * [25] G. Brooijmans et al., _Les Houches 2019 Physics at TeV Colliders: New Physics Working Group Report_ , in _11th Les Houches Workshop on Physics at TeV Colliders: PhysTeV Les Houches_ , 2, 2020, 2002.12220.
  * [26] C. Arina, _Review on Dark Matter Tools_ , _PoS_ TOOLS2020 (2021) 001, [2012.09462].
  * [27] P. Gondolo, J. Edsjo, P. Ullio, L. Bergstrom, M. Schelke and E. A. Baltz, _DarkSUSY: Computing supersymmetric dark matter properties numerically_ , _JCAP_ 07 (2004) 008, [astro-ph/0406204].
  * [28] T. Bringmann, J. Edsjö, P. Gondolo, P. Ullio and L. Bergström, _DarkSUSY 6 : An Advanced Tool to Compute Dark Matter Properties Numerically_ , _JCAP_ 07 (2018) 033, [1802.03399].
  * [29] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _MicrOMEGAs: A Program for calculating the relic density in the MSSM_ , _Comput. Phys. Commun._ 149 (2002) 103–120, [hep-ph/0112278].
  * [30] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _micrOMEGAs: Version 1.3_ , _Comput. Phys. Commun._ 174 (2006) 577–604, [hep-ph/0405253].
  * [31] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _MicrOMEGAs 2.0: A Program to calculate the relic density of dark matter in a generic model_ , _Comput. Phys. Commun._ 176 (2007) 367–382, [hep-ph/0607059].
  * [32] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _Dark matter direct detection rate in a generic model with micrOMEGAs 2.2_ , _Comput. Phys. Commun._ 180 (2009) 747–767, [0803.2360].
  * [33] G. Belanger, F. Boudjema, P. Brun, A. Pukhov, S. Rosier-Lees, P. Salati et al., _Indirect search for dark matter with micrOMEGAs2.4_ , _Comput. Phys. Commun._ 182 (2011) 842–856, [1004.1092].
  * [34] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _micrOMEGAs $\\_$3: A program for calculating dark matter observables_, _Comput. Phys. Commun._ 185 (2014) 960–985, [1305.0237].
  * [35] G. Bélanger, F. Boudjema, A. Pukhov and A. Semenov, _micrOMEGAs4.1: two dark matter candidates_ , _Comput. Phys. Commun._ 192 (2015) 322–329, [1407.6129].
  * [36] M. Backovic, K. Kong and M. McCaskey, _MadDM v.1.0: Computation of Dark Matter Relic Abundance Using MadGraph5_ , _Physics of the Dark Universe_ 5-6 (2014) 18–28, [1308.4955].
  * [37] M. Backović, A. Martini, O. Mattelaer, K. Kong and G. Mohlabeng, _Direct Detection of Dark Matter with MadDM v.2.0_ , _Phys. Dark Univ._ 9-10 (5, 2015) 37–50, [1505.04190].
  * [38] F. Ambrogi, C. Arina, M. Backovic, J. Heisig, F. Maltoni, L. Mantani et al., _MadDM v.3.0: a Comprehensive Tool for Dark Matter Studies_ , _Phys. Dark Univ._ 24 (2019) 100249, [1804.00044].
  * [39] C. Arina, J. Heisig, F. Maltoni, L. Mantani, D. Massaro, O. Mattelaer et al., _Studying dark matter with MadDM 3.1: a short user guide_ , _PoS_ TOOLS2020 (2021) 009, [2012.09016].
  * [40] A. Arbey and F. Mahmoudi, _SuperIso Relic: A Program for calculating relic density and flavor physics observables in Supersymmetry_ , _Comput. Phys. Commun._ 181 (2010) 1277–1292, [0906.0369].
  * [41] A. Arbey and F. Mahmoudi, _SuperIso Relic v3.0: A program for calculating relic density and flavour physics observables: Extension to NMSSM_ , _Comput. Phys. Commun._ 182 (2011) 1582–1583.
  * [42] A. Arbey, F. Mahmoudi and G. Robbins, _SuperIso Relic v4: A program for calculating dark matter and flavour physics observables in Supersymmetry_ , _Comput. Phys. Commun._ 239 (2019) 238–264, [1806.11489].
  * [43] GAMBIT Dark Matter Workgroup collaboration, T. Bringmann et al., _DarkBit: A GAMBIT module for computing dark matter observables and likelihoods_ , _Eur. Phys. J. C_ 77 (2017) 831, [1705.07920].
  * [44] GAMBIT collaboration, P. Athron et al., _Global analyses of Higgs portal singlet dark matter models using GAMBIT_ , _Eur. Phys. J. C_ 79 (2019) 38, [1808.10465].
  * [45] F. Bishara, J. Brod, B. Grinstein and J. Zupan, _DirectDM: a tool for dark matter direct detection_ , 1708.02678.
  * [46] J. Brod, A. Gootjes-Dreesbach, M. Tammaro and J. Zupan, _Effective Field Theory for Dark Matter Direct Detection up to Dimension Seven_ , _JHEP_ 10 (2018) 065, [1710.10218].
  * [47] N. A. Kozar, A. Caddell, L. Fraser-Leach, P. Scott and A. C. Vincent, _Capt’n General: A generalized stellar dark matter capture and heat transport code_ , in _Tools for High Energy Physics and Cosmology_ , 5, 2021, 2105.06810.
  * [48] IceCube collaboration, P. Scott et al., _Use of event-level neutrino telescope data in global fits for theories of new physics_ , _JCAP_ 11 (2012) 057, [1207.0810].
  * [49] IceCube collaboration, M. G. Aartsen et al., _Improved limits on dark matter annihilation in the Sun with the 79-string IceCube detector and implications for supersymmetry_ , _JCAP_ 04 (2016) 022, [1601.00653].
  * [50] P. Stöcker, M. Krämer, J. Lesgourgues and V. Poulin, _Exotic energy injection with ExoCLASS: Application to the Higgs portal model and evaporating black holes_ , _JCAP_ 03 (2018) 018, [1801.01871].
  * [51] D. Blas, J. Lesgourgues and T. Tram, _The Cosmic Linear Anisotropy Solving System (CLASS). Part II: Approximation schemes_ , _JCAP_ 2011 (July, 2011) 034, [1104.2933].
  * [52] DarkMachines High Dimensional Sampling Group collaboration, C. Balázs et al., _A comparison of optimisation algorithms for high-dimensional particle and astrophysics applications_ , _JHEP_ 05 (2021) 108, [2101.04525].
  * [53] S. S. AbdusSalam et al., _Simple and statistically sound strategies for analysing physical theories_ , 2012.09874.
  * [54] P. Bechtle et al., _Constrained Supersymmetry after two years of LHC data: a global view with Fittino_ , _JHEP_ 06 (2012) 098, [1204.4199].
  * [55] R. Ruiz de Austri, R. Trotta and L. Roszkowski, _A Markov chain Monte Carlo analysis of the CMSSM_ , _JHEP_ 05 (2006) 002, [hep-ph/0602028].
  * [56] C. Han, K.-i. Hikasa, L. Wu, J. M. Yang and Y. Zhang, _Status of CMSSM in light of current LHC Run-2 and LUX data_ , _Phys. Lett. B_ 769 (2017) 470–476, [1612.02296].
  * [57] J. De Blas et al., _HEPfit : a code for the combination of indirect and direct constraints on high energy physics models_, _Eur. Phys. J. C_ 80 (2020) 456, [1910.14012].
  * [58] C.-W. Chiang, G. Cottin and O. Eberhardt, _Global fits in the Georgi-Machacek model_ , _Phys. Rev. D_ 99 (2019) 015001, [1807.10660].
  * [59] D. Chowdhury and O. Eberhardt, _Update of Global Two-Higgs-Doublet Model Fits_ , _JHEP_ 05 (2018) 161, [1711.02095].
  * [60] A. Kvellestad, P. Scott and M. White, _GAMBIT and its Application in the Search for Physics Beyond the Standard Model_ , 1912.04079.
  * [61] GAMBIT Flavour Workgroup collaboration, F. U. Bernlochner et al., _FlavBit: A GAMBIT module for computing flavour observables and likelihoods_ , _Eur. Phys. J. C_ 77 (2017) 786, [1705.07933].
  * [62] GAMBIT Models Workgroup collaboration, P. Athron et al., _SpecBit, DecayBit and PrecisionBit: GAMBIT modules for computing mass spectra, particle decay rates and precision observables_ , _Eur. Phys. J. C_ 78 (2018) 22, [1705.07936].
  * [63] M. Chrzaszcz, M. Drewes, T. E. Gonzalo, J. Harz, S. Krishnamurthy and C. Weniger, _A frequentist analysis of three right-handed neutrinos with GAMBIT_ , _Eur. Phys. J. C_ 80 (2020) 569, [1908.02302].
  * [64] GAMBIT Cosmology Workgroup collaboration, J. J. Renk et al., _CosmoBit: A GAMBIT module for computing cosmological observables and likelihoods_ , _JCAP_ 02 (2021) 022, [2009.03286].
  * [65] GAMBIT collaboration, G. D. Martinez, J. McKay, B. Farmer, P. Scott, E. Roebber, A. Putze et al., _Comparison of statistical sampling methods with ScannerBit, the GAMBIT scanning module_ , _Eur. Phys. J. C_ 77 (2017) 761, [1705.07959].
  * [66] GAMBIT collaboration, P. Athron et al., _Global fits of GUT-scale SUSY models with GAMBIT_ , _Eur. Phys. J. C_ 77 (2017) 824, [1705.07935].
  * [67] GAMBIT collaboration, P. Athron et al., _A global fit of the MSSM with GAMBIT_ , _Eur. Phys. J. C_ 77 (2017) 879, [1705.07917].
  * [68] GAMBIT collaboration, P. Athron et al., _Combined collider constraints on neutralinos and charginos_ , _Eur. Phys. J. C_ 79 (2019) 395, [1809.02097].
  * [69] GAMBIT collaboration, P. Athron et al., _Status of the scalar singlet dark matter model_ , _Eur. Phys. J. C_ 77 (2017) 568, [1705.07931].
  * [70] P. Athron, J. M. Cornell, F. Kahlhoefer, J. Mckay, P. Scott and S. Wild, _Impact of vacuum stability, perturbativity and XENON1T on global fits of $\mathbb{Z}_{2}$ and $\mathbb{Z}_{3}$ scalar singlet dark matter_, _Eur. Phys. J. C_ 78 (2018) 830, [1806.11281].
  * [71] GAMBIT collaboration, P. Athron et al., _Thermal WIMPs and the Scale of New Physics: Global Fits of Dirac Dark Matter Effective Field Theories_ , 2106.02056.
  * [72] S. Hoof, F. Kahlhoefer, P. Scott, C. Weniger and M. White, _Axion global fits with Peccei-Quinn symmetry breaking before inflation using GAMBIT_ , _JHEP_ 03 (2019) 191, [1810.07192].
  * [73] P. Athron et al., _Global fits of axion-like particles to XENON1T and astrophysical data_ , _JHEP_ 05 (2021) 159, [2007.05517].
  * [74] J. Bhom, M. Chrzaszcz, F. Mahmoudi, M. T. Prim, P. Scott and M. White, _A model-independent analysis of $b\to s\mu^{+}\mu^{-}$ transitions with GAMBIT’s FlavBit_, 2006.03489.
  * [75] GAMBIT Cosmology Workgroup collaboration, P. Stöcker et al., _Strengthening the bound on the mass of the lightest neutrino with terrestrial and cosmological experiments_ , _Phys. Rev. D_ 103 (2021) 123508, [2009.03287].
  * [76] N. D. Christensen and C. Duhr, _FeynRules - Feynman rules made easy_ , _Comput. Phys. Commun._ 180 (2009) 1614–1641, [0806.4194].
  * [77] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, _FeynRules 2.0 - A complete toolbox for tree-level phenomenology_ , _Comput. Phys. Commun._ 185 (2014) 2250–2300, [1310.1921].
  * [78] F. Staub, _SARAH_ , 0806.0538.
  * [79] F. Staub, _SARAH 3.2: Dirac Gauginos, UFO output, and more_ , _Comput. Phys. Commun._ 184 (2013) 1792–1809, [1207.0906].
  * [80] F. Staub, _SARAH 4 : A tool for (not only SUSY) model builders_ , _Comput. Phys. Commun._ 185 (2014) 1773–1790, [1309.7223].
  * [81] A. V. Semenov, _LanHEP: A Package for automatic generation of Feynman rules in gauge models_ , hep-ph/9608488.
  * [82] A. V. Semenov, _LanHEP: A Package for automatic generation of Feynman rules in field theory. Version 2.0_ , hep-ph/0208011.
  * [83] A. Semenov, _LanHEP: A Package for the automatic generation of Feynman rules in field theory. Version 3.0_ , _Comput. Phys. Commun._ 180 (2009) 431–454, [0805.0555].
  * [84] R. Mertig, M. Bohm and A. Denner, _FEYN CALC: Computer algebraic calculation of Feynman amplitudes_ , _Comput. Phys. Commun._ 64 (1991) 345–359.
  * [85] J. Kublbeck, M. Bohm and A. Denner, _Feyn Arts: Computer Algebraic Generation of Feynman Graphs and Amplitudes_ , _Comput. Phys. Commun._ 60 (1990) 165–180.
  * [86] T. Hahn, _Generating Feynman diagrams and amplitudes with FeynArts 3_ , _Comput. Phys. Commun._ 140 (2001) 418–431, [hep-ph/0012260].
  * [87] A. Belyaev, N. D. Christensen and A. Pukhov, _CalcHEP 3.4 for collider physics within and beyond the Standard Model_ , _Comput. Phys. Commun._ 184 (2013) 1729–1769, [1207.6082].
  * [88] G. Uhlrich, F. Mahmoudi and A. Arbey, _MARTY - Modern ARtificial Theoretical phYsicist A C++ framework automating symbolic calculations Beyond the Standard Model_ , _Comput. Phys. Commun._ 264 (2021) 107928, [2011.02478].
  * [89] S. Bloor, T. E. Gonzalo, P. Scott, C. Chang, A. Raklev, J. E. Camargo-Molina et al., _The GAMBIT Universal Model Machine: from Lagrangians to Likelihoods_ , 2107.00030.