# COMPOSING PARTIAL DIFFERENTIAL EQUATIONS
## WITH PHYSICS-AWARE NEURAL NETWORKS

**Anonymous authors**
Paper under double-blind review

ABSTRACT

We introduce a compositional physics-aware neural network (FINN) for learning
spatiotemporal advection-diffusion processes. FINN implements a new way of
combining the learning abilities of artificial neural networks with physical and
structural knowledge from numerical simulation by modeling the constituents of
partial differential equations (PDEs) in a compositional manner. Results on both
one- and two-dimensional PDEs (Burger’s, diffusion-sorption, diffusion-reaction,
Allen-Cahn) demonstrate FINN’s superior modeling accuracy and excellent outof-distribution generalization ability beyond initial and boundary conditions. With
only one tenth of the number of parameters on average, FINN outperforms pure
machine learning and other state-of-the-art physics-aware models in all cases—
often even by multiple orders of magnitude. Moreover, FINN outperforms a calibrated physical model when approximating sparse real-world data in a diffusionsorption scenario, confirming its generalization abilities and showing explanatory
potential by revealing the unknown retardation factor of the observed process.

1 INTRODUCTION

Artificial neural networks (ANNs) are considered universal function approximators (Cybenko,
1989). Their effective learning ability, however, greatly depends on domain and task-specific prestructuring and methodological modifications referred to as inductive biases (Battaglia et al., 2018).
Typically, inductive biases limit the space of possible models by reducing the opportunities for computational shortcuts, which can lead to erroneous implications derived from a potentially limited
dataset (overfitting). The recently evolving field of physics-informed machine learning employs
physical knowledge as inductive bias providing huge generalization advantages in contrast to pure
machine learning (ML) in physical domains (Raissi et al., 2019). While numerous approaches have
been introduced to augment ANNs with physical knowledge, these methods either do not allow the
incorporation of explicitly defined physical equations (Long et al., 2018; Seo et al., 2019; Guen
& Thome, 2020; Li et al., 2020a; Sitzmann et al., 2020) or cannot generalize to other initial and
boundary conditions than those encountered during training (Raissi et al., 2019).

In this work, we present the finite volume neural network (FINN) model—a physics-aware neural
network structure adhering to the idea of spatial and temporal discretization in numerical simulation.
FINN consists of multiple neural network modules that interact in a distributed, compositional manner (Battaglia et al., 2018; Lake et al., 2017; Lake, 2019). The modules are designed to account for
specific parts of advection-diffusion equations, a class of partial differential equations (PDEs). This
modularization allows to combine two advantages that are not yet met by state-of-the-art models:
the explicit incorporation of physical knowledge and the generalization over initial and boundary
conditions. To the best of our knowledge, FINN’s ability to adjust to different initial and boundary
conditions and to explicitly learn constitutive relationships and reaction terms is unique, yielding
excellent out-of-distribution generalization. The core contributions of this work are:

-  Introduction of FINN, a physics-aware neural network model, explicitly designed to generalize over initial and boundary conditions, demonstrating excellent generalization ability.

-  Evaluation of state-of-the-art pure ML and physics-aware models, contrasted to FINN on
one- and two-dimensional benchmarks, demonstrating benefits of explicit model design.

-  Application of FINN to a real-world contamination-diffusion problem, verifying its applicability to real, spatially and temporally constrained training data.


-----

2 RELATED WORK

**Non-physics-aware ANN architectures** Pure ML models that are designed for spatiotemporal
data processing can be separated into temporal convolution (TCN, Kalchbrenner et al., 2016) and
recurrent neural networks. While the former perform convolutions over space and time, representatives of the latter, e.g., convolutional LSTM (ConvLSTM, Shi et al., 2015) or DISTANA (Karlbauer
et al., 2019), aggregate spatial neighbor information to further process the temporal component with
recurrent units. Since pure ML models do not adhere to physical principles, they require enormous
amounts of training data and parameters in order to approximate a desired physical process; but still
are not guaranteed to behave consistently outside the regime of the training data.

**Physics-aware ANN architectures** When designed to satisfy physical equations, ANNs are reported to have greater robustness in terms of physical plausibility. For example, the physicsinformed neural network (PINN, Raissi et al., 2019) consists of an MLP that satisfies an explicitly
defined PDE with specific initial and boundary conditions, using automatic differentiation. However, the remarkable results beyond the time steps encountered during training are limited to the
very particular PDE and its conditions. A trained PINN can actually not be applied to other different
initial and boundary conditions, which limits its applicability in real-world scenarios.

Other methods by Long et al. (PDENet, 2018), Guen & Thome (PhyDNet, 2020), or Sitzmann et al.
(SIREN, 2020) learn the first n derivatives to achieve a physically plausible behavior, but lack the
option to include physical equations. The same limitation applies when operators are learned to approximate PDEs (Li et al., 2020a;b), or when physics-aware graph neural networks are applied (Seo
et al., 2019). Yin et al. (APHYNITY, 2020) propose to approximate equations with an appropriate
physical model and to augment the result by an ANN, preventing the ANN to approximate a distinct
part within the physical model. For more comparison to related work, please refer to subsection A.1.

In summary, none of the above methods can explicitly learn particular constitutive relationships or
reaction terms while simultaneously generalizing beyond different initial and boundary conditions.

3 FINITE VOLUME NEURAL NETWORK (FINN)

**Problem formulation** Here, we focus on modeling spatiotemporal physical processes. Specifically, we consider systems governed by advection-diffusion type equations (Smolarkiewicz, 1983):
_∂u_

(1)

_∂t_ [=][ D][(][u][)] _[∂]∂x[2][u][2][ −]_ _[v][(][u][)]_ _[∂u]∂x_ [+][ q][(][u][)][,]

where u is the quantity of interest, t is time, x is the spatial coordinate, D is the diffusion coefficient, v is the advection velocity, and q is the source/sink term. Eq. 1 can be partitioned into three
parts: the storage term, the flux terms, and the source/sink term. The storage term _[∂u]∂t_ [describes]

the change of the quantity u over time. The flux terms are the advective flux v(u) _[∂u]∂x_ [and the diffu-]

sive flux D(u) _[∂]∂x[2][u][2][ . Both calculate the amount of][ u][ exchanged between neighboring positions. The]_

source/sink term q(u) describes the generation or elimination of the quantity u. Eq. 1 is a general
form of PDE with up to second order spatial derivative, but it has a wide range of applicability due
to the flexibility of defining D(u), v(u), and q(u) as different functions of u, as is shown by the
numerical experiments in this work.

The finite volume method (FVM, Moukalled et al., 2016) discretizes a simulation domain into control volumes (i = 1, . . ., Nx), where exchange fluxes are calculated using a surface integral (Riley
et al., 2009). In order to match this structure in FINN, we introduce two different kernels, which are
(spatially) casted across the discretized control volumes: the flux kernel, modeling the flux terms,
and the state kernel, modeling the source/sink term as well as the storage term. The overall FINN
architecture is shown in Figure 1.

**Flux kernel** The flux kernel F approximates the surface integral for each control volume i with
boundary Ω by a composition of multiple subkernels fj, each representing the flux through a discretized surface element j:


_Nsi_

_fj_
_j=1_ _≈_

X


_D(u)_ _[∂][2][u]_

_∂x[2][ −]_ _[v][(][u][)]_ _[∂u]∂x_



_· ˆn dΓ,_ (2)


_Fi =_


_ω⊆Ω_


-----

_ϕN+i−_
_ϕNi+_

_∂[2]ui_

_∂x[2][ −]_


_ϕorNi−_
_ϕNi+_

_∂ui_

_∂x_


_ϕD_

_D (ui)_


_R(ϕA)_

_v (ui)_


Φψ

_q (ui)_


_∂ui_

_∂t_ [=]


_Fi = {fi−, fi+}_

_∂u∂ti(t)_ ODE solver _L_ _uˆ[(]i[t][+1)], u[(]i[t][+1)]_

uˆ[(]i[t][+1)] 

+

+ _fi−_ _Si_ _fi+_ +

_R_ _×_ _×_ _×_ _×_ _R_

_ϕA_ _ϕNi−_ _ϕD_ Φψ _ϕD_ _ϕNi+_ _ϕA_

_· · ·_ _u[(]i−[t][)]1_ _u[(]i[t][)]_ _u[(]i+1[t][)]_ _· · ·_


Figure 1: Flux and state kernels in FINN for the one-dimensional case (left); detailed assignment of
the individual modules with their contribution to Eq. 1 (right). Read lines indicate gradient flow.

where Nsi is the number of discrete surface elements of control volume i, ω is a continuous surface
element (a subset of Ω), fj are subkernels (which are realized as feedforward network modules), and
_nˆ is the unit normal vector pointing outwards of ω. In our exemplary one-dimensional arrangement,_
two subkernels fi− and fi+ (see Figure 1) contain the modules ϕD, ϕA, and ϕN (ϕNi− or ϕNi+ for
_fi_ and fi+, respectively). The module ϕ is a linear layer with the purpose to learn the numerical
_−_ _N_
FVM stencil. More specifically, the module ϕNi− takes the inputs ui and its neighbor ui−1, while
_ϕNi+ takes the inputs ui and ui+1, and output the approximation of the spatial derivative_ _[∂u]∂x_ [. This]

signifies that the weight of ϕ should amount to [ 1, 1] with respect to [ui, ui 1] and [ui, ui+1]
_N_ _−_ _−_
(i.e. simple difference between neighboring control volumes) in ideal one-dimensional problems.

The other modules ϕ and ϕ are responsible for advective and diffusive fluxes, respectively. Both
_A_ _D_
modules receive only ui as input. Moreover, the module R applies an upwind differencing scheme,
which prevents numerical instability in the first order spatial derivative calculation (Versteeg &
Malalasekera, 1995). It ensures that calculation of the advective fluxes using ϕ is performed only
_A_
on one control volume surface (either left or right, and not both at the same time), whereas calculation of the diffusive fluxes using ϕ is performed on both surfaces. This means, for the advective
_D_
flux calculation, the summation of the numerical stencil from both fi and fi+ will only lead to
_−_

[ 1, 1], which corresponds to either [ui, ui 1] when ϕ _> 0 or [ui, ui+1] when ϕ_ _< 0 (i.e. only_
_−_ _−_ _A_ _A_
first order spatial derivative). For the diffusive flux calculation, the summation of the numerical
stencil will lead to the classical one-dimensional numerical Laplacian with [1, −2, 1] corresponding
to [ui 1, ui, ui+1], that is, the second order spatial derivative, because the calculation of ϕ is per_−_ _D_
formed on both surfaces (see subsection A.3 for more details). Accordingly, module R generates
the inductive bias to make ϕ only approximate the advective, and ϕ the diffusive flux.
_A_ _D_

For the advective flux calculation, if the velocity v is a function of u, the module ϕ approximates
_A_
this dependence using a feedforward neural network, i.e. ϕ (u) _v(u). Otherwise, ϕ_ is a scalar
_A_ _≈_ _A_
value ϕ _v, which can be also set as a learnable parameter. The output of ϕ_ is activated by the
following case-sensitive ReLUA ≡ _R:_ _A_

(ϕ ) = ReLU(ϕA), on fi−, (3)
_R_ _A_ − ReLU(−ϕA), on fi+,

and then multiplied with the output of ϕN (i.e. _[∂u]∂x_ [) to obtain the corresponding advective flux at the]

surface. Similarly for the diffusive fluxes, if the diffusion coefficient D depends on u, the module
_ϕ_ takes u as an input to learn the function ϕ (u) _D(u) using a feedforward neural network._
_D_ _D_ _≈_
Otherwise, ϕ is a scalar value ϕ _D, which can also be set as a learnable parameter. The output_
of ϕ is not activated by the moduleD _D ≡_, but is directly multiplied with the output of ϕ to obtain
_D_ _R_ _N_
the corresponding diffusive flux at the surface. Consequently, integration of subkernels fj in Eq. 2
and substitution of Eq. 3 leads to the following representation:

_ϕ_ (ui, ui 1)[ϕ _ϕ_ ]fi + ϕ (ui, ui+1)[ϕ ]fi+ _,_ if ϕ _> 0,_
_i =_ _N_ _−_ _D −_ _A_ _−_ _N_ _D_ _A_ (4)
_F_ _ϕ_ (ui, ui 1)[ϕ ]fi + ϕ (ui, ui+1)[ϕ _ϕ_ ]fi+ _,_ otherwise.
 _N_ _−_ _D_ _−_ _N_ _D −_ _A_

**Boundary conditions** A means of applying boundary conditions in the model is essential when
solving PDEs. Currently available models mostly adopt a convolutional structure for modelling
spatiotemporal processes. However, a convolutional structure only allows a constant value to be


-----

padded at the domain boundaries (e.g. zero-padding or mirror-padding), which is only appropriate
for the implementation of Dirichlet or periodic boundary condition types. Other types of frequently
used boundary conditions are Neumann and Cauchy. They are defined as a derivative of the quantity
of interest, and hence cannot be easily implemented in convolutional models. With certain predefined boundary condition types, the flux kernels at the boundaries are adjusted accordingly to
allow for straightforward boundary condition implementation. For Dirichlet boundary condition, a
constant value u = ub is set as the input ui 1 (for the flux kernel fi ) or ui+1 (for fi+) at the
_−_ _−_
corresponding boundary. For Neumann boundary condition ν, the output of the flux kernel fi− or
_fi+ at the corresponding boundary is set to be equal to ν. With Cauchy boundary condition, the_
solution-dependent derivative is calculated and set as ui 1 or ui+1 at the corresponding boundary.
_−_

**State kernel** The state kernel S calculates the source/sink and storage terms of Eq. 1. The
source/sink (if required) is learned using the module Φψ _q(u), which takes u as input. The_
_≈_
storage, _[∂u]∂t_ [, is then calculated using the output of the flux kernel and module][ Φ][ψ][ of the state kernel:]

_Si = Fi + Φψ(ui) ≈_ _[∂u]∂t [i]_ _[.]_ (5)

By doing so, the PDE in Eq. 1 is now reduced to a system of coupled ordinary differential equations
(ODEs), which are functions of ui, ui−1, ui+1, and t. Thus, the solutions of the coupled ODE system
can be found using a numerical integration over time. Since first order explicit approaches, such as
the Euler method (Butcher, 2008), suffer from numerical instability (Courant et al., 1967; Isaacson
& Keller, 1994), we employ the neural ordinary differential equation method (Neural ODE, Chen
et al., 2018) to reduce numerical instability via the Runge-Kutta adaptive time-stepping strategy.
When using the Neural ODE method, the predicted values of u at time t are fed back into the
network as direct inputs to predict u at time t + 1. Therefore, the entire training is performed in
closed-loop, improving stability and accuracy of the prediction compared to networks trained with
teacher forcing, i.e. one-step-ahead prediction (Praditia et al., 2020). The weight update is realized
by applying backpropagation through time (indicated by red arrows in Figure 1). In short, FINN
takes only the initial condition u at time t = 0 and propagates the dynamics forward.

4 EXPERIMENTS, RESULTS & DISCUSSION

4.1 SYNTHETIC DATASET

To demonstrate FINN’s performance in comparison to other models, four different equations are
considered as applications. First, Burger’s equation (Basdevant et al., 1986) is chosen as a challenging function, as it is a non-linear PDE with v(u) = u that could lead to a shock in the solution u(x, t). Second, the diffusion-sorption equation (Nowak & Guthke, 2016) is selected with the
non-linear retardation factor R(u) as coefficient for the storage term, which contains a singularity
_R(u) →∞_ for u → 0 due to the parameter choice. Third, the two-dimensional Fitzhugh-Nagumo
equation (Klaasen & Troy, 1984) as candidate for a diffusion-reaction equation (Turing, 1952) is
selected which is challenging because it consists of two non-linearly coupled PDEs to solve two
main unknowns: the activator u1 and the inhibitor u2. Fourth, the Allen-Cahn equation with a cubic reaction term was chosen, leading to multiple jumps in the solution u(x, t). Details on all four
equations, data generation and architecture designs can be found in subsection A.4, subsection A.5,
subsection A.6, subsection A.7 of the appendix, respectively.

For each problem, three different datasets are generated (by conventional numerical simulation):
_train, used to train the models, in-distribution test (in-dis-test), being the train data simulated with_
a longer time span to test the models’ generalization ability (extrapolation), and out-of-distribution
test (out-dis-test). Out-dis-test data are used to test a trained ML model under conditions that are far
away from training conditions, not only in terms of querying outputs for unseen inputs. Instead, out_dis-test data query outputs with regards to changes not captured by the inputs. These are changes_
that the ML tool per definition cannot be made aware of during training. In this work, they are
represented by data generated with different initial or boundary condition, to test the generalization
ability of the models outside the training distributions. FINN is trained and compared with both
spatiotemporal deep learning models such as TCN, ConvLSTM, DISTANA and physics-aware neural network models such as PINN and PhyDNet. All models are trained with ten different random


-----

Table 1: Comparison of MSE and according standard deviation scores across ten repetitions between
different deep learning (above dashed line) and physics-aware neural network (below dashed line)
methods on the different equations. Best results are reported in bold.

Dataset


Eqn. Model Params _Train_ _In-dis-test_ _Out-dis-test_


TCN 38 500 (1.6 ± 3.4) × 10[−][1] (1.8 ± 3.1) × 10[−][1] (1.6 ± 3.3) × 10[−][1]

ConvLSTM 13 200 (6.8 ± 9.9) × 10[−][2] (1.2 ± 1.1) × 10[−][1] (7.3 ± 9.5) × 10[−][2]

DISTANA 25 126 (1.8 ± 1.0) × 10[−][4] (4.0 ± 3.1) × 10[−][3] (1.5 ± 1.6) × 10[−][3]

Burger PINN 3 021 (5.1 ± 0.5) × 10[−][4] (5.0 ± 8.4) × 10[−][3] - 

PhyDNet 37 718 (7.2 ± 2.2) × 10[−][5] (1.8 ± 1.6) × 10[−][1] (4.5 ± 2.5) × 10[−][2]

FINN **421** (2.8 ± 2.9) × 10[−][6] (2.5 ± 3.1) × 10[−][6] (2.8 ± 2.9) × 10[−][6]

TCN 3 834 (9.7 ± 13.5) × 10[−][2] (1.2 ± 1.7) × 10[−][1] (1.1 ± 1.4) × 10[−][1]

ConvLSTM 3 960 (3.2 ± 2.9) × 10[−][2] (3.0 ± 2.3) × 10[−][2] (5.8 ± 4.0) × 10[−][2]

DISTANA 3 739 (4.6 ± 2.5) × 10[−][5] (2.4 ± 2.5) × 10[−][3] (4.6 ± 4.6) × 10[−][3]

sorption PINN 3 042 (4.7 8.4) 10[−][5] (4.1 8.7) 10[−][3] - 
Diffusion- _±_ _×_ _±_ _×_

PhyDNet 37 815 (3.5 ± 1.7) × 10[−][5] (9.1 ± 15.4) × 10[−][3] (1.7 ± 0.9) × 10[−][2]

FINN **528** (4.7 ± 4.9) × 10[−][5] (1.3 ± 1.3) × 10[−][4] (4.1 ± 4.0) × 10[−][5]

TCN 31 734 (1.4 ± 0.9) × 10[−][2] (4.7 ± 2.1) × 10[−][1] (1.5 ± 0.8) × 10[−][1]

ConvLSTM 24 440 (8.7 ± 21.3) × 10[−][3] (9.3 ± 4.9) × 10[−][2] (1.5 ± 1.3) × 10[−][2]

DISTANA 75 629 (4.0 ± 3.4) × 10[−][3] (1.8 ± 0.6) × 10[−][1] (1.3 ± 0.9) × 10[−][2]

reaction PINN 3 062 (2.7 1.8) 10[−][4] (7.0 6.3) 10[−][2] - 
Diffusion- _±_ _×_ _±_ _×_

PhyDNet 185 589 (7.5 ± 0.9) × 10[−][5] (7.8 ± 1.8) × 10[−][2] (3.5 ± 1.3) × 10[−][2]

FINN **882** (1.3 ± 0.3) × 10[−][4] (2.1 ± 0.5) × 10[−][3] (6.1 ± 0.3) × 10[−][3]

seeds using PyTorch’s default weight initialization, and mean and standard deviation of the prediction errors are summarized in Table 1 for train, in-dis-test and out-dis-test. The details of each run
are reported in the appendix. It is noteworthy that PINN cannot be tested on the out-dis-test dataset,
since PINN assumes that the unknown variable u is an explicit function of x and t, and hence, when
the initial or boundary condition is changed, the function will also be different and no longer valid.


**Burger’s** The predictions of the best trained model of each method for the in-dis-test and the out_dis-test data are shown in Figure 2 and Figure 3, respectively. Both TCN and ConvLSTM fail to_
produce reasonable predictions, but qualitatively the other models manage to capture the shape of
the data sufficiently, even towards the end of the in-dis-test period, where closed loop prediction has
been applied for 380 time steps (after 20 steps of teacher forcing, see subsection A.4 for details).
DISTANA, PINN, and FINN stand out in particular, but FINN produces more consistent and accurate predictions, evidenced by the mean value of the prediction errors. When tested against data
generated with a different initial condition (out-dis-test), all models except for TCN and PhyDNet

Figure 2: Plots of Burger’s data (red) and in-dis-test (blue) prediction using different models. The
plots in the first row show the solution over x and t (the red lines mark the transition from train to
_in-dis-test), the second row visualizes the best model’s solution distributed in x at t = 2._


-----

Figure 3: Plots of Burger’s data (red) and prediction (blue) of out-dis-test data using different models. The plots in the first row show the solution over x and t, the second row visualizes the best
model’s solution over x at t = 1.

perform well. However, FINN still outperforms the other models with a significantly lower prediction error. The advective velocity learned by FINN’s module ϕ is shown in Figure 8 (left) and
_A_
verifies that it successfully learned the advective velocity to be described by an identity function.

**Diffusion-sorption** The predictions of the best trained model of each method for the concentration
_u from the diffusion-sorption equation are shown in Figure 4 and Figure 5. TCN and ConvLSTM are_
shown to perform poorly even on the train data, evidenced by the high mean value of the prediction
errors. On in-dis-test data, all models successfully produce relatively accurate predictions. However, when tested against different boundary conditions (out-dis-test), only FINN is able to capture
the modifications and generalize well. The other models are shown to still overfit to the different
boundary condition used in the train data (as detailed in subsection A.5). The retardation factor
_R(u) learned by FINN’s ϕ_ module is shown in Figure 8 (second from left). The plot shows that
_D_
the module learned the Freundlich retardation factor with reasonable accuracy.

**Diffusion-reaction** The predictions of the best trained model of each method for activator u1 in
the diffusion-reaction equation are shown in Figure 6 and Figure 7. TCN is again shown to fail to
learn sufficiently from the train data. On in-dis-test data, DISTANA and the physics-aware models
all predict with reasonable accuracy. When tested against data with different initial condition (out_dis-test), however, DISTANA and PhyDNet produce predictions with lower accuracy, and we find_
that FINN is the only model producing relatively low prediction errors. The reaction functions
learned by FINN’s Φψ module are shown in Figure 8. The plots show that the module successfully
learned the Fitzhugh-Nagumo reaction function, both for the activator and inhibitor.

**Discussion** Even with a high number of parameters, the prediction errors obtained using the pure
ML methods (TCN, ConvLSTM, DISTANA) are worse compared to the physics-aware models, as

Figure 4: Plots of diffusion-sorption’s dissolved concentration (red) and in-dis-test prediction (blue)
using different models. The first row shows the solution over x and t (red lines mark the transition
from train to in-dis-test), the second row visualizes the solution distributed in x at t = 10 000.


-----

Figure 5: Plots of diffusion-sorption’s dissolved concentration data (red) and prediction (blue) of
_out-dis-test data using different models. The plots in the first row show the solution over x and t,_
the second row visualizes the solution distributed in x at t = 10 000.

shown in Table 1. As a physics-aware model, PhyDNet also possesses a high number of parameters.
However, most of the parameters are allocated to the data-driven part (i.e. ConvLSTM branch),
compared to the physics-aware part (i.e. PhyCell branch). In contrast to the other pure ML methods,
DISTANA predicts with higher accuracy. This could act as an evidence that appropriate structural
design of a neural network is as essential as physical knowledge to regularize the model training.

Among the physics-aware models, PINN and PhyDNet lie on different extremes. On one side, PINN
requires complete knowledge of the modelled system in form of the equation. This means that all
the functions, such as the advective velocity in the Burger’s equation, the retardation factor in the
diffusion-sorption equation, and the reaction function in the diffusion-reaction equation have to be
pre-defined in order to train the model. This leaves less room for learning from data and could be
error-prone if the designer’s assumption is incorrect. On the other side, PhyDNet relies more heavily
on the data-driven part and, therefore, could overfit the train data. This can be shown by the fact that
PhyDNet reaches the lowest training errors for the diffusion-sorption and diffusion-reaction equation predictions compared to the other models, but its performance significantly deteriorates when
applied to in- and out-dis-test data. Our proposed model, FINN, lies somewhere in the middle, compromising between putting physical knowledge into the model and leaving room for learning from
data. As a consequence, FINN shows excellent generalization ability. It significantly outperforms
the other models up to multiple orders of magnitude, especially on out-dis-test data, when tested
with different initial and boundary conditions; which is considered a particularly challenging task
for ML models. Furthermore, the structure of FINN allows the extractions of learned functions such
as the advective velocity, retardation factor, and reaction functions, showing good interpretability of
the model.

We also show that FINN properly handles the provided numerical boundary condition, as evidenced
when applied to the test data that is generated with a different left boundary condition value, visual
Figure 6: Plots of diffusion-reaction’s activator data u1 (red) and extrapolated prediction (blue) using
different models. The plots in the first row show the solution distributed over x and y at t = 50, and
the plots in the second row show the solution distributed in x at y = 0 and t = 50.


-----

Figure 7: Plots of diffusion-reaction’s activator data (red) and prediction (blue) of unseen dataset
using different models. The plots in the first row show the solution distributed over x and y at
_t = 10, and the plots in the second row show the solution distributed in x at y = 0 and t = 10._

Figure 8: Plots of the learned functions (blue) as a function of u compared to the data (red) for
Burger’s (left) and diffusion-sorption (second from left). The learned activator u1 and inhibitor u2
reaction functions (right plots of the pairs in center and right) in the diffusion-reaction equation are
contrasted to the corresponding ground truth (left plots of the pairs).

ized in Figure 5. Here, the test data is generated with a Dirichlet boundary condition u(0, t) = 0.7,
which is lower than the value used in the train data, u(0, t) = 1.0. However, FINN is the only model
that appropriately processes this different boundary condition value so that the prediction fits the test
data nicely. The other models overestimate their prediction by consistently showing a tendency to
still fit the prediction to the higher boundary condition value encountered during training.

Even though the spatial resolution used for the synthetic data generation is relatively coarse, leading
to sparse data availability, FINN generalizes well. PINN, on the other hand, slightly suffers from
the low resolution of the train data, although it still shows reasonable performance for the three test
cases. Nevertheless, we conducted an experiment showing that PINN performs slightly better and
more consistently when trained on higher resolution data (see appendix, subsection A.10), albeit
still worse than FINN on coarse data. Therefore, we conclude that FINN is also applicable to real
observation data that are often available only in low resolution, and/or in limited quantity. We
demonstrate this further in the next section, when we apply FINN to real experimental data. FINN’s
generalization ability is superior to PINN, due to the fact that it is not possible to apply a trained
PINN model to predict data with different initial or boundary condition.

In terms of interpretability, FINN allows the extraction of functions learned by its dedicated modules.
These learned functions can be compared with the real functions that generated the data (at least in
the synthetic data case); examples are shown in Figure 8. The learned functions are the main datadriven discovery part of FINN and can also be used as a physical plausibility check. PhyDNet also
comprises of a physics-aware and a data-driven part. However, it is difficult, if not impossible, to
infer the learned equation from the model. Furthermore, the data-driven part of PhyDNet does not
possess comparable interpretability, and could lead to overfitting as discussed earlier.

4.2 EXPERIMENTAL DATASET

Real observation data are often available in limited amount, and they only provide partial information about the system of interest. To demonstrate how FINN can cope with real observation data,
we use experimental data for the diffusion-sorption problem. The experimental data are collected
from three different core samples that are taken from the same geographical area: #1,#2, and #2B


-----

Figure 9: Breakthrough curve prediction of FINN (blue line) during training using data from core
sample #2 (left), during testing using data from core sample #1 (second from left) and total concentration profile prediction using data from core sample #2B (second from right). The predictions are
compared with the experimental data (red circles) and the results obtained using the physical model
(orange dashed line). The learned retardation factor R(u) is shown in the right-most plot.

(see subsection A.12 in the appendix for details). The objective of this experiment is to learn the
retardation factor function from one of the core samples that concurrently applies to all the other
samples. For this particular purpose, we only implement FINN since the other models have no
means of learning the retardation factor explicitly. Here, we use the module ϕ to learn the retarda_D_
tion factor function, and we assume that the diffusion coefficient values of all the core samples are
known. FINN is trained with the breakthrough curve of u, which is the dissolved concentration only
at x = L|0≤t≤tend (i.e. only 55 data points).

FINN reaches a higher accuracy for the training with core sample #2, with MSE = 4.84 × 10[−][4]
compared to a calibrated physical model from Nowak & Guthke (2016) with MSE = 1.06 × 10[−][3],
because the latter has to assume a specific function R(u) with a few parameters for calibration. Our
learned retardation factor is then applied and tested to core samples #1 and #2B. Figure 9 shows
that FINN’s prediction accuracy is competitive compared to the calibrated physical model. For
core sample #1, FINN’s prediction has an accuracy of 1.37 × 10[−][3] compared to the physical model
that underestimates the breakthough curve (i.e. concentration profile) with MSE = 2.50 × 10[−][3].
Core sample #2B has significantly longer length than the other samples, and therefore a no-flow
Neumann boundary condition was assumed at the bottom of the core. Because there is no breakthrough curve data available for this specific sample, we compare the prediction against the so-called
total concentration profile u(x, tend) at the end of the experiment. FINN produced a prediction
with an accuracy of 1.16 × 10[−][3], whereas the physical model overestimates the concentration with
MSE = 2.73 × 10[−][3]. To briefly summarize, FINN is able to learn the retardation factor from a
sparse experimental dataset and apply it to other core samples with similar soil properties with reasonable accuracy, even applying to a different boundary condition type.

5 CONCLUSION

Spatiotemporal dynamics often can be described by means of advection-diffusion type equations,
such as Burger’s, diffusion-sorption, or diffusion-reaction equations. When modeling those dynamics with ANNs, large efforts must be taken to prevent the models from overfitting (given the model
is able to learn the problem at all). The incorporation of physical knowledge can be seen as regularization, yielding robust predictions beyond the training data distribution.

With FINN, we have introduced a modular, physics-aware neural network with excellent generalization abilities beyond different initial and boundary conditions, when contrasted to pure ML models
and other physics-aware models. FINN is able to model and extract unknown constituents of differential equations, allowing high interpretability and an assessment of the plausibility of the model’s
out-of-distribution predictions. As next steps we seek to apply FINN beyond second order spatial derivatives, improve its scalability to large datasets, and make it applicable to heterogeneously
distributed data (i.e. represented as graphs) by modifying the module ϕ to approximate variable
_N_
and location-specific stencils (for more details on limitations of FINN, the reader is referred to subsection A.2). Another promising future direction is the application of FINN to real-world weather
data.


-----

AUTHOR CONTRIBUTIONS

Hidden for anonymization reasons.

ACKNOWLEDGMENTS

We thank the reviewers for suggesting the evaluation of polynomial fitting and least squares regression as additional baselines as well as modeling another challenging equation (we chose Allen-Cahn
with its cubic term) to demonstrate the wide applicability of the method, which we have added
to the appendix. The results confirm our claims and further suggestions helped us to increase the
comprehensibility and clarity of the method. Additional parts are hidden for anonymization reasons

REPRODUCIBILITY STATEMENT

Code and data that are used for this paper can be found in the repository ***** [hidden for
_anonymization reasons. Instead, the supplementary material of this submission contains an anony-_
_mous version of the repository’s README.md file along with the according data and model scripts.]_

REFERENCES

Y. Bar-Sinai, S. Hoyer, J. Hickey, and M.P. Brenner. Learning data-driven discretizations for partial
differential equations. Proceedings of the National Academy of Sciences, 116(31):15344–15349,
2019.

C. Basdevant, M. Deville, P. Haldenwang, J.M. Lacroix, J. Ouazzani, R. Peyret, P. Orlandi, and A.T.
Patera. Spectral and finite difference solutions of the burgers equation. Computers & Fluids, 14
(1):23–41, 1986. doi: https://doi.org/10.1016/0045-7930(86)90036-8.

P.W. Battaglia, J.B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, et al. Relational inductive biases, deep learning, and
graph networks. arXiv preprint arXiv:1806.01261, 2018.

J.C. Butcher. _Numerical Methods for Ordinary Differential Equations._ Wiley, 2008. ISBN
9780470753750.

R.T.Q. Chen, Y. Rubanova, J. Bettencourt, and D.K. Duvenaud. Neural ordinary differential equations. arXiv preprint arXiv:1806.07366, 2018.

R. Courant, K. Friedrichs, and H. Lewy. On the partial difference equations of mathematical physics.
_IBM Journal of Research and Development, 11(2):215–234, 1967. doi: 10.1147/rd.112.0215._

G. Cybenko. Approximation by superpositions of a sigmoidal function. Mathematics of control,
_signals and systems, 2(4):303–314, 1989._

B. Fornberg. Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of
_computation, 51(184):699–706, 1988._

V.L. Guen and N. Thome. Disentangling physical dynamics from unknown factors for unsupervised
video prediction. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern
_Recognition, pp. 11474–11484, 2020._

E. Isaacson and H.B. Keller. Analysis of Numerical Methods. Dover Books on Mathematics. Dover
Publications, 1994. ISBN 9780486680293.

N. Kalchbrenner, L. Espeholt, K. Simonyan, A. Oord, A. Graves, and K. Kavukcuoglu. Neural
machine translation in linear time. arXiv preprint arXiv:1610.10099, 2016.

M. Karlbauer, S. Otte, H.P.A. Lensch, T. Scholten, V. Wulfmeyer, and M.V. Butz. A distributed
neural network architecture for robust non-linear spatio-temporal prediction. _arXiv preprint_
_arXiv:1912.11141, 2019._


-----

G.A. Klaasen and W.C. Troy. Stationary wave solutions of a system of reaction-diffusion equations
derived from the fitzhugh–nagumo equations. SIAM Journal on Applied Mathematics, 44(1):
96–110, 1984. doi: 10.1137/0144008.

D. Kochkov, J.A. Smith, A. Alieva, Q. Wang, M.P. Brenner, and S. Hoyer. Machine learning–
accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences, 118
(21), 2021.

B.M. Lake. Compositional generalization through meta sequence-to-sequence learning. In Advances
_in Neural Information Processing Systems 32, pp. 9791–9801. 2019._

B.M. Lake, T.D. Ullman, J.B. Tenenbaum, and S.J. Gershman. Building machines that learn and
think like people. Behavioral and Brain Sciences, 2017. doi: 10.1017/S0140525X16001837.

Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar. Fourier neural operator for parametric partial differential equations. _arXiv preprint_
_arXiv:2010.08895, 2020a._

Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar. Neural operator: Graph kernel network for partial differential equations. arXiv preprint
_arXiv:2003.03485, 2020b._

Z. Long, Y. Lu, X. Ma, and B. Dong. Pde-net: Learning pdes from data. In International Conference
_on Machine Learning, pp. 3208–3216. PMLR, 2018._

F. Moukalled, L. Mangani, and M. Darwish. The Finite Volume Method in Computational Fluid
_Dynamics. Springer, 1 edition, 2016. doi: 10.1007/978-3-319-16874-6._

W. Nowak and A. Guthke. Entropy-based experimental design for optimal model discrimination in
the geosciences. Entropy, 18(11), 2016. doi: 10.3390/e18110409.

T. Praditia, T. Walser, S. Oladyshkin, and W. Nowak. Improving thermochemical energy storage dynamics forecast with physics-inspired neural network architecture. Energies, 13(15):3873, 2020.
doi: 10.3390/en13153873.

T. Praditia, M. Karlbauer, S. Otte, S. Oladyshkin, M.V. Butz, and W. Nowak. Finite volume neural
network: Modeling subsurface contaminant transport. arXiv preprint arXiv:2104.06010, 2021.

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, 378:686–707, 2019. doi: 10.1016/j.jcp.2018.10.045.

K.F. Riley, M.P. Hobson, and S.J. Bence. Mathematical Methods for Physics and Engineering.
Cambridge University Press, 2009. ISBN 9780521139878.

S. Seo, C. Meng, and Y. Liu. Physics-aware difference graph networks for sparsely-observed dynamics. In International Conference on Learning Representations, 2019.

X. Shi, Z. Chen, H. Wang, D.Y. Yeung, W.K. Wong, and W.C. Woo. Convolutional lstm network: A
machine learning approach for precipitation nowcasting. arXiv preprint arXiv:1506.04214, 2015.

V. Sitzmann, J. Martel, A Bergman, D. Lindell, and G. Wetzstein. Implicit neural representations
with periodic activation functions. Advances in Neural Information Processing Systems, 33, 2020.

P.K. Smolarkiewicz. A simple positive definite advection scheme with small implicit diffusion.
_Monthly Weather Review, 111(3):479 – 486, 1983. doi: 10.1175/1520-0493._

A. Turing. The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society
_B, 237:37–72, 1952._

H.K. Versteeg and W. Malalasekera. An Introduction to Computational Fluid Dynamics: The Finite
_Volume Method. New York, 1995. ISBN 9780470235157._


-----

Y. Yin, V.L. Guen, J. Dona, I. Ayed, E. de B´ezenac, N. Thome, and P. Gallinari. Augmenting physical
models with deep networks for complex dynamics forecasting. arXiv preprint arXiv:2010.04456,
2020.

J. Zhuang, D. Kochkov, Y. Bar-Sinai, M.P. Brenner, and S. Hoyer. Learned discretizations for passive
scalar advection in a two-dimensional turbulent flow. Physical Review Fluids, 6(6):064605, 2021.

A APPENDIX

This appendix is structured as follows: The detailed differences with the benchmark models used
in this work are presented in subsection A.1. Then, limitations of FINN are discussed briefly in
subsection A.2. Detailed numerical derivation of the first and second order spatial derivative is
shown in subsection A.3. Additional details on the equations used as well as model specifications
and in-depth results are presented in subsection A.4 (Burger’s), subsection A.5 (diffusion-sorption),
subsection A.6 (diffusion-reaction), and subsection A.7 (Allen-Cahn). Additional polynomial-fitting
baselines to account for the equations are reported in subsection A.8, including an ablation where
we replace the neural network modules of FINN by polynomials. A robustness analysis of FINN can
be found in subsection A.9, where our method is evaluated on noisy data from all equations. Results
of training PINN with high-resolution data are reported in subsection A.10, whereas results of the
original PhyDNet architecture are outlined in subsection A.11. The soil parameters and simulation
domain used in the diffusion-sorption laboratory experiment are presented in subsection A.12.

A.1 DISTINCTION TO RELATED WORK

**Pure ML models** Originally, FINN was inspired by DISTANA, a pure ML model proposed by
Karlbauer et al. (2019). While DISTANA has large similarities to Shi et al. (2015)’s ConvLSTM,
it propagates lateral information through an additional latent map—instead of reading lateral information from the input map directly, as done in ConvLSTM—and transforms that lateral information
by means of a user-defined combination of arbitrary layers. Additionally, the processing of the twopoint flux approximation in FINN using the lateral information, is more akin to DISTANA, albeit
motivated and augmented by physical knowledge.

Accordingly, the lateral information flow is guaranteed to behave in a physically plausible manner.
That is, quantity can either be locally generated by a source (increased), locally absorbed by a
sink (decreased), or spatially distributed (move to neighboring cells). Terminology separates the
spatially distribution of quantity into diffusion and advection. Diffusion describes the equalization
of quantity from high to low concentration levels, whereas advection is defined as the bulk motion
of a large group of particles/atoms caused by external forces. FINN ensures the conservation of
laterally propagating quantity, i.e. what flows from left to right will effectively cause a decrease
left and an increase right. Pure ML models (without physical constraints) cannot be guaranteed to
adhere to these fundamental rules.

**PINN** Since ML models guided by physical knowledge not only behave empirically plausible
but also require fewer parameters and generalize to much broader extent, numerous approaches
have been proposed recently to combine artificial neural networks with physical knowledge. Raissi
et al. (2019) introduced the physics-informed neural network (PINN), a concrete and outstanding
model that explicitly learns a provided PDE. As a result, PINN mimics e.g. Burger’s equation (see
equation 7) by learning the quantity function u(x, t) for defined initial and boundary conditions with
a feedforward network. The partial derivatives are implicitly realized by automatic differentiation
of the feedforward network (representing u) with respect to the input variable x or t. The learned
neural network thus satisfies the constraints defined by the partial derivatives.

This method has the advantage that it explicitly provides a solution of the desired function u(x, t)
and, correspondingly, predictions can be generated for an arbitrary combination of input values,
circumventing the need for simulating the entire domain with e.g. a carefully chosen simulation
step size. In contrast, FINN does not learn the explicit function u(x, t) defined for particular initial
and boundary conditions, but approximates the distinct components of the function. These components are combined as suggested by the physical equation to result in a compositional model that is
more universally applicable. That is, in stark contrast to PINN, the compositional function learned


-----

by FINN can be applied to varying initial and boundary conditions, since the learned individual
components provide the same functionality as the corresponding components of the PDE when processed by a numerical solver. Applying the numerical solver, i.e. the finite volume method (FVM),
however, requires either complete knowledge of the equation or careful calibration of the unknowns
by choosing equations from a library of possible solutions and tuning the parameters. FINN’s data
driven component can reveal unknown relations, such as the retardation factor of a function, without
the need for subjective prior assumptions.

**Learning derivatives via convolution** An alternative and much addressed approach of implanting
physical plausibilty into artificial neural networks is to implicitly learn the first n derivatives using
appropriately designed convolution kernels (Long et al., 2018; Li et al., 2020a;b; Guen & Thome,
2020; Yin et al., 2020; Sitzmann et al., 2020). These methods exploit the link that most PDEs, such
as u(t, x, y, . . . ), can be reformulated as

_∂u_ _t, x, y, . . ., u, [∂u]_ _∂u_ _._ (6)

_∂t_ [=][ F] _∂x_ _[, ∂u]y [,]_ _∂x∂y [, ∂]∂x[2][u][2][, . . .]_
 


When learning partial derivatives up to order n and setting irrelevant features to zero, these methods
have, in principle, the capacity to represent most PDEs. However, the degrees of freedom in these
methods are still very high and can fail to safely guide the ML algorithm. FINN accounts for the first
and second derivative by learning an according stencil in the module ϕ and combining this stencil
_N_
with the case-sensitive ReLU module R, which allows a precise control of the information flow
resulting from the first and second derivative. More importantly, convolutional structure only allows
implementation of Dirichlet and periodic boundary condition (by means of zero- or mirror-padding),
and is not appropriate for implementation of other boundary condition types.

**Learning ODE coefficients** The group around Brenner and Hoyer (Bar-Sinai et al., 2019;
Kochkov et al., 2021; Zhuang et al., 2021) follow another line of research about learning the coefficients of ordinary differential equations (ODEs); note that PDEs can be transformed into a set
of coupled ODEs by means of polynomial expansion or spectral differentiation as shown in BarSinai et al. (2019). The physical constraint in these works is mostly realized by using the temporal
derivative as loss.

While this approach share similarities with our method by incorporating physically motivated inductive bias into the model, it uses this bias mainly to improve interpolation from coarse to finer
grid resolutions and thus to accelerate simulation. Our work focuses on discovering unknown relationships/laws (or re-discovering laws in the case of the synthetic examples), such as the advective
velocity in the Burger’s example, the retardation factor function in the diffusion-sorption example,
and the reaction function in the diffusion-reaction example. Additionally, the works by Brenner and
Hoyer employ a convolutional structure, which is only applicable to Dirichlet or periodic boundary
conditions, and it suffers from a slight instability during training when the training data trajectory
is unrolled for a longer period. In contrast, FINN employs the flux kernel, calculated at all control
volume surfaces, which enables the implementation and discovery of various boundary conditions.
Furthermore and in contrast to Brenner and Hoyer, FINN employs the Neural ODE method as the
time integrator to reduce numerical instability during training with long time series. However and in
accord with this line of research, FINN is also able to generalize well when trained with a relatively
sparse dataset (coarse resolution), reducing the computational burden.

**Numerical ODE solvers** Traditional numerical solvers for ordinary differential equations (ODEs)
can be seen as the pure-physics contrast to FINN. However, these do not have a learning capacity to
reveal unknown dependencies or functions from data. Also, as shown in (Yin et al., 2020), a pure
Neural ODE without physical inductive bias does not reach the same level of accuracy as physicsaware neural network models. This is underlined by our field experiment where FINN reached lower
errors on a real-world dataset compared to a conventional FVM model.

A.2 LIMITATIONS OF FINN

While the largest limitation of our current method can be seen in the capacity to only represent
first and second order spatial derivatives, this is an issue that we will address in follow up work.
Still, FINN can already be applied to a very wide range of problems as most equations in fact only


-----

depend on up to the second spatial derivative. So far, FINN is only applicable to homogeneously
distributed data—we intend to extend it to heterogeneous data from graphs. Although we have
successfully applied FINN to 2D diffusion reaction data, the training time is considerable. While,
according to our observations, this appears to be a common issue for physics-aware neural networks,
the implementation of custom-convolution layers could widen this bottleneck with today’s hardwareaccelerated computation of convolution operations.

A.3 LEARNING THE NUMERICAL STENCIL

Semantically, the ϕ module learns the numerical stencil, that is the geometrical arrangement of a
_N_
group of neighboring cells to approximate the derivative numerically, effectively learning the first
spatial derivativemodule, respectively.[∂u]∂x [from both][ [][u][i][−][1][, u][i][]][ and][ [][u][i][, u][i][+1][]][, which are the inputs to the][ ϕ][N][i][−] [and][ ϕ][N][i][+]

The lateral information flowing from ui 1 and ui+1 toward ui is controlled by the ϕ (advective
_−_ _A_
flux, i.e. bulk motion of many particles/atoms that can either move to the left or to the right) and
the ϕ (diffusive flux, i.e. drive of particles/atoms to equilibrium from regions of high to low
_D_
concentration) modules. Since the advective flux can only move either to the left or to the right, it
will be considered only in the left flux kernel (fi−) or in the right flux kernel (fi+), and not both at
the same time. The case-sensitive ReLU module R (Eq. 3) decides on this, by setting the advective
flux in the irrelevant flux kernel to zero (effectively depending on the sign of the output of ϕ ).
_A_
Thus, the advective flux is only considered from either ui 1 or ui+1 to ui, which amounts to the
_−_
first order spatial derivative.

The diffusive flux, on the other hand, can propagate from both sides towards the control volume
of interest ui and, hence, the second order spatial derivative, accounting for the difference between ui 1 and ui+1, has to be applied. In our method, this is realized through the combina_−_
tion of the ϕ and ϕ modules, calculating the diffusive fluxes δ = ϕ (ui 1, ui)ϕ (ui) and
_N_ _D_ _−_ _N_ _−_ _D_
_δ+ = ϕ_ (ui, ui+1)ϕ (ui) inside of the respective left and right flux kernel. The combination of
_N_ _D_
these two deltas in the state kernel ensures the consideration of the diffusive fluxes from left (including ui−1) and from right (including ui+1), resulting in the ability to account for the second order
spatial derivative.

Technically, the first, i.e. [−1, 1], and second, i.e. [1, −2, 1], order spatial differentiation schemes
are common definitions and a derivation can be found, for example, in Fornberg (1988). However,
a quick derivation of the Laplace scheme [1, −2, 1] can be formulated as follows. Define the second
order spatial derivative as
_∂[2]u_
_∂x[2][ ≈]_ [(][∂u/∂x][)][|][i][−]∆[−]x[(][∂u/∂x][)][|][i][+]

with (∂u/∂x)|i− _≈_ (ui−1 − _ui)/∆x and (∂u/∂x)|i+ ≈_ (ui − _ui+1)/∆x. Then_

_∂[2]u_

= _[u][i][−][1][ −]_ [2][u][i][ +][ u][i][+1]

_∂x[2][ ≈]_ [(][u][i][−][1][ −] _[u][i][)]∆[ −]x[2][(][u][i][ −]_ _[u][i][+1][)]_ ∆x[2]

hence the [+1, −2, +1] as coefficients.

A.4 BURGER

The Burger’s equation is commonly employed in various research areas, including fluid mechanics.

**Data** The 1D generalized Burger’s equation is written as

_∂u_

(7)

_∂t_ [=][ −][v][(][u][)] _[∂u]∂x_ [+][ D∂]∂x[2][u][2][,]

where the main unknown is u, the advective velocity is denoted as v(u) which is a function of u
and the diffusion coefficient is D = 0.01/π. In the current work, the advective velocity function
is chosen to be an identity function v(u) = u to reproduce the experiment conducted in the PINN
paper (Raissi et al., 2019).

The simulation domain for the train data is defined with x = [−1, 1], t = [0, 1] and is discretized
with Nx = 49 spatial locations, and Nt = 201 simulation steps. The initial condition is defined as
_u(x, 0) = −_ sin(πx), and the boundary condition is defined as u(−1, t) = u(1, t) = 0.


-----

Figure 10: Prediction mean over ten different trained models (with 95% confidence interval) of the
Burger’s equation at t = 2 for the in-dis-test dataset.

**_In-dis-test data is simulated with x = [−1, 1] and a time span of t = [1, 2] and Nt = 401. Initial_**
condition is taken from the train data at t = 1 and boundary conditions are also similar to the train
data.

The simulation domain for the out-dis-test data is identical with the train data, except for the initial
condition that is defined as u(x, 0) = sin(πx).

**Model architectures** Both TCN and ConvLSTM are designed to have one input neuron, one
hidden layer of size 32, and one output neuron. The lateral and dynamic input- and output sizes
of the DISTANA model are set to one and a hidden layer of size 32 is used. The pure ML models
were trained on the first 150 time steps and validated on the remaining 50 time steps of the train data
(applying early stopping). Also, to prevent the pure ML models from diverging too much in closed
loop, the boundary data are fed into the models as done during teacher forcing. PINN was defined
as a feedforward network with the size of [2, 20, 20, 20, 20, 20, 20, 20, 20, 1] (8 hidden layers, each
contains 20 hidden neurons), as reported in the original work by Raissi et al. (2019). PhyDNet was
defined with the PhyCell containing 32 input dimensions, 7 hidden dimensions, 1 hidden layer, and
the ConvLSTM containing 32 input dimensions, 32 hidden dimensions, 1 hidden layer. For FINN,
the modules ϕ, ϕ, and ϕ were used, with ϕ defined as a feedforward network with the
_N_ _D_ _R_ _A_ _A_
size of [1, 10, 20, 10, 1] that takes u as an input and outputs the advective velocity v(u), and ϕ as
_D_
a learnable scalar that learns the diffusion coefficient D. All models are trained until convergence
using the L-BFGS optimizer, except for PhyDNet, which is trained with the Adam optimizer and a
learning rate of 1 × 10[−][3] due to stability issues when training with the L-BFGS optimizer.

**Additional results** Individual errors of the ten different training runs as well as visualizations for
the train (Table 2), in-dis-test (Table 3 and Figure 10), and out-dis-test (Table 4 and Figure 11)
datasets.

Table 2: Closed loop MSE on the train data from ten different training runs for each model for the
Burger’s equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 4.1 × 10[−][2] 7.9 × 10[−][2] 1.2 × 10[−][4] 6.1 × 10[−][4] 8.0 × 10[−][5] 1.9 × 10[−][6]

2 5.1 × 10[−][2] 1.7 × 10[−][1] 1.8 × 10[−][4] 4.1 × 10[−][4] 4.7 × 10[−][5] 1.5 × 10[−][6]

3 3.6 × 10[−][2] 1.3 × 10[−][2] 6.4 × 10[−][5] 4.9 × 10[−][4] 6.8 × 10[−][5] 2.2 × 10[−][6]

4 4.3 × 10[−][2] 2.7 × 10[−][4] 2.9 × 10[−][4] 5.6 × 10[−][4] 1.2 × 10[−][4] 2.5 × 10[−][6]

5 5.4 × 10[−][2] 2.1 × 10[−][4] 2.2 × 10[−][4] 5.1 × 10[−][4] 6.7 × 10[−][5] 1.6 × 10[−][6]

6 3.9 × 10[−][2] 1.5 × 10[−][3] 1.3 × 10[−][4] 5.2 × 10[−][4] 5.5 × 10[−][5] 9.9 × 10[−][7]

7 6.6 × 10[−][2] 4.1 × 10[−][4] 5.4 × 10[−][5] 5.2 × 10[−][4] 5.1 × 10[−][5] 2.2 × 10[−][6]

8 1.2 × 10[0] 3.6 × 10[−][2] 2.7 × 10[−][4] 4.8 × 10[−][4] 5.6 × 10[−][5] 3.0 × 10[−][6]

9 3.4 × 10[−][2] 5.5 × 10[−][2] 3.9 × 10[−][4] 5.5 × 10[−][4] 9.9 × 10[−][5] 1.1 × 10[−][5]

10 6.8 × 10[−][2] 3.2 × 10[−][1] 1.1 × 10[−][4] 5.0 × 10[−][4] 8.1 × 10[−][5] 4.5 × 10[−][7]


-----

Table 3: Closed loop MSE on in-dis-test data from ten different training runs for each model for the
Burger’s equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 3.0 × 10[−][2] 1.2 × 10[−][1] 2.6 × 10[−][3] 4.1 × 10[−][3] 1.1 × 10[−][1] 1.6 × 10[−][6]

2 7.0 × 10[−][2] 3.7 × 10[−][1] 1.1 × 10[−][3] 1.3 × 10[−][3] 2.5 × 10[−][1] 1.3 × 10[−][6]

3 2.4 × 10[−][2] 2.0 × 10[−][1] 8.9 × 10[−][4] 2.9 × 10[−][3] 3.4 × 10[−][2] 1.9 × 10[−][6]

4 5.1 × 10[−][2] 5.6 × 10[−][3] 6.2 × 10[−][3] 3.7 × 10[−][3] 4.7 × 10[−][1] 2.1 × 10[−][6]

5 5.4 × 10[−][2] 2.6 × 10[−][3] 2.6 × 10[−][3] 2.4 × 10[−][3] 1.1 × 10[−][2] 1.3 × 10[−][6]

6 3.2 × 10[−][2] 1.0 × 10[−][2] 1.2 × 10[−][2] 1.9 × 10[−][3] 3.2 × 10[−][1] 6.8 × 10[−][7]

7 8.8 × 10[−][2] 5.5 × 10[−][3] 1.1 × 10[−][3] 2.1 × 10[−][3] 3.3 × 10[−][1] 1.8 × 10[−][6]

8 1.1 × 10[0] 1.6 × 10[−][1] 4.0 × 10[−][3] 4.0 × 10[−][4] 2.7 × 10[−][1] 2.4 × 10[−][6]

9 3.3 × 10[−][2] 1.9 × 10[−][1] 4.0 × 10[−][3] 3.0 × 10[−][2] 2.6 × 10[−][2] 1.2 × 10[−][5]

10 3.6 × 10[−][1] 9.8 × 10[−][2] 6.0 × 10[−][3] 1.5 × 10[−][3] 8.0 × 10[−][3] 2.8 × 10[−][7]


Table 4: Closed loop MSE on out-dis-test data from ten different training runs for each model for
the Burger’s equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 3.6 × 10[−][2] 7.2 × 10[−][2] 5.9 × 10[−][3] -  5.3 × 10[−][2] 1.9 × 10[−][6]

2 3.4 × 10[−][2] 1.9 × 10[−][1] 4.4 × 10[−][4] -  2.1 × 10[−][2] 1.5 × 10[−][6]

3 3.3 × 10[−][2] 2.9 × 10[−][3] 2.9 × 10[−][3] -  4.5 × 10[−][2] 2.2 × 10[−][6]

4 3.7 × 10[−][2] 8.9 × 10[−][4] 1.4 × 10[−][3] -  2.0 × 10[−][2] 2.5 × 10[−][6]

5 4.2 × 10[−][2] 1.2 × 10[−][3] 1.0 × 10[−][3] -  3.1 × 10[−][2] 1.6 × 10[−][6]

6 2.7 × 10[−][2] 1.4 × 10[−][3] 3.0 × 10[−][4] -  1.1 × 10[−][1] 9.9 × 10[−][7]

7 6.7 × 10[−][2] 3.3 × 10[−][2] 6.2 × 10[−][4] -  2.8 × 10[−][2] 2.2 × 10[−][6]

8 1.1 × 10[0] 4.0 × 10[−][2] 7.0 × 10[−][4] -  4.7 × 10[−][2] 3.0 × 10[−][6]

9 3.0 × 10[−][2] 8.5 × 10[−][2] 7.4 × 10[−][4] -  4.5 × 10[−][2] 1.1 × 10[−][5]

10 1.2 × 10[−][1] 3.1 × 10[−][1] 7.0 × 10[−][4] -  4.5 × 10[−][2] 4.5 × 10[−][7]


Figure 11: Prediction mean over ten different trained models (with 95% confidence interval) of the
Burger’s equation at t = 1 for the out-dis-test data.


-----

A.5 DIFFUSION-SORPTION

The diffusion-sorption equation is another widely applied equation in fluid mechanics. A practical
example of the equation is to model contaminant transport in groundwater. Its retardation factor R
can be modelled using different closed parametric relations known as sorption isotherms that should
be calibrated to observation data.

**Data** The 1D diffusion-sorption equation is written as the following coupled system of equations:


_D_ _∂[2]u_ (8) _∂ut_ (9)

_R(u)_ _∂x[2][,]_ _∂t_ [=][ Dφ∂]∂x[2][u][2][,]


_∂u_

_∂t_ [=]


where the dissolved concentration u and total concentration ut are the main unknowns. The effective
diffusion coefficient is denoted by D = 5 × 10[−][4], the retardation factor is R(u), a function of u,
and the porosity is denoted by φ = 0.29.

In this work, the Freundlich sorption isotherm, see details in (Nowak & Guthke, 2016), was chosen
to define the retardation factor:


_R(u) = 1 + [1][ −]_ _[φ]_


_ρsknf_ _u[n][f][ −][1],_ (10)


where ρs = 2 880 is the bulk density, k = 3.5 _×_ 10[−][4] is the Freundlich’s parameter, and nf = 0.874
is the Freundlich’s exponent.

The simulation domain for the train data is defined with x = [0, 1], t = [0, 2 500] and is discretized
with Nx = 26 spatial locations, and Nt = 501 simulation steps. The initial condition is defined as
_u(x, 0) = 0, and the boundary condition is defined as u(0, t) = 1.0 and u(1, t) = D_ _[∂u]∂x_ [.]

The in-dis-test data was simulated with x = [0, 1] and with the time span of t = [2 500, 10 000] and
_Nt = 2 001. Initial condition is taken from the train data at t = 2 500 and boundary conditions are_
also identical to the train data.

The simulation domain for the out-dis-test data was identical with the in-dis-test data, except for
the boundary condition that was defined as u(0, t) = 0.7.

**Model architectures** **TCN is designed to have two input neurons, four hidden layers of size [4,**
8, 8, 8], and two output neurons. ConvLSTM has two input- and output neurons and one hidden
layer with 16 neurons. The lateral and dynamic input- and output sizes of the DISTANA model
are set to one and two, respectively, while a hidden layer of size 12 is used. The pure ML models
were trained on the first 400 time steps and validated on the remaining 100 time steps of the train
data (applying early stopping). Also, to prevent the pure ML models from diverging too much in
closed loop, the boundary data are fed into the models as done during teacher forcing. PINN was

Table 5: Closed loop MSE on the train data from ten different training runs for each model for the
diffusion-sorption equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 3.6 × 10[−][2] 4.8 × 10[−][5] 2.1 × 10[−][5] 2.2 × 10[−][4] 3.1 × 10[−][5] 1.1 × 10[−][4]

2 2.0 × 10[−][1] 2.3 × 10[−][2] 1.6 × 10[−][5] 2.3 × 10[−][6] 2.5 × 10[−][5] 1.4 × 10[−][5]

3 8.1 × 10[−][2] 4.4 × 10[−][2] 2.7 × 10[−][5] 1.2 × 10[−][6] 3.2 × 10[−][5] 1.1 × 10[−][4]

4 2.0 × 10[−][1] 1.0 × 10[−][2] 9.3 × 10[−][5] 2.1 × 10[−][4] 2.9 × 10[−][5] 4.1 × 10[−][5]

5 1.8 × 10[−][3] 8.7 × 10[−][2] 6.2 × 10[−][5] 2.5 × 10[−][6] 7.7 × 10[−][5] 5.3 × 10[−][6]

6 2.7 × 10[−][4] 5.9 × 10[−][2] 4.0 × 10[−][5] 1.7 × 10[−][6] 4.4 × 10[−][5] 1.9 × 10[−][5]

7 3.2 × 10[−][4] 6.4 × 10[−][2] 7.7 × 10[−][5] 8.6 × 10[−][6] 1.8 × 10[−][5] 1.5 × 10[−][5]

8 1.7 × 10[−][2] 1.4 × 10[−][4] 2.5 × 10[−][5] 2.0 × 10[−][6] 2.4 × 10[−][5] 1.2 × 10[−][4]

9 4.3 × 10[−][1] 2.2 × 10[−][3] 5.9 × 10[−][5] 9.0 × 10[−][6] 2.1 × 10[−][5] 1.2 × 10[−][5]

10 6.6 × 10[−][4] 2.7 × 10[−][2] 3.9 × 10[−][5] 1.1 × 10[−][5] 5.2 × 10[−][5] 1.7 × 10[−][5]


-----

Table 6: Closed loop MSE on in-dis-test data from ten different training runs for each model for the
diffusion-sorption equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 1.9 × 10[−][2] 7.7 × 10[−][3] 6.9 × 10[−][4] 3.0 × 10[−][2] 3.8 × 10[−][3] 3.2 × 10[−][4]

2 2.4 × 10[−][1] 1.0 × 10[−][2] 8.5 × 10[−][5] 8.6 × 10[−][4] 8.4 × 10[−][3] 3.5 × 10[−][5]

3 4.7 × 10[−][2] 2.1 × 10[−][2] 1.7 × 10[−][3] 9.3 × 10[−][5] 2.9 × 10[−][4] 3.2 × 10[−][4]

4 2.4 × 10[−][1] 2.8 × 10[−][2] 3.5 × 10[−][3] 5.2 × 10[−][3] 3.6 × 10[−][4] 1.1 × 10[−][4]

5 3.5 × 10[−][2] 7.9 × 10[−][2] 7.3 × 10[−][3] 7.9 × 10[−][6] 2.8 × 10[−][2] 9.8 × 10[−][6]

6 2.4 × 10[−][2] 5.4 × 10[−][2] 3.1 × 10[−][4] 1.4 × 10[−][5] 4.8 × 10[−][2] 4.9 × 10[−][5]

7 1.4 × 10[−][3] 4.2 × 10[−][2] 6.5 × 10[−][3] 1.8 × 10[−][4] 1.0 × 10[−][3] 3.7 × 10[−][5]

8 3.9 × 10[−][2] 2.2 × 10[−][4] 9.4 × 10[−][4] 6.9 × 10[−][4] 2.5 × 10[−][4] 3.4 × 10[−][4]

9 5.6 × 10[−][1] 1.4 × 10[−][2] 1.8 × 10[−][3] 4.5 × 10[−][3] 2.4 × 10[−][4] 2.8 × 10[−][5]

10 3.4 × 10[−][2] 4.2 × 10[−][2] 9.9 × 10[−][4] 3.5 × 10[−][4] 2.8 × 10[−][4] 4.1 × 10[−][5]


Figure 12: Prediction mean over ten different trained models (with 95% confidence interval) of the
dissolved concentration in the diffusion-sorption equation at t = 10 000 for the in-dis-test dataset.

Table 7: Closed loop MSE on out-dis-test data from ten different training runs for each model for
the diffusion-sorption equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 4.7 × 10[−][2] 2.3 × 10[−][2] 2.9 × 10[−][3] -  1.3 × 10[−][2] 9.7 × 10[−][5]

2 2.1 × 10[−][1] 4.3 × 10[−][2] 1.6 × 10[−][3] -  1.1 × 10[−][2] 1.5 × 10[−][5]

3 9.2 × 10[−][2] 5.5 × 10[−][2] 3.4 × 10[−][3] -  1.4 × 10[−][2] 9.5 × 10[−][5]

4 2.0 × 10[−][1] 1.9 × 10[−][2] 2.0 × 10[−][3] -  1.1 × 10[−][2] 3.4 × 10[−][5]

5 4.4 × 10[−][2] 1.3 × 10[−][1] 4.6 × 10[−][4] -  2.1 × 10[−][2] 4.9 × 10[−][6]

6 3.9 × 10[−][4] 1.0 × 10[−][1] 8.7 × 10[−][3] -  4.4 × 10[−][2] 1.9 × 10[−][5]

7 9.3 × 10[−][4] 9.1 × 10[−][2] 1.6 × 10[−][2] -  1.5 × 10[−][2] 1.5 × 10[−][5]

8 5.0 × 10[−][3] 2.6 × 10[−][3] 9.5 × 10[−][4] -  1.2 × 10[−][2] 1.0 × 10[−][4]

9 4.6 × 10[−][1] 3.4 × 10[−][2] 7.1 × 10[−][3] -  1.3 × 10[−][2] 1.2 × 10[−][5]

10 4.1 × 10[−][2] 8.0 × 10[−][2] 2.6 × 10[−][3] -  1.3 × 10[−][2] 1.6 × 10[−][5]


Figure 13: Prediction mean over ten different trained models (with 95% confidence interval) of the
dissolved concentration in the diffusion-sorption equation at t = 10 000 for the out-dis-test dataset.


-----

defined as a feedforward network with the size of [2, 20, 20, 20, 20, 20, 20, 20, 20, 2]. PhyDNet
was defined with the PhyCell containing 32 input dimensions, 7 hidden dimensions, 1 hidden layer,
and the ConvLSTM containing 32 input dimensions, 32 hidden dimensions, 1 hidden layer. For
**FINN, the modules ϕ** and ϕ were used, with ϕ defined as a feedforward network with the size
_N_ _D_ _D_
of [1, 10, 20, 10, 1] that takes u as an input and outputs the retardation factor R(u). All models are
trained until convergence using the L-BFGS optimizer, except for PhyDNet, which is trained with
the Adam optimizer and a learning rate of 1 × 10[−][3] due to stability issues when training with the
L-BFGS optimizer.

**Additional results** Individual errors of the ten different training runs as well as visualizations for
the train (Table 5), in-dis-test (Table 6 and Figure 12), and out-dis-test (Table 7 and Figure 13)
datasets. Results for the total concentration ut were omitted due to high similarity to the concentration of the reported contamination solution u.

A.6 DIFFUSION-REACTION

The diffusion-reaction equation is applicable in physical and biological systems, for example in
pattern formation (Turing, 1952).

**Data** In the current paper, we consider the 2D diffusion-reaction for that class of problems:


_∂2u1_

_∂x[2][ +][ ∂]∂y[2][u][2][1]_




_∂2u2_

_∂x[2][ +][ ∂]∂y[2][u][2][2]_




_∂u1_

_∂t_ [=][ R][1][(][u][1][, u][2][) +][ D][1]


_∂u2_

_∂t_ [=][ R][2][(][u][1][, u][2][) +][ D][2]


(11)


(12)


Here, D1 = 10[−][3] and D2 = 5 × 10[−][3] are the diffusion coefficient for the activator and inhibitor, respectively. The system of equations is coupled through the reaction terms R1(u1, u2) and
_R2(u1, u2) which are both dependent on u1 and u2. In this work, the Fitzhugh-Nagumo (Klaasen_
& Troy, 1984) system was considered to define the reaction function:

_R1(u1, u2) = u1_ _u[3]1_ (13) _R2(u1, u2) = u1_ _u2,_ (14)
_−_ _[−]_ _[k][ −]_ _[u][2][,]_ _−_

with k = 5 × 10[−][3].

The simulation domain for the train data is defined with x = [−1, 1], y = [−1, 1], t = [0, 10] and
is discretized with Nx = 49 and Ny = 49 spatial locations, and Nt = 101 simulation steps. The
initial condition was defined as u1(x, 0) = u2(x, 0) = sin(π(x + 1)/2) sin(π(y + 1)/2), and the
corresponding boundary condition was defined as _[∂u]∂x[1]_ [(][−][1][, t][) = 0][,][ ∂u]∂x[1] [(1][, t][) = 0][,][ ∂u]∂x[2] [(][−][1][, t][) = 0][,]

and _[∂u]∂x[2]_ [(1][, t][) = 0][.]

Table 8: Closed loop MSE on train data from ten different training runs for each model for the
diffusion-reaction equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN


1 8.7 × 10[−][3] 1.7 × 10[−][3] 1.7 × 10[−][3] 2.7 × 10[−][4] 8.5 × 10[−][5] 1.0 × 10[−][4]

2 4.4 × 10[−][3] 1.6 × 10[−][3] 4.2 × 10[−][3] 5.6 × 10[−][4] 7.1 × 10[−][5] 1.6 × 10[−][4]

3 3.2 × 10[−][2] 1.4 × 10[−][3] 9.1 × 10[−][3] 2.1 × 10[−][4] 5.8 × 10[−][5] 1.4 × 10[−][4]

4 8.7 × 10[−][3] 7.8 × 10[−][4] 6.5 × 10[−][4] 2.8 × 10[−][4] 6.8 × 10[−][5] 1.2 × 10[−][4]

5 2.5 × 10[−][2] 1.5 × 10[−][3] 4.7 × 10[−][4] 2.8 × 10[−][4] 7.6 × 10[−][5] 1.0 × 10[−][4]

6 8.0 × 10[−][3] 4.8 × 10[−][4] 2.1 × 10[−][4] 2.5 × 10[−][4] 6.5 × 10[−][5] 1.7 × 10[−][4]

7 2.4 × 10[−][2] 3.7 × 10[−][3] 8.6 × 10[−][3] 6.1 × 10[−][5] 8.2 × 10[−][5] 1.6 × 10[−][4]

8 1.1 × 10[−][2] 1.0 × 10[−][3] 2.4 × 10[−][3] 1.0 × 10[−][4] 7.8 × 10[−][5] 1.0 × 10[−][4]

9 1.3 × 10[−][2] 2.2 × 10[−][3] 3.8 × 10[−][3] 3.2 × 10[−][5] 8.8 × 10[−][5] 1.5 × 10[−][4]

10 5.6 × 10[−][3] 7.3 × 10[−][2] 8.6 × 10[−][3] 6.2 × 10[−][4] 7.4 × 10[−][5] 1.1 × 10[−][4]


-----

Table 9: Closed loop MSE on in-dis-test data from ten different training runs for each model for the
diffusion-reaction equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 1.8 × 10[−][1] 8.0 × 10[−][2] 1.6 × 10[−][1] 7.7 × 10[−][2] 6.2 × 10[−][2] 2.9 × 10[−][3]

2 5.9 × 10[−][1] 5.7 × 10[−][2] 1.9 × 10[−][1] 2.3 × 10[−][1] 8.8 × 10[−][2] 1.6 × 10[−][3]

3 6.5 × 10[−][1] 9.2 × 10[−][2] 3.0 × 10[−][1] 3.5 × 10[−][2] 7.8 × 10[−][2] 1.7 × 10[−][3]

4 3.1 × 10[−][1] 6.7 × 10[−][2] 2.0 × 10[−][1] 2.6 × 10[−][2] 6.8 × 10[−][2] 2.9 × 10[−][3]

5 7.3 × 10[−][1] 4.9 × 10[−][2] 1.9 × 10[−][1] 1.3 × 10[−][1] 5.3 × 10[−][2] 1.7 × 10[−][3]

6 1.9 × 10[−][1] 5.4 × 10[−][2] 2.8 × 10[−][2] 8.9 × 10[−][2] 9.9 × 10[−][2] 2.0 × 10[−][3]

7 6.0 × 10[−][1] 1.0 × 10[−][1] 2.2 × 10[−][1] 5.4 × 10[−][2] 6.7 × 10[−][2] 1.5 × 10[−][3]

8 7.1 × 10[−][1] 8.6 × 10[−][2] 1.9 × 10[−][1] 1.3 × 10[−][2] 9.1 × 10[−][2] 2.3 × 10[−][3]

9 5.5 × 10[−][1] 2.3 × 10[−][1] 1.6 × 10[−][1] 2.3 × 10[−][2] 6.7 × 10[−][2] 1.8 × 10[−][3]

10 2.2 × 10[−][1] 1.2 × 10[−][1] 1.7 × 10[−][1] 2.6 × 10[−][2] 1.1 × 10[−][1] 2.1 × 10[−][3]


Figure 14: Prediction mean over ten different trained models (with 95% confidence interval) of the
activator in the diffusion-reaction equation at t = 50 for the in-dis-test dataset.

Table 10: Closed loop MSE on out-dis-test data from ten different training runs for each model for
the diffusion-reaction equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 2.3 × 10[−][2] 3.2 × 10[−][3] 9.4 × 10[−][3] -  3.9 × 10[−][2] 6.5 × 10[−][3]

2 1.7 × 10[−][1] 1.3 × 10[−][2] 1.3 × 10[−][2] -  5.4 × 10[−][2] 5.5 × 10[−][3]

3 2.5 × 10[−][1] 1.3 × 10[−][2] 2.1 × 10[−][2] -  2.5 × 10[−][2] 5.9 × 10[−][3]

4 7.3 × 10[−][2] 7.2 × 10[−][3] 7.9 × 10[−][3] -  4.2 × 10[−][2] 6.5 × 10[−][3]

5 2.6 × 10[−][1] 5.9 × 10[−][3] 2.6 × 10[−][3] -  2.8 × 10[−][2] 6.1 × 10[−][3]

6 6.4 × 10[−][2] 1.5 × 10[−][2] 2.1 × 10[−][3] -  2.9 × 10[−][2] 6.0 × 10[−][3]

7 2.2 × 10[−][1] 1.0 × 10[−][2] 1.9 × 10[−][2] -  2.8 × 10[−][2] 5.7 × 10[−][3]

8 2.0 × 10[−][1] 3.9 × 10[−][3] 1.2 × 10[−][2] -  2.1 × 10[−][2] 6.1 × 10[−][3]

9 2.3 × 10[−][1] 4.3 × 10[−][2] 1.3 × 10[−][2] -  2.1 × 10[−][2] 6.0 × 10[−][3]

10 3.6 × 10[−][2] 3.7 × 10[−][2] 3.3 × 10[−][2] -  5.9 × 10[−][2] 6.2 × 10[−][3]


Figure 15: Prediction mean over ten different trained models (with 95% confidence interval) of the
activator in the diffusion-reaction equation at t = 10 for the out-dis-test dataset.


-----

The in-dis-test data is simulated with with x = [−1, 1], y = [−1, 1] and with the time span of
_t = [10, 50] and Nt = 501. Initial condition is taken from the train data at t = 10 and boundary_
conditions are also identical to the train data.

The simulation domain for the out-dis-test data was identical with the train data, except for the
initial condition that was defined as u1(x, 0) = u2(x, 0) = sin(π(x + 1)/2) sin(π(y + 1)/2) 0.5,
_−_
i.e. subtracting 0.5 from the original initial condition.

**Model architectures** **TCN is designed to have two input- and output neurons, and one hidden**
layer of size 32. ConvLSTM has two input- and output neurons and one hidden layer of size 24.
The lateral and dynamic input- and output sizes of the DISTANA model are set to one and two,
respectively, while a hidden layer of size 32 is used. The pure ML models were trained on the first
70 time steps and validated on the remaining 30 time steps of the train data (applying early stopping).
Also, to prevent the pure ML models from diverging too much in closed loop, the boundary data are
fed into the models as done during teacher forcing. PINN is defined as a feedforward network with
the size of [3, 20, 20, 20, 20, 20, 20, 20, 20, 2]. PhyDNet is defined with the PhyCell containing
32 input dimensions, 49 hidden dimensions, 1 hidden layer, and the ConvLSTM containing 32 input
dimensions, 32 hidden dimensions, 1 hidden layer. For FINN, the modules ϕ, ϕ and Φψ are
_N_ _D_
used, with ϕ set as two learnable scalars that learn the diffusion coefficients D1 and D2, and Φψ
_D_
defined as a feedforward network with the size of [2, 20, 20, 20, 2] that takes u1 and u2 as inputs and
outputs the reaction functions R1(u1, u2) and R2(u1, u2). All models are trained until convergence
using the L-BFGS optimizer, except for PhyDNet, which is trained with the Adam optimizer and a
learning rate of 1 × 10[−][3] due to stability issues when training with the L-BFGS optimizer.

**Additional results** Individual errors of the ten different training runs as well as visualizations
for the train (Table 8), in-dis-test (Table 9 and Figure 14), and out-dis-test (Table 10 and Figure 15)
datasets. Results for the total inhibitor u2 were omitted due to high similarity to the reported activator
_u1._

A.7 ALLEN-CAHN

The Allen-Cahn equation is commonly employed in reaction-diffusion systems, in particular to
model phase separation in multi-component alloy systems (Raissi et al., 2019).

**Data** The 1D Allen-Cahn equation is written as

_∂u_

(15)

_∂t_ [=][ D∂]∂x[2][u][2][ +][ R][(][u][)][,]


where the main unknown is u, the reaction term is denoted as R(u) which is a function of u and the
diffusion coefficient is D = 10[−][4]. In the current work, the reaction term is defined as:

_R(u) = 5u −_ 5u[3], (16)

Figure 16: Plots of Allen-Cahn’s data and prediction of in-dis-test data using different models. The
plots in the first row show the solution over x and t (the red lines mark the transition from train to
_in-dis-test), the second row visualizes the best model’s solution distributed in x at t = 1._


-----

Figure 17: Plots of Allen-Cahn’s data and prediction of out-dis-test data using different models.
The plots in the first row show the solution over x and t, the second row visualizes the solution
distributed in x at t = 1. Right wagga

to reproduce the experiment conducted in the PINN paper (Raissi et al., 2019).

The simulation domain for the train data is defined with x = [−1, 1], t = [0, 0.5] and is discretized
with Nx = 49 spatial locations, and Nt = 201 simulation steps. The initial condition is defined as
_u(x, 0) = x[2]_ cos(πx), and periodic boundary condition is used, i.e. u(−1, t) = u(1, t).

**_In-dis-test data is simulated with x = [−1, 1] and a time span of t = [0.5, 1] and Nt = 401. Initial_**
condition is taken from the train data at t = 0.5 and boundary conditions are also similar to the train
data.

The simulation domain for the out-dis-test data is identical with the train data, except for the initial
condition that is defined as u(x, 0) = sin(πx/2).

**Model architectures** Both TCN and ConvLSTM are designed to have one input neuron, one
hidden layer of size 32, and one output neuron. The lateral and dynamic input- and output sizes
of the DISTANA model are set to one and a hidden layer of size 32 is used. The pure ML models
were trained on the first 150 time steps and validated on the remaining 50 time steps of the train
data (applying early stopping). Also, to prevent the pure ML models from diverging too much in
closed loop, the boundary data are fed into the models as done during teacher forcing. PINN was
defined as a feedforward network with the size of [2, 20, 20, 20, 20, 20, 20, 20, 20, 1] (8 hidden
layers, each contains 20 hidden neurons). PhyDNet was defined with the PhyCell containing 32
input dimensions, 7 hidden dimensions, 1 hidden layer, and the ConvLSTM containing 32 input
dimensions, 32 hidden dimensions, 1 hidden layer. For FINN, the modules ϕ, ϕ, and Φψ were
_N_ _D_
used, with ϕ defined as a learnable scalar that learns the diffusion coefficient D, and Φψ defined
_D_

Table 11: Closed loop MSE on the train data from ten different training runs for each model for the
Allen-Cahn equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 1.1 × 10[−][2] 7.1 × 10[−][2] 2.9 × 10[−][4] 1.6 × 10[−][5] 1.5 × 10[−][4] 1.8 × 10[−][5]

2 4.5 × 10[−][2] 3.0 × 10[−][2] 4.7 × 10[−][4] 4.2 × 10[−][5] 6.4 × 10[−][5] 8.2 × 10[−][6]

3 2.9 × 10[−][2] 1.2 × 10[−][2] 8.0 × 10[−][4] 4.2 × 10[−][5] 2.5 × 10[−][4] 5.4 × 10[−][6]

4 8.9 × 10[−][2] 6.5 × 10[−][3] 4.9 × 10[−][4] 1.6 × 10[−][6] 4.8 × 10[−][5] 4.7 × 10[−][7]

5 1.4 × 10[−][2] 8.0 × 10[−][2] 2.6 × 10[−][4] 1.3 × 10[−][5] 5.1 × 10[−][5] 2.3 × 10[−][5]

6 1.1 × 10[−][2] 4.6 × 10[−][3] 5.7 × 10[−][4] 6.4 × 10[−][6] 4.7 × 10[−][5] 2.5 × 10[−][6]

7 9.6 × 10[−][3] 7.1 × 10[−][3] 4.2 × 10[−][4] 1.4 × 10[−][5] 1.7 × 10[−][4] 1.1 × 10[−][5]

8 3.6 × 10[−][2] 5.3 × 10[−][4] 8.5 × 10[−][4] 2.1 × 10[−][5] 3.7 × 10[−][5] 2.9 × 10[−][7]

9 6.9 × 10[−][1] 4.0 × 10[−][1] 7.6 × 10[−][4] 2.0 × 10[−][5] 9.8 × 10[−][5] 3.0 × 10[−][7]

10 2.8 × 10[−][2] 1.8 × 10[−][4] 1.5 × 10[−][4] 3.8 × 10[−][5] 9.8 × 10[−][5] 2.9 × 10[−][7]


-----

Table 12: Closed loop MSE on in-dis-test data from ten different training runs for each model for
the Allen-Cahn equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 9.9 × 10[−][2] 1.2 × 10[0] 4.8 × 10[−][2] 7.0 × 10[−][3] 1.4 × 10[−][1] 5.6 × 10[−][5]

2 2.9 × 10[−][1] 8.6 × 10[−][1] 9.1 × 10[−][2] 1.8 × 10[−][1] 8.2 × 10[−][2] 2.4 × 10[−][5]

3 2.3 × 10[−][1] 3.6 × 10[−][1] 7.9 × 10[−][2] 4.3 × 10[−][1] 1.0 × 10[−][1] 1.6 × 10[−][5]

4 4.7 × 10[−][1] 3.4 × 10[−][1] 7.4 × 10[−][2] 4.6 × 10[−][3] 9.9 × 10[−][2] 2.1 × 10[−][6]

5 1.9 × 10[−][1] 5.5 × 10[−][1] 7.3 × 10[−][2] 2.2 × 10[−][3] 8.5 × 10[−][2] 9.3 × 10[−][5]

6 1.2 × 10[−][1] 3.2 × 10[−][1] 1.0 × 10[−][1] 2.0 × 10[−][3] 1.2 × 10[−][1] 1.1 × 10[−][5]

7 1.2 × 10[−][1] 5.4 × 10[−][1] 1.3 × 10[−][1] 5.0 × 10[−][4] 2.3 × 10[−][1] 3.4 × 10[−][5]

8 3.0 × 10[−][1] 5.6 × 10[−][1] 9.5 × 10[−][2] 1.3 × 10[−][1] 1.0 × 10[−][1] 1.3 × 10[−][6]

9 1.4 × 10[0] 1.3 × 10[0] 1.6 × 10[−][1] 1.2 × 10[−][2] 1.0 × 10[−][1] 1.3 × 10[−][6]

10 2.3 × 10[−][1] 5.5 × 10[−][2] 4.6 × 10[−][2] 3.3 × 10[−][2] 9.4 × 10[−][2] 1.3 × 10[−][6]


Figure 18: Prediction mean over ten different trained models (with 95% confidence interval) of the
Allen-Cahn equation at t = 1 for the in-dis-test dataset.

Table 13: Closed loop MSE on out-dis-test data from ten different training runs for each model for
the Allen-Cahn equation.

Run TCN ConvLSTM DISTANA PINN PhyDNet FINN

1 4.1 × 10[−][2] 5.3 × 10[−][2] 1.4 × 10[−][2] -  6.2 × 10[−][1] 7.4 × 10[−][5]

2 2.0 × 10[−][1] 3.8 × 10[−][1] 3.7 × 10[−][2] -  2.4 × 10[−][1] 2.9 × 10[−][5]

3 1.3 × 10[−][1] 1.5 × 10[−][1] 6.2 × 10[−][2] -  6.2 × 10[−][1] 1.8 × 10[−][5]

4 2.0 × 10[−][1] 2.9 × 10[−][1] 1.1 × 10[−][2] -  7.3 × 10[−][1] 3.0 × 10[−][6]

5 1.7 × 10[−][1] 8.3 × 10[−][2] 4.0 × 10[−][2] -  6.8 × 10[−][1] 1.2 × 10[−][4]

6 8.7 × 10[−][2] 2.9 × 10[−][1] 3.1 × 10[−][2] -  7.8 × 10[−][1] 1.5 × 10[−][5]

7 1.0 × 10[−][1] 1.5 × 10[−][1] 5.9 × 10[−][2] -  1.0 × 10[0] 4.3 × 10[−][5]

8 1.5 × 10[−][1] 5.3 × 10[−][1] 9.1 × 10[−][2] -  5.4 × 10[−][1] 1.9 × 10[−][6]

9 1.2 × 10[0] 1.4 × 10[0] 1.0 × 10[−][1] -  6.2 × 10[−][1] 1.8 × 10[−][6]

10 9.6 × 10[−][2] 7.5 × 10[−][2] 1.5 × 10[−][2] -  5.9 × 10[−][1] 1.9 × 10[−][6]


Figure 19: Prediction mean over ten different trained models (with 95% confidence interval) of the
Allen-Cahn equation at t = 1 for the out-dis-test data. The right-most plot shows the accurately
learned reaction function from the data samples.


-----

as a feedforward network with the size of [1, 10, 20, 10, 1] that takes u as an input and outputs
the reaction function R(u). All models are trained until convergence using the L-BFGS optimizer,
except for PhyDNet, which is trained with the Adam optimizer and a learning rate of 1 × 10[−][3] due
to stability issues when training with the L-BFGS optimizer.

**Additional results** Individual errors of the ten different training runs are reported as well as visualizations are generated for the train (Table 11), in-dis-test (Table 12, Figure 16 and Figure 18), and
_out-dis-test (Table 13, Figure 17 and Figure 19) datasets. Moreover, the right-most plot in Figure 19_
shows the dataset’s

A.8 BASELINE ASSESSMENT WITH POLYNOMIAL REGRESSION

In this section, we apply polynomial regression to show that the example problems chosen in this
work (i.e. Burger’s, diffusion-sorption, diffusion-reaction, and Allen-Cahn) are not easy to solve.
First, we use polynomial regression to fit the unknown variable u = f (x, t), similar to PINN. Figure 20 shows the prediction of u for each example, obtained using the fitted polynomial coefficients.
For the Burger’s equation, the train and in-dis-test predictions have MSE values of 5.0 × 10[−][2] and
4.1 × 10[−][2], respectively. For the diffusion-sorption equation, the train and in-dis-test predictions
have MSE values of 3.0 × 10[−][3] and 4.1 × 10[1], respectively. For the diffusion-reaction equation,
the train and in-dis-test predictions have MSE values of 1.7 × 10[−][2] and 3.9 × 10[4], respectively.



|Col1|0 0|
|---|---|

|Col1|0 0|
|---|---|

|Col1|0 0|
|---|---|


Burgers', polynomial (d=5)

Training Extrapolation

1 1

0.5 0.5

x x

0 0.0 0 0.0

0.5 0.5

1 1

0.0 0.5 1.0 0 1 2

t t

t = 1 t = 2

0.5

0.5

u 0.0 u 0.0

0.5

0.5

1 0 1 1 0 1

x x

Data Prediction

Diffusion reaction, polynomial (d=5)

Training Extrapolation

1 1

0.5

y 0 0.0 y 0 0.0

0.5

-1 -1

-1 0 1 -1 0 1

x x

t = 10, y = 0 t = 50, y = 0

0.5

700

u1 0 u1 350

-0.5 0

-1 0 1 -1 0 1

x x

Data Prediction


Diffusion sorption, polynomial (d=3)

Training Extrapolation

1.0 1.0

x 0.5 0.5 x 0.5 0.5

0.0 0.0 0.0 0.0

0 2500 0 10000

t t

t = 2500 t = 10000

1

20

u 0.5 u

10

0 0

0 0.5 1 0 0.5 1

x x

Data Prediction

Allen-Cahn, polynomial (d=10)

Training Extrapolation

1 0.5 1

0.5

x 0.0 x

0 0 0.0

0.5 0.5

1 1

0.00 0.25 0.50 0.0 0.5 1.0

t t

t = 0.5 t = 1

0.5 300

0.0

u 0.5 u 150

1.0

0

1 0 1 1 0 1

x x

Data Prediction


Figure 20: Prediction pairs for train and in-dis-test (left and right figure columns of the pairs) of the
Burger’s (left top, polynomial order 5), diffusion-sorption (top right, polynomial order 3), diffusionreaction (bottom left, polynomial order 5 since higher orders did not converge), and Allen-Cahn
(bottom right, polynomial order 10) equations using polynomial regression.


-----

Burgers'

1

0.5

x 0 0.0

0.5

1

0.0 0.5 1.0

t

t = 1

0.5

u 0.0

0.5

1 0 1

x

Data Prediction


Diffusion reaction

1

y 0 0.0

-1

-1 0 1

x

t = 10, y = 0

0.5

u1 0

-0.5

-1 0 1

x

Data Prediction


Allen-Cahn

1 0.5

0.0

x 0

0.5

1

0.00 0.25 0.50

t

t = 0.5

0.5

0.0

u

0.5

1.0

1 0 1

x

Data Prediction


(a)


(b)


(c)


Figure 21: Plots of the train prediction in the Burger’s (left), diffusion-reaction (center), and AllenCahn equations (right) using FINN with polynomial fitting. Due to instability issues, the diffusionsorption equation could not be solved with the polynomial FINN. The plots in the first row show the
solution over x and t, and the plots in the second row show the solution distributed in x.

For the Allen-Cahn equation, the train and in-dis-test predictions have MSE values of 8.8 × 10[−][3]
and 2.2 × 10[3], respectively. The simple polynomial fitting fails to obtain accurate predictions of the
solution for all example problems. The results also show that the polynomials overfit the data, evidenced by the significant deterioration of performance during extrapolation (prediction of in-dis-test
data). The diffusion-reaction and the Allen-Cahn equations are particularly the most difficult to fit,
because they require higher order polynomials to obtain reasonable accuracy. With the high order,
they still fail to even fit the train data well.

Next, we also consider using polynomial fitting in lieu of ANNs (namely the modules ϕ, ϕ, and
_A_ _D_
Φψ) in FINN. With this method, the model successfully obtain accurate prediction of the Burger’s
equation (Figure 21a). The MSE values are 1.9 × 10[−][4], 1.0 × 10[−][3], and 1.3 × 10[−][4] for train, in_dis-test, and out-dis-test data, respectively. However, the model fails to complete the training for the_
diffusion-sorption equation due to major instabilities (the polynomials can produce negative output
and therefore, negative diffusion coefficient, leading to numerical instability). Moreover, the model
also fails to learn the diffusion-reaction (Figure 21b) and the Allen-Cahn (Figure 21c) equations with
good accuracy. For the diffusion-reaction equation, the MSE values are 2.5 _×_ 10[−][2], 1.7 _×_ 10[−][1], and
4.8 _×_ 10[−][2] for train, in-dis-test, and out-dis-test data, respectively. For the Allen-Cahn equation, the
MSE values are 5.6 × 10[−][2], 2.5 × 10[−][1], and 4.3 × 10[−][1] for train, in-dis-test, and out-dis-test data,
respectively. Even though the unknown equations do not seem too complicated, it is still difficult to
solve them together with the PDE as a whole. The results show that for these particular problems,
ANNs serve better because they allow better control during training, for example with constraints,
and they produce more regularized outputs than high order polynomials. However, ANNs are not
unique in their selection, but they are more convenient for our implementation.

A.9 ROBUSTNESS TEST OF FINN

In this section, we test the robustness of FINN when trained using noisy data. All the synthetic data
is generated with the same parameters, only added with noise with the distribution N (0.0, 0.05). For
the Burger’s equation (Figure 22a and Figure 23a), the average MSE values are 2.5 × 10[−][3] _± 4.1 ×_
10[−][6], 2.4 × 10[−][3] _± 6.6 × 10[−][6], and 2.5 × 10[−][3]_ _± 4.0 × 10[−][6]_ for the train, in-dis-test, and out-dis_test prediction, respectively. For the diffusion-sorption equation (Figure 22b and Figure 23b), the_
average MSE values are 2.5×10[−][3]±4.5×10[−][6], 2.5×10[−][3]±3.7×10[−][6], and 2.5×10[−][3]±3.7×10[−][6]
for the train, in-dis-test, and out-dis-test prediction, respectively. For the diffusion-reaction equation


-----

(a) Burger’s (b) Diffusion-sorption

(c) Diffusion-reaction (d) Allen-Cahn

Figure 22: Paired plots of the in-dis-test and out-dis-test prediction (left and right of the pairs,
respectively) in the Burger’s (top left), diffusion-sorption (top right), diffusion-reaction (bottom left),
and Allen-Cahn (bottom right) equations obtained by training FINN with noisy data. The plots in
the first row of the pairs show the solution over x and t (red line marks the transition from train to
_in-dis-test), and the plots in the second row of the pairs show the best model’s solution distributed_
in x.

(Figure 22c and Figure 23c), the average MSE values are 3.2 × 10[−][3] _± 5.3 × 10[−][4], 1.6 × 10[−][2]_ _±_
8.8 × 10[−][3], and 8.3 × 10[−][3] _± 4.9 × 10[−][4]_ for the train, in-dis-test, and out-dis-test prediction,
respectively. For the Allen-Cahn equation (Figure 22d and Figure 23d), the average MSE values are
2.5 × 10[−][3] _± 1.1 × 10[−][6], 2.5 × 10[−][3]_ _± 9.1 × 10[−][6], and 2.5 × 10[−][3]_ _± 6.3 × 10[−][6]_ for the train,
_in-dis-test, and out-dis-test prediction, respectively. These results show that even though FINN is_
trained with noisy data, it is still able to capture the essence of the equation and generalize well to
different initial and boundary conditions. Additionally, the prediction is consistent, shown by the
low values of the MSE standard deviation, as well as the very narrow confidence interval in the plots.

A.10 TRAINING PINN WITH FINER SPATIAL RESOLUTION

In order to determine whether the reduced accuracy of PINN in our experiments was caused by a
coarse spatial resolution (we only used 49, 26 and 49 × 49 spatial positions at Burger, diffusionsorption, and diffusion-reaction, respectively), another experiment was conducted per target equation where the spatial resolution was increased to Nx = 999, Nx = 251, and Nx = 99, Ny = 99,
respectively. As reported in Table 14, Figure 24, and Figure 25a, the performance increased slightly
but by far did not reach FINN’s accuracy. Identical results were achieved in the Burger’s and
diffusion-sorption equations but are omitted due to high conceptual similarity.


-----

(a) Burger’s (b) Diffusion-sorption

(c) Diffusion-reaction (d) Allen-Cahn

Figure 23: Prediction mean over ten different trained FINN (with 95% confidence interval) of the
Burger’s (top left), diffusion-sorption (top right), diffusion-reaction (bottom left), and Allen-Cahn
(bottom right) equations obtained by training FINN with noisy data for the in-dis-test and out-dis_test (top and bottom, accordingly) prediction._


-----

(a) In-dis-test prediction (left) and out-dis-test (right) (b) In-dis-test prediction (left) and out-dis-test (right)

Figure 24: Both left: plots of the diffusion-reaction equation’s activator u using PINN trained with
finer resolution dataset. Both right: plots of the diffusion-reaction equation’s activator u using
PhyDNet with the original network size.


Table 14: MSE of PINN trained using data with
finer resolution from ten different training runs
for the diffusion-reaction equation.

Run _Train_ _In-dis-test_ _Out-dis-test_

1 1.4×10[−][4] 7.8×10[−][2] - 
2 1.8×10[−][4] 5.6×10[−][2] - 
3 8.6×10[−][5] 1.5×10[−][2] - 
4 1.0×10[−][3] 1.5×10[−][1] - 
5 5.2×10[−][5] 2.7×10[−][2] - 
6 2.9×10[−][4] 8.1×10[−][1] - 
7 1.3×10[−][4] 3.4×10[−][2] - 
8 8.9×10[−][5] 2.9×10[−][2] - 
9 3.9×10[−][5] 1.3×10[−][2] - 
10 3.3×10[−][5] 1.7×10[−][2] - 


Table 15: MSE of PhyDNet using the original
network size from ten different training runs for
the diffusion-reaction equation.

Run _Train_ _In-dis-test_ _Out-dis-test_

1 9.3×10[−][5] 1.8×10[−][1] 7.6 × 10[−][2]

2 2.9×10[−][5] 6.5×10[−][2] 3.5 × 10[−][2]

3 2.2×10[−][5] 6.7×10[−][2] 4.0 × 10[−][2]

4 4.6×10[−][5] 7.2×10[−][2] 3.9 × 10[−][2]

5 3.5×10[−][5] 6.5×10[−][2] 8.8 × 10[−][2]

6 5.7×10[−][5] 1.2×10[−][1] 1.0 × 10[−][1]

7 3.5×10[−][5] 5.7×10[−][2] 9.0 × 10[−][2]

8 3.5×10[−][5] 5.9×10[−][2] 2.6 × 10[−][2]

9 2.3×10[−][5] 6.3×10[−][2] 4.1 × 10[−][2]

10 3.5×10[−][5] 6.1×10[−][2] 8.1 × 10[−][2]


A.11 PHYDNET WITH ORIGINAL AMOUNT OF PARAMETERS

To verify whether our reduction of parameters and the removal of the encoder and decoder layers
caused PhyDNet to perform worse, we repeated the experiments for the three equations of interest
using the original PhyDNet architecture as proposed in Guen & Thome (2020). However, our results
indicate no significant changes in performance, as reported in Table 15, Figure 24b, and Figure 25b.
Again, results for the inhibitor u2 as well as for the Burger’s and diffusion-reaction equations were
omitted due to high conceptual similarity.

A.12 SOIL PARAMETERS AND SIMULATION DOMAINS FOR THE EXPERIMENTAL DATASET

Identical to Praditia et al. (2021), the soil parameters and simulation (and experimental) domain
used in the real-world diffusion-sorption experiment are summarized in Table 16 for core samples
#1, #2, and #2B.

For all experiments, the core samples are subjected to a constant contaminant concentration at the
top us, which can be treated as a Dirichlet boundary condition numerically. Notice that, for core
sample #2, we set us to be slightly higher to compensate for the fact that there might be fractures at
the top of core sample #2, so that the contaminant can break through the core sample faster.


-----

(a) PINN fine resolution. (b) PhyDNet original architecture

Figure 25: Prediction mean (with 95% confidence interval) of the activator (top) and inhibitor (bottom) in the diffusion-reaction equation at t = 50 compared with the in-dis-test data.

For core samples #1 and #2, Q is the flow rate of clean water at the bottom reservoir that determines
the Cauchy boundary condition at the bottom of the core samples. For core sample #2B, note that
the sample length is significantly longer than the other samples, and by the end of the experiment,
no contaminant has broken through the core sample. Therefore, we assume the bottom boundary
condition to be a no-flow Neumann boundary condition; see (Praditia et al., 2021) for details.

Table 16: Soil and experimental parameters of core samples #1, #2, and #2B. D is the diffusion
coefficient, φ is the porosity, ρs is the bulk density, L and r are the length and radius of the sample,
_tend is the simulation time, Q is the flow rate in the bottom reservoir and us is the total concentration_
of trichloroethylene in the sample.

Soil parameters

Parameter Unit Core #1 Core #2 Core #2B

_D_ m[2]/day 2.00 × 10[−][5] 2.00 × 10[−][5] 2.78 × 10[−][5]

_φ_ -  0.288 0.288 0.288
_ρs_ kg/m[3] 1957 1957 1957

Simulation domain

Parameter Unit Core #1 Core #2 Core #2B

_L_ m 0.0254 0.02604 0.105
_r_ m 0.02375 0.02375 N/A
_tend_ days 38.81 39.82 48.88
_Q_ m[3]/day 1.01 × 10[−][4] 1.04 × 10[−][4] N/A
_us_ kg/m[3] 1.4 1.6 1.4


-----

Neural ODE runtime comparison

1.2

CPU
GPU

1.0

0.8

0.6

Runtime [s]

0.4

0.2

0 10000 20000 30000 40000 50000 60000

Number of batches


Figure 26: Runtime comparison of Neural ODE with a single hidden layer consisting of 50 hidden
nodes run on CPU and GPU for 1 000 time steps. The benefit of using GPU starts to be seen when
larger batch sizes (> 20 000) are used.

A.13 ON TIME STEP ADAPTIVITY AND RUNTIME

We compare the runtime for each model, run on a CPU with i9-9900K core, a clock speed of 3.60
GHz, and 32 GB RAM. Additionally, we also perform the comparison of GPU runtime on a GTX
1060 (i.e. with 6GB VRAM). The results are summarized in Table 17. Note that the purpose of this
comparison is not an optimized benchmark, but only to show that the runtime of FINN is comparable
with the other model, especially when run in CPU. When run on GPU, however, FINN runs slightly
slower. This is caused by the fact that FINN’s implementation is not yet optimized for GPU. More
importantly, the Neural ODE package we use benefits only from a larger batch size. As shown in
Figure 26, GPU is faster for batch size larger than 20 000, whereas the maximum size that we use in
the example is 2 401. With smaller batch size, CPU usage is faster.

In general, only the PINN model has a benefit when computed on the GPU, since the function is
only called once on all batches. This is different for all other models that have to unroll a prediction
of the sequence into the future recurrently (except for TCN which is a convolution approach that
is faster on GPU). Accordingly, The overhead of copying the tensors to GPU outweighs the GPU’s
parallelism benefit, compared to directly processing the sequence on the CPU iteratively. On the
two-dimensional diffusion-reaction benchmark, the GPU’s speed-up comes into play, since in here,
the simulation domain is discretized into 49 × 49 = 2401 volumes (compared to 49 for the Burger’s
and Allen-Cahn, and 26 for the diffusion-sorption equations). Note that we have observed significantly varying runtimes on different GPUs (i.e. up to three seconds for FINN on Burger’s on an
RTX 3090), which might be caused by lack of support of certain packages for a particular hardware,
but further investigation is required.

Additionally, we want to emphasize that the higher runtime of FINN on GPU is not caused by the
time step adaptivity. In fact, employing the adaptive time stepping strategy is cheaper than choosing
all time steps to be small enough (to guarantee numerical stability). As we learn a PDE, we have
no exact knowledge to derive a dedicated time integration scheme, but the adaptive Runge-Kutta
method is one of the best generic choices. As our PDE and its characteristics change during training, time step adaptivity is a real asset, because for example the Courant–Friedrichs–Lewy (CFL)


-----

Table 17: Comparison of runtime (in seconds) of single forward passes between different deep
learning (above dashed line) and physics-aware neural network (below dashed line) methods on the
different equations.

Eqn. Model CPU GPU (GTX 1060)


TCN 0.423 0.130
ConvLSTM 0.052 0.079
DISTANA 0.059 0.098

PINN 0.036 0.007
PhyDNet 0.107 0.192
FINN 0.066 0.161

TCN 1.228 0.393
ConvLSTM 0.119 0.194
DISTANA 0.145 0.223

PINN 0.073 0.010
PhyDNet 0.263 0.475
FINN 0.676 1.638

TCN 14.15 0.800
ConvLSTM 0.230 0.052
DISTANA 0.190 0.051

PINN 1.647 0.787
PhyDNet 0.159 0.113
FINN 0.342 0.330

TCN 0.442 0.128
ConvLSTM 0.050 0.082
DISTANA 0.060 0.097

PINN 0.035 0.007
PhyDNet 0.108 0.191
FINN 0.028 0.071


Burger

Diffusion-sorption

Diffusion-reaction

Allen-Cahn


condition (Courant et al., 1967) would consistently change throughout the training. Therefore, the
time step adaptivity is not a bottleneck, but rather a solution for a more efficient computation.

Furthermore, in relation to FINN’s limitation (see subsection A.2) topics like numerical stabilization
schemes for larger time steps, adaptive spatial grid, optimized implementation in an High Performance Computing (HPC) setting, parallelization, etc. are still very interesting for future works.


-----

