summaries: This paper presents a compositional physics-aware neural network (FINN) for learning spatiotemporal advection-diffusion processes. It claims that the FINN outperforms pure machine learning and other state-of-the-art physics-aware models in all cases—often even by multiple orders of magnitude. However, the design of the network depends too much on the form of the equation, which leads to a very narrow application of the method. In this work, the authors propose the finite volume neural network to solve advection-diffusion partial differential equations. The authors define a flux kernel as the sum of sub-kernel f_i on each element. And the authors define specific modules \phi_D, \phi_A, \phi_N according to the form of the advection-diffusion equation. The model requires very few parameters to achieve the state of art results. The paper introduces Finite Volume neural network for modelling fluid dynamics inspired by the finite volume method. The paper models the velocity and the spacial derivatives as the neural networks. The work demonstrates more precise fluid simulation than the related physics-inspired models. The authors propose to model advection-diffusion partial differential equations as a composition of multiple neural networks. According to the authors, this leads to better generalization over different initial and boundary conditions for the physics system, and also the ability to learn different factors of the process as modeled by the terms in the PDE. Extensive experimental results are presented to support the author’s claims that their framework performs better than state-of-the-art.
abstractText: 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 oneand 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.
references.authors: Y. Bar-Sinai S. Hoyer J. Hickey M.P. Brenner
references.title: Learning data-driven discretizations for partial differential equations
references.venue: Proceedings of the National Academy of Sciences,
references.authors: C. Basdevant M. Deville P. Haldenwang J.M. Lacroix J. Ouazzani R. Peyret P. Orlandi A.T. Patera
references.title: Spectral and finite difference solutions of the burgers equation
references.venue: Computers & Fluids,
references.authors: P.W. Battaglia J.B. Hamrick V. Bapst A. Sanchez-Gonzalez V. Zambaldi M. Malinowski A. Tacchetti D. Raposo A. Santoro R. Faulkner
references.title: Relational inductive biases, deep learning, and graph networks
references.venue: arXiv preprint arXiv:1806.01261,
references.authors: J.C. Butcher
references.title: Numerical Methods for Ordinary Differential Equations
references.authors: R.T.Q. Chen Y. Rubanova J. Bettencourt D.K. Duvenaud
references.title: Neural ordinary differential equations
references.venue: arXiv preprint arXiv:1806.07366,
references.authors: R. Courant K. Friedrichs H. Lewy
references.title: On the partial difference equations of mathematical physics
references.venue: IBM Journal of Research and Development,
references.authors: G. Cybenko
references.title: Approximation by superpositions of a sigmoidal function
references.venue: Mathematics of control, signals and systems,
references.authors: B. Fornberg
references.title: Generation of finite difference formulas on arbitrarily spaced grids
references.venue: Mathematics of computation,
references.authors: V.L. Guen N. Thome
references.title: Disentangling physical dynamics from unknown factors for unsupervised video prediction
references.venue: In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition,
references.authors: E. Isaacson H.B. Keller
references.title: Analysis of Numerical Methods
references.venue: Dover Books on Mathematics. Dover Publications,
references.authors: N. Kalchbrenner L. Espeholt K. Simonyan A. Oord A. Graves K. Kavukcuoglu
references.title: Neural machine translation in linear time
references.venue: arXiv preprint arXiv:1610.10099,
references.authors: M. Karlbauer S. Otte H.P.A. Lensch T. Scholten V. Wulfmeyer M.V. Butz
references.title: A distributed neural network architecture for robust non-linear spatio-temporal prediction
references.venue: arXiv preprint arXiv:1912.11141,
references.authors: G.A. Klaasen W.C. Troy
references.title: Stationary wave solutions of a system of reaction-diffusion equations derived from the fitzhugh–nagumo equations
references.venue: SIAM Journal on Applied Mathematics,
references.authors: D. Kochkov J.A. Smith A. Alieva Q. Wang M.P. Brenner S. Hoyer
references.title: Machine learning– accelerated computational fluid dynamics
references.venue: Proceedings of the National Academy of Sciences,
references.authors: B.M. Lake
references.title: Compositional generalization through meta sequence-to-sequence learning
references.venue: In Advances in Neural Information Processing Systems
references.authors: B.M. Lake T.D. Ullman J.B. Tenenbaum S.J. Gershman
references.title: Building machines that learn and think like people
references.venue: Behavioral and Brain Sciences,
references.authors: Z. Li N. Kovachki K. Azizzadenesheli B. Liu K. Bhattacharya A. Stuart A. Anandkumar
references.title: Fourier neural operator for parametric partial differential equations
references.venue: arXiv preprint arXiv:2010.08895,
references.authors: Z. Li N. Kovachki K. Azizzadenesheli B. Liu K. Bhattacharya A. Stuart A. Anandkumar
references.title: Neural operator: Graph kernel network for partial differential equations
references.venue: arXiv preprint arXiv:2003.03485,
references.authors: Z. Long Y. Lu X. Ma B. Dong
references.title: Pde-net: Learning pdes from data
references.venue: In International Conference on Machine Learning,
references.authors: F. Moukalled L. Mangani M. Darwish
references.title: The Finite Volume Method in Computational Fluid Dynamics. Springer, 1 edition, 2016
references.authors: W. Nowak A. Guthke
references.title: Entropy-based experimental design for optimal model discrimination in the geosciences
references.authors: T. Praditia T. Walser S. Oladyshkin W. Nowak
references.title: Improving thermochemical energy storage dynamics forecast with physics-inspired neural network architecture. Energies
references.authors: T. Praditia M. Karlbauer S. Otte S. Oladyshkin M.V. Butz W. Nowak
references.title: Finite volume neural network: Modeling subsurface contaminant transport
references.venue: arXiv preprint arXiv:2104.06010,
references.authors: M. Raissi P. Perdikaris G.E. Karniadakis
references.title: Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations
references.venue: Journal of Computational Physics,
references.authors: K.F. Riley M.P. Hobson S.J. Bence
references.title: Mathematical Methods for Physics and Engineering
references.authors: S. Seo C. Meng Y. Liu
references.title: Physics-aware difference graph networks for sparsely-observed dynamics
references.venue: In International Conference on Learning Representations,
references.authors: X. Shi Z. Chen H. Wang D.Y. Yeung W.K. Wong W.C. Woo
references.title: Convolutional lstm network: A machine learning approach for precipitation nowcasting
references.venue: arXiv preprint arXiv:1506.04214,
references.authors: V. Sitzmann J. Martel A Bergman D. Lindell G. Wetzstein
references.title: Implicit neural representations with periodic activation functions
references.venue: Advances in Neural Information Processing Systems,
references.authors: P.K. Smolarkiewicz
references.title: A simple positive definite advection scheme with small implicit diffusion
references.venue: Monthly Weather Review,
references.authors: A. Turing
references.title: The chemical basis of morphogenesis
references.venue: Philosophical Transactions of the Royal Society B,
references.authors: H.K. Versteeg W. Malalasekera
references.title: An Introduction to Computational Fluid Dynamics: The Finite Volume Method
references.venue: New York,
references.authors: Zhuang
references.title: 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)
references.venue: Kochkov et al.,
sections.heading: 1 INTRODUCTION
sections.text: 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.
sections.heading: 2 RELATED WORK
sections.text: 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.
sections.heading: 3 FINITE VOLUME NEURAL NETWORK (FINN)
sections.text: Problem formulation Here, we focus on modeling spatiotemporal physical processes. Specifically, we consider systems governed by advection-diffusion type equations (Smolarkiewicz, 1983):
  ∂u ∂t = D(u)
  ∂2u
  ∂x2 − v(u)∂u ∂x + q(u), (1)
  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 diffusive flux D(u)∂
  2u ∂x2 . 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:
  Fi = Nsi∑
  j=1
  fj ≈ ∮
  ω⊆Ω
  ( D(u) ∂2u
  ∂x2 − v(u)∂u ∂x
  ) · n̂ dΓ, (2)
  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 ϕN is a linear layer with the purpose to learn the numerical 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 ϕN should amount to [−1, 1] with respect to [ui, ui−1] and [ui, ui+1] (i.e. simple difference between neighboring control volumes) in ideal one-dimensional problems.
  The other modules ϕA and ϕD are responsible for advective and diffusive fluxes, respectively. Both modules receive only ui as input. Moreover, the moduleR 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 ϕA is performed only on one control volume surface (either left or right, and not both at the same time), whereas calculation of the diffusive fluxes using ϕD is performed on both surfaces. This means, for the advective 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 ϕA > 0 or [ui, ui+1] when ϕA < 0 (i.e. only 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 ϕD is performed on both surfaces (see subsection A.3 for more details). Accordingly, module R generates the inductive bias to make ϕA only approximate the advective, and ϕD the diffusive flux.
  For the advective flux calculation, if the velocity v is a function of u, the module ϕA approximates this dependence using a feedforward neural network, i.e. ϕA(u) ≈ v(u). Otherwise, ϕA is a scalar value ϕA ≡ v, which can be also set as a learnable parameter. The output of ϕA is activated by the following case-sensitive ReLUR:
  R(ϕA) = {
  ReLU(ϕA), on fi−, −ReLU(−ϕA), on fi+,
  (3)
  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 ϕD takes u as an input to learn the function ϕD(u) ≈ D(u) using a feedforward neural network. Otherwise, ϕD is a scalar value ϕD ≡ D, which can also be set as a learnable parameter. The output of ϕD is not activated by the module R, but is directly multiplied with the output of ϕN to obtain 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:
  Fi = { ϕN (ui, ui−1)[ϕD − ϕA]fi− + ϕN (ui, ui+1)[ϕD]fi+ , if ϕA > 0, ϕN (ui, ui−1)[ϕD]fi− + ϕN (ui, ui+1)[ϕD − ϕA]fi+ , otherwise.
  (4)
  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) ≈ ∂ui ∂t . (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.
sections.heading: 4 EXPERIMENTS, RESULTS & DISCUSSION
sections.heading: 4.1 SYNTHETIC DATASET
sections.text: 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, outdis-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
  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 outdis-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
  perform well. However, FINN still outperforms the other models with a significantly lower prediction error. The advective velocity learned by FINN’s module ϕA is shown in Figure 8 (left) and 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 ϕD module is shown in Figure 8 (second from left). The plot shows that 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 (outdis-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
  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-
  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.
sections.heading: 4.2 EXPERIMENTAL DATASET
sections.text: 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
  (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 ϕD to learn the retardation 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.
sections.heading: 5 CONCLUSION
sections.text: 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 ϕN to approximate variable 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.
sections.heading: AUTHOR CONTRIBUTIONS
sections.text: Hidden for anonymization reasons.
sections.heading: ACKNOWLEDGMENTS
sections.text: 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
sections.heading: REPRODUCIBILITY STATEMENT
sections.text: 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 anonymous version of the repository’s README.md file along with the according data and model scripts.]
sections.heading: A APPENDIX
sections.text: 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.
sections.heading: A.1 DISTINCTION TO RELATED WORK
sections.text: 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 = F
  ( t, x, y, . . . , u, ∂u
  ∂x , ∂u y , ∂u ∂x∂y , ∂2u ∂x2 , . . .
  ) . (6)
  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 ϕN and combining this stencil 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.
sections.heading: A.2 LIMITATIONS OF FINN
sections.text: 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.
sections.heading: A.3 LEARNING THE NUMERICAL STENCIL
sections.text: Semantically, the ϕN module learns the numerical stencil, that is the geometrical arrangement of a group of neighboring cells to approximate the derivative numerically, effectively learning the first spatial derivative ∂u∂x from both [ui−1, ui] and [ui, ui+1], which are the inputs to the ϕNi− and ϕNi+ module, respectively.
  The lateral information flowing from ui−1 and ui+1 toward ui is controlled by the ϕA (advective flux, i.e. bulk motion of many particles/atoms that can either move to the left or to the right) and the ϕD (diffusive flux, i.e. drive of particles/atoms to equilibrium from regions of high to low 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 combination of the ϕN and ϕD modules, calculating the diffusive fluxes δ− = ϕN (ui−1, ui)ϕD(ui) and δ+ = ϕN (ui, ui+1)ϕD(ui) inside of the respective left and right flux kernel. The combination of 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
  ∂2u ∂x2 ≈ (∂u/∂x)|i− − (∂u/∂x)|i+
  ∆x with (∂u/∂x)|i− ≈ (ui−1 − ui)/∆x and (∂u/∂x)|i+ ≈ (ui − ui+1)/∆x. Then
  ∂2u ∂x2 ≈ (ui−1 − ui)− (ui − ui+1) ∆x2 = ui−1 − 2ui + ui+1 ∆x2
  hence the [+1,−2,+1] as coefficients.
sections.heading: A.4 BURGER
sections.text: 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
  ∂t = −v(u)∂u ∂x +D
  ∂2u ∂x2 , (7)
  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.
  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 ϕN , ϕD, R and ϕA were used, with ϕA defined as a feedforward network with the size of [1, 10, 20, 10, 1] that takes u as an input and outputs the advective velocity v(u), and ϕD as 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.
sections.heading: A.5 DIFFUSION-SORPTION
sections.text: 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:
  ∂u ∂t = D R(u)
  ∂2u ∂x2 , (8) ∂ut
  ∂t = Dφ
  ∂2u ∂x2 , (9)
  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− φ φ ρsknfu nf−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
  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 ϕN and ϕD were used, with ϕD defined as a feedforward network with the size 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.
sections.heading: A.6 DIFFUSION-REACTION
sections.text: 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:
  ∂u1 ∂t = R1(u1, u2) +D1 ( ∂2u1 ∂x2 + ∂2u1 ∂y2 ) ,
  (11)
  ∂u2 ∂t = R2(u1, u2) +D2 ( ∂2u2 ∂x2 + ∂2u2 ∂y2 ) .
  (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 − u31 − k − u2, (13) R2(u1, u2) = u1 − u2, (14)
  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 ∂u1∂x (−1, t) = 0, ∂u1∂x (1, t) = 0, ∂u2∂x (−1, t) = 0, and ∂u2∂x (1, t) = 0.
  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 ϕN , ϕD and Φψ are used, with ϕD set as two learnable scalars that learn the diffusion coefficients D1 and D2, and Φψ 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.
sections.heading: A.7 ALLEN-CAHN
sections.text: 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 ∂t = D
  ∂2u ∂x2 +R(u), (15)
  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− 5u3, (16)
  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) = x2 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 ϕN , ϕD, and Φψ were used, with ϕD defined as a learnable scalar that learns the diffusion coefficient D, and Φψ defined
  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
sections.heading: A.8 BASELINE ASSESSMENT WITH POLYNOMIAL REGRESSION
sections.text: 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 × 101, respectively. For the diffusion-reaction equation, the train and in-dis-test predictions have MSE values of 1.7 × 10−2 and 3.9 × 104, respectively.
  For the Allen-Cahn equation, the train and in-dis-test predictions have MSE values of 8.8 × 10−3 and 2.2× 103, 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 ϕA, ϕD, and Φψ) 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, indis-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.
sections.heading: A.9 ROBUSTNESS TEST OF FINN
sections.text: 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 distributionN (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-distest 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
  (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.
sections.heading: A.10 TRAINING PINN WITH FINER SPATIAL RESOLUTION
sections.text: 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.
sections.heading: A.11 PHYDNET WITH ORIGINAL AMOUNT OF PARAMETERS
sections.text: 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.
sections.heading: A.12 SOIL PARAMETERS AND SIMULATION DOMAINS FOR THE EXPERIMENTAL DATASET
sections.text: 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.
  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.
sections.heading: A.13 ON TIME STEP ADAPTIVITY AND RUNTIME
sections.text: 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)
  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.
title: COMPOSING PARTIAL DIFFERENTIAL EQUATIONS