2023

[1,2]Shubhaditya Burela

[1]Institute of Mathematics, Technische Universität Berlin, Straße des 17.
Juni 136, Berlin, 10623, Berlin, Germany 2]Institute of Fluid Mechanics and
Technical Acoustics, Technische Universität Berlin, Müller-Breslau-Straße 15,
Berlin, 10623, Berlin, Germany 3]Institut de Mathématiques de Marseille (I2M),
Aix-Marseille Université, 39 Rue Joliot-Curie, Marseille, 13453, France

# Parametric model order reduction for a wildland fire model via the shifted
POD based deep learning method

<EMAIL_ADDRESS>Philipp Krah<EMAIL_ADDRESS>Julius Reiss
<EMAIL_ADDRESS>*  [  [

###### Abstract

Parametric model order reduction techniques often struggle to accurately
represent transport-dominated phenomena due to a slowly decaying Kolmogorov
n-width. To address this challenge, we propose a non-intrusive, data-driven
methodology that combines the shifted proper orthogonal decomposition (POD)
with deep learning. Specifically, the shifted POD technique is utilized to
derive a high-fidelity, low-dimensional model of the flow, which is
subsequently utilized as input to a deep learning framework to forecast the
flow dynamics under various temporal and parameter conditions. The efficacy of
the proposed approach is demonstrated through the analysis of one- and two-
dimensional wildland fire models with varying reaction rates, and its
performance is evaluated using multiple error measures. The results indicate
that the proposed approach yields highly accurate results within the percent
range, while also enabling rapid prediction of system states within seconds.

###### keywords:

Model order reduction, shifted proper orthogonal decomposition, data-driven
models, deep learning, artificial neural networks, wildland fires

## Acknowledgement

We gratefully acknowledge the support of the Deutsche Forschungsgemeinschaft
(DFG) as part of GRK2433 DAEDALUS. The authors were granted access to the HPC
resources of IDRIS under allocation No. AD012A01664R1 attributed by Grand
Équipement National de Calcul Intensif (GENCI). Centre de Calcul Intensif
d’Aix-Marseille is acknowledged for granting access to its high-performance
computing resources. We also thank Tobias Breiten for his valuable comments
and feedback.

##  1 Introduction

Forecasting the spread of forest fires has become an important aspect of civil
protection european_commission_joint_research_centre_european_2020 . The
intensification of extreme droughts and heat waves due to climate change has
increased the frequency and severity of large forest fires worldwide.
Consequently, predicting the risk and identifying the causes behind these
events has become crucial in comprehending the connection between climate and
land surface and helping manage forest fires kondylatos_wildfire_2022 .
Numerous studies such as those referenced in sayad_predictive_2019 ;
vilar_modelling_2021 ; valero_multifidelity_2021 ; mandel_wildland_2008  have
been undertaken to model and forecast forest fires utilizing wildland fire
models that can be simulated under various scenarios. Through this paper, we
present a novel approach for predicting the spread of forest fires through a
combination of dimensionality reduction and deep learning (DL).

The phenomena of wildland fires are usually governed by parametric non-linear
time-dependent partial differential equations (PDEs) with different scenarios
being achieved by changing physical parameters in the system of equations.
However, solving such equations for modeling large-scale wildland fires for
the numerous parameters is unaffordable in real-time emergency situations
requiring significant computing power and time. To speed up the numerical
simulations reduced-order modeling provides a promising strategy.

Reduced order modeling (ROM) for the parameterized PDEs usually relies on an
offline-online computational splitting benner_model_2017 . The expensive task
of building the low-dimensional subspace from the full order model (FOM)
snapshots is performed once in the so-called offline stage, and the reduced
order model (ROM) approximation corresponding to any new parameter value is
computed in the so-called online stage. Classical projection-based methods
such as proper orthogonal decomposition (POD) perform the dimensionality
reduction by constructing a low-dimensional trial subspace by the leading
modes of a singular value decomposition (SVD) and then projecting the system
dynamics onto the constructed low-dimensional subspace using Petrov-Galerkin
methods hesthaven_certified_2016 ; quarteroni_reduced_2016 . In many
applications, namely with faster decay of the Kolmogorov $n$-width
greif_decay_2019  (in turn the Hankel singular values unger_kolmogorov_2019 ),
a low-dimensional subspace with low approximation error if found with
conventional methods like POD. But the problem arises when the dimension of
the trial subspace needs to be large for a desired approximation quality. This
is almost always the case with transport-dominated fluid systems (TDFS) like
propagating flame fronts and traveling acoustic or shock waves krah_front_2022
; huang_challenges_2018 ; ohlberger_nonlinear_2013 .

To overcome the problem of slowly decaying singular values, model order
reduction (MOR) methods for TDFS have emerged over the last couple of years.
The work of rim_transport_2018  uses transport reversal where they do template
fitting by posing a minimization problem. Along similar lines
cagniart_model_2019  tries to transform/twist the set of solutions so that the
combination of a proper shift and appropriate linear combination recovers
accurate approximation. Similar approaches are followed in
welper_transformed_2020 ; black_projection-based_2020 . Transport maps were
explored in nonino_overcoming_2019 , approximating the field variable with a
front shape function and a level set function for efficient model reduction is
explored in krah_model_2021 ; krah_front_2022 . In general, most of these
methods rely on an offline-online splitting where the computationally
expensive offline step comprises a non-linear projection framework to capture
the high-dimensional space and the online step constructs the ROM
approximation intrusively for any unseen parameter value. Intrusive parameter
predictions have already been utilized for the wildland fire model in
black_efficient_2021 , where the shifted proper orthogonal decomposition
(sPOD) reiss_shifted_2018  is combined with residual minimization
lee_model_2020 ; black_projection-based_2020  to obtain an efficient ROM. More
powerful formulations of the sPOD are discussed in reiss_optimization-
based_2021 ; krah_non-linear_nodate . The sPOD has also been used in
mendible_dimensionality_2020  for predicting the detonation rotation waves
using UnTWIST algorithm and in papapicco_neural_2022  for the prediction of
shifts in nonlinear hyperbolic equations.

However, the non-linearity of the dimensional reduction for transport
phenomena hinders the projection of even a linear PDE on a small basis, and
thus it needs to be evaluated in the original dimension of the FOM in the
online stage. This prevents any significant speedup during the simulation. In
order to avoid scaling with the FOM dimension hyper-reduction techniques
sarna_hyper-reduction_2021 ; carlberg_gnat_2013 ; barrault_empirical_2004
have been developed, one such being discrete empirical interpolation method
(DEIM) jain_hyper-reduction_2017 ; chaturantabut_nonlinear_2010  which
approximates non-linear terms by evaluating them at a few, carefully selected
interpolation points and approximating all other components via interpolation
in a low-dimensional space. This method also serves as the starting point for
further extension of the method called shifted-DEIM by black_efficient_2021
for treating the wildland fire model.

Another class of methods called online adaptive basis methods
peherstorfer_model_2020 ; koch_dynamical_2007 ; etter_online_2020 ;
dihlmann_model_2011  have also been introduced to circumvent the problem of
slow-decaying Kolmogorov n-width by exploiting the time-local low-rank
structure of the TDFS. One such method called adaptive bases and adaptive
sampling discrete empirical interpolation method (AADEIM) was introduced by
Peherstorfer et.al (2020) in peherstorfer_model_2020 . It builds upon the
argument that the solutions of TDFS are typically low-dimensional if
considered locally in time. This notion is then exploited by approximating the
FOM solution in local low-dimensional spaces that are adapted via a low-rank
basis updates over time zimmermann_geometric_2018 ; peherstorfer_online_2015 .
However, the problem with these adaptive versions is that they usually do not
deliver speedups comparable to a traditional POD-DEIM approach, as shown for
shallow water flows in koellermeier_macro-micro_2023 . To circumvent these
problems, various non-intrusive methods have emerged in the recent past. These
methods rely in general on projecting the high-dimensional dynamics onto a
low-dimensional subspace usually by using the POD and then employing deep
learning (DL) framework for modeling the reduced dynamics.

A few attempts have been made to construct a low-dimensional subspace with
artificial neural networks (ANNs) gonzalez_deep_2018 ; lee_model_2020 ,
shallow masked encoder kim_fast_2022  and methods like kernel POD
salvador_non_2021  as well which allow for a non-linear representation of the
projected high-dimensional dynamics. Subsequently DL techniques such as feed-
forward neural network (FNN) based regression model san_neural_2018 ;
salvador_non_2021 ; hesthaven_non-intrusive_2018 ; wang_non-intrusive_2019 ;
berzins_standardized_2021  and Gaussian process regression kast_non-
intrusive_2020 ; guo_data-driven_2019  have been employed to model the reduced
dynamics. Model order reduction techniques based on a deep convolutional
autoencoder are proposed in gonzalez_deep_2018 ; bhattacharya_model_2021 .
Other DL techniques in the literature include long short-term memory (LSTMs)
mohan_deep_2018  and recurrent neural networks (RNNs) j_nagoor_kani_reduced-
order_2019 .

In our research, we expand upon the non-intrusive parameter prediction
framework introduced by Fresca et al. fresca_comprehensive_2021 ; fresca_pod-
dl-rom_2022 . This involves constructing a trial subspace using Proper
Orthogonal Decomposition (POD) and modeling the reduced dynamics using a
convolutional autoencoder coupled with a Feedforward Neural Network (FNN).
However, instead of the POD, we employ the sPOD for dimensionality reduction,
as it avoids the issue of slowly decaying singular values. On this base, we
employ a deep FNN to learn the time and parameter-dependent amplitudes and
shift transformations. This allows us to efficiently predict the states of the
system for unseen parameter values. A similar technique has been presented in
Kovarnova et al. kovarnova_shifted_2022 , where they only predict in time. We
apply our proposed approach to one- and two-dimensional wildland fire models,
and present computational results to demonstrate its effectiveness.

The structure of the paper is as follows. In Section 2 we lay out the
theoretical foundation concerning the mode decomposition methods, namely POD
and sPOD along with the formulation describing the data-driven approach. A
brief definition of the errors considered in our studies for model performance
evaluation is also mentioned in Section 2. In Section 3 we present the
numerical results for all the example test cases along with a timing study.
Finally, the conclusions are drawn in Section 4, and auxiliary tests are
provided in the Appendix.

##  2 Theory and background

In this section, we introduce the theoretical aspects of the presented
methods. We study the relevant mode decomposition approaches followed by the
data-driven approach. Furthermore, we scrutinize the types of errors that may
arise from the implementation of the proposed method.

###  2.1 Mode decomposition

We specifically consider POD and sPOD for our analysis.

####  2.1.1 Proper Orthogonal Decomposition (POD)

POD is a method to extract optimal basis sets from a collection of snapshots.
The snapshot matrix $Q\in\mathbb{R}^{M\times N}$ comprises of snapshots
$q(\bm{x}_{i},t_{j},\bm{\mu}_{j})$ arranged in a column-wise fashion for each
time step $t_{j}$ where $j=1,\ldots,N$, $\bm{x}_{i}\in\Omega$ are the spatial
grid points inside the domain $\Omega$ and $\bm{\mu}_{j}$ being the parameter
dependency. The POD method approximates with the help of SVD the snapshot
matrix $Q$

$Q\approx\tilde{Q}=U_{r}\Sigma_{r}(V_{r})^{\top}\,.$ (1)

Here, $r\ll\min(N,M)$ is the truncation rank,
$\Sigma_{r}=\mathrm{diag}{(\sigma_{1},\dots,\sigma_{r})}$ is a diagonal matrix
containing the singular values $\sigma_{1}\geq\sigma_{2}\geq\dots\sigma_{r}$
and $U_{r}\in\mathbb{R}^{M\times{r}}$, $V_{r}\in\mathbb{R}^{N\times{r}}$ are
orthogonal matrices containing the left and right singular vectors. The
reconstruction is optimal in the sense that the time-averaged least-square
error of the POD approximation

$\epsilon_{r}=\left\|Q-\tilde{Q}\right\|^{2}_{\mathrm{F}}=\sum^{\min(N,M)}_{k=r+1}\sigma_{k}$
(2)

is minimized.

####  2.1.2 Shifted proper orthogonal decomposition (sPOD)

The sPOD was introduced in reiss_shifted_2018  and further algorithmic
developments were presented in reiss_optimization-based_2021 ; krah_non-
linear_nodate  based on the optimization of singular values and in
black_projection-based_2020 ; black_efficient_2021  based on the optimization
of dyadic pairs and their shifts. The sPOD method aims to decompose the
snapshot matrix
${Q}=[q(\bm{x}_{i},t_{j},\bm{\mu}_{j})]_{ij}\in\mathbb{R}^{M\times N}$ into
multiple co-moving fields $\\{Q^{k}\in\mathbb{R}^{M\times N}\\}_{k=1,\dots,f}$
(i.e. $\\{q^{k}(\bm{x}_{i},t_{j},\bm{\mu}_{j})\\}_{k=1,\dots,f}$). The
decomposition follows as:

$\displaystyle q(\bm{x},t,\bm{\mu})$
$\displaystyle\approx{\tilde{q}(\bm{x},t,\bm{\mu})}:=\sum_{k=1}^{f}\mathcal{T}^{{\Delta}^{k}}q^{k}(\bm{x},t,\bm{\mu})\,,\quad$
(continuous) (3) $\displaystyle Q$
$\displaystyle\approx{\tilde{Q}}:=\sum_{k=1}^{f}T^{\bm{\Delta}^{k}}Q^{k}\,,\quad$
(discrete) (4)

The interpolation based discrete transformation operators are given by
$\\{T^{\bm{\Delta^{k}}}\\}_{k=1,\dots,f}$ where
$\bm{\Delta}^{k}(t,\mu)={(\Delta^{k}(\bm{x}_{1},t,\bm{\mu}),\dots,\Delta^{k}(\bm{x}_{M},t,\bm{\mu}))}$
is the time dependent shift such that
$T^{\pm\bm{\Delta}^{k}}(Q)_{ij}=q(\bm{x}_{i}\mp{\Delta}^{k}(\bm{x}_{i},t,\bm{\mu})_{i},t_{j},\bm{\mu}_{j})$,
and $\tilde{Q}$ is the approximate reconstruction of $Q$. The assumption is
that for traveling wave systems the superposition Equation 4 can decompose the
data more efficiently than the POD. This is because the traveling wave is only
changing slowly in the co-moving data frame and can therefore be well
approximated with only a few modes calculated with a truncated SVD:

$\displaystyle Q^{k}\approx
U_{r^{k}}^{k}\Sigma_{r^{k}}^{k}(V_{r^{k}}^{k})^{\top}\,\quad k=1,\dots,f\,.$
(5)

Here, $r_{k}\ll N$ is the truncation rank of each co-moving field,
$\Sigma^{k}_{r_{k}}=\mathrm{diag}{(\sigma_{1}^{k},\dots,\sigma_{r_{k}}^{k})}$
is a diagonal matrix containing the singular values
$\sigma_{1}^{k}\geq\sigma_{2}^{k}\geq\dots\sigma_{r_{k}}^{k}$ and
$U_{r_{k}}^{k}\in\mathbb{R}^{M\times{r_{k}}}$,
$V_{r_{k}}^{k}\in\mathbb{R}^{N\times{r_{k}}}$ are orthogonal matrices
containing the left and right singular vectors. In this work we use the sPOD
algorithm based on the minimization of the nuclear-one norm, which is
presented in krah_non-linear_nodate  algorithm 8, and the corresponding tuning
parameters for the algorithm are stated in the Appendix C.

######  Remark 1.

In order to be more efficient than the POD, when disregarding additional
degrees of freedom introduced by shifts, we assume that for the total number
of degrees of freedom $r=\sum_{k=1}^{f}r_{k}$ the truncation error of the sPOD
is smaller than the POD truncation error eq. 2.

###  2.2 Data-driven approach

We consider a non-linear parameterized dynamical system which often stems from
a discretized PDE system like:

$\displaystyle\bm{\dot{q}}(t,\bm{\mu})$
$\displaystyle=\bm{f}(t,\bm{q}(t,\bm{\mu}),\mu)\quad t\in[0,T]$ (6)
$\displaystyle\bm{q}(0,\bm{\mu})$ $\displaystyle=\bm{q_{0}}(\bm{\mu})$

where, $\bm{q}\in\mathbb{R}^{M}$ is the parameterized solution to the problem,
$\bm{q_{0}}$ is the initial data, $\bm{f}$ is the non-linear function
describing the system dynamics and $\bm{\mu}$ is the parameter vector. For
simplicity we only consider a single parameter $\mu$ for our study, however,
real-world scenarios may depend upon multiple parameters.

####  2.2.1 Non-intrusive predictions with POD and sPOD

To develop the new method, we closely follow the POD-DL-ROM method described
by Fresca et. al fresca_pod-dl-rom_2022 . The said method first builds a low-
dimensional subspace $Col(U_{r})$ which is spanned by the first $r$ singular
vectors of the parameter snapshot matrix $Q$:

$Q=\left[\begin{array}[]{ccccccc}|&&|&&|&&|\\\
\bm{q}(t_{1},\mu_{1})&\ldots&\bm{q}(t_{N_{t}},\mu_{1})&\ldots&\bm{q}(t_{1},\mu_{N_{p}})&\ldots&\bm{q}(t_{N_{t}},\mu_{N_{p}})\\\
|&&|&&|&&|\end{array}\right]\in\mathbb{R}^{M\times N}.$ (7)

The matrix $Q$ is the collection of $N=N_{p}N_{t}$ number of FOM snapshots
computed for different parameter instances $\mu_{1},\ldots,\mu_{N_{p}}$ that
are sampled over different time instances $t_{1},\ldots,t_{N_{t}}$. Using POD
on $Q$ and truncating to a rank $r\ll M$, we have

$\tilde{Q}=U_{r}\Sigma_{r}(V_{r})^{\top}\ \quad\text{where},\quad
A=\Sigma_{r}(V_{r})^{\top}\ $

and

$A=\left[\begin{array}[]{ccccccc}|&|&&&|\\\
\bm{a}(t_{1},\mu_{1})&\bm{a}(t_{2},\mu_{1})&\ldots&\ldots&\bm{a}(t_{N_{t}},\mu_{N_{p}})\\\
|&|&&&|\end{array}\right]\in\mathbb{R}^{r\times N}.$ (8)

A deep learning model is used to approximate the mapping
$(t_{1},\mu_{1})\to\bm{a}$, by using $A$ as training data along with $P$:

$P=\left[\begin{array}[]{ccccccc}(t_{1},\mu_{1})&\ldots&(t_{N_{t}},\mu_{1})&\ldots&(t_{1},\mu_{N_{p}})&\ldots&(t_{N_{t}},\mu_{N_{p}})\end{array}\right]\in\mathbb{R}^{(n_{\mu}+1)\times
N}\,.$ (9)

The matrix $P$ corresponds to $n_{\mu}$ number of parameters along with the
contribution of time $t$ as an additional parameter thus making
$P\in\mathbb{R}^{(n_{\mu}+1)\times N}$. By this, the network learns the
mapping from parameter and time to amplitude, which is the mapping
traditionally constructed by the Galerkin approach. Subsequently, once the
network is trained, the prediction step is performed in order to generate the
time amplitude for any $(t,\mu)$

$\displaystyle\bm{a}(t,\mu)=\left[\begin{array}[]{c}a_{1}(t,\mu)\\\
a_{2}(t,\mu)\\\ \vdots\\\
a_{r}(t,\mu)\end{array}\right]\approx\mathbf{N_{\mathrm{POD}}}(t,\mu)$ (14)

and the state $\bm{\tilde{q}^{\prime}}$ for those unseen parameters is
reconstructed as:

$\bm{\tilde{q}^{\prime}}(t,\mu)=U_{r}\bm{a}(t,\mu)\,.$ (15)

More in-depth analysis of the POD-DL-ROM method along with the network
architectures, working examples, and the error study can be found in
fresca_comprehensive_2021 ; fresca_pod-dl-rom_2022 .

Since the POD is not suitable for transport-dominated systems, we present a
novel non-intrusive model order reduction technique named sPOD-NN that uses
sPOD for constructing the desired low-dimensional subspace. The sPOD is
applied on the parameter snapshot matrix $Q$ from Equation 7. The sPOD
decomposes the snapshot matrix $Q$ into co-moving frames $Q^{k}$. The
algorithm also outputs the basis vectors $U^{k}_{r^{k}}$ for every frame with
which the time amplitude matrix $A^{k}$ is extracted for all $k=1,\ldots,f$.
Following Equation 5 the time amplitudes are defined,

$A^{k}=(U^{k}_{r^{k}})^{\top}\ Q^{k}\in\mathbb{R}^{r^{k}\times N},$ (16)

where

$A^{k}=\left[\begin{array}[]{ccccccc}|&|&&&|\\\
\bm{a}^{k}(t_{1},\mu_{1})&\bm{a}^{k}(t_{2},\mu_{1})&\ldots&\ldots&\bm{a}^{k}(t_{N_{t}},\mu_{N_{p}})\\\
|&|&&&|\end{array}\right].$ (17)

We now have the time amplitude matrices $A^{k}$ available for all $k$ along
with the shifts. However, we note that the shifts for any new unseen parameter
set $(t,\mu)$ are not available upfront. This poses a challenge because, in
due course of reconstructing the final state for unseen parameter values, the
shifts must be computed first.

##### Shifts:

The shifts
$\bm{\Delta}^{k}(t,\mu)=(\Delta^{k}(\bm{x}_{1},t,\mu),\dots,\Delta^{k}(\bm{x}_{M},t,\mu))\in\mathbb{R}^{M}$,
used to encode the transport, are part of the description. They are in general
dependent on the model parameters $\mu$, time $t$ and space $\bm{x}$. The
flame propagation speed changes with the Arrhenius factor for the wildland
fire model, and it needs to be predicted by the neural network. To have a low-
dimensional description, the description of the shifts needs to be low-
dimensional as well. We basically study two possible scenarios, where the
first can be seen as a special case of the second:

##### (a) Low-dimensional shifts:

We look at a case where the shifts are independent of the spatial coordinates
$\bm{x}\in\mathbb{R}^{d}$:
${\Delta}^{k}(\bm{x},t,\mu)=\sum_{i=1}^{n_{\Delta}^{k}}\bm{e}_{i}\>\underline{\Delta}^{k}_{i}(t,\mu)\,,\quad\text{where}\quad
n_{\Delta}^{k}\leq d\,.$ (18) Here, $\bm{e}_{i}$ denotes the $i$th standard
basis vector. For our one- and two-dimensional examples $d=1,2$. Since
$n_{\Delta}^{k}\leq d$ the coefficient vectors of the shifts
$\underline{\bm{\Delta}}^{k}(t,\mu)=(\underline{\Delta}_{1}^{k}(t,\mu),\dots,\underline{{\Delta}}^{k}_{n_{\Delta}^{k}}(t,\mu))$
are already low-dimensional and only depend on time and parameter.

##### (b) Low-rank description of high-dimensional shifts:

The sPOD usually assumes a low-dimensional description of the shifts. However,
for complicated systems, the shifts might depend on the spatial position
$\bm{x}$ itself and are thus high-dimensional. Nevertheless, the shifts often
have a low-rank structure and we thus assume, they can be well represented
with the help of the POD:
${\Delta}^{k}(\bm{x},t,\mu)\approx\sum_{n=1}^{n_{\Delta}^{k}}\Upsilon_{n}^{k}(\bm{x})\underline{\Delta}_{n}^{k}(t,\mu)\,,\quad\text{with}\quad
n_{\Delta}^{k}\ll N\,.$ (19) Thus after we assemble the shift matrix
$[\Delta^{k}(\bm{x}_{i},t_{j},\mu_{j})]_{ij}\in\mathbb{R}^{M\times N}$ using a
threshold algorithm (see section 3.1.3), we decompose it with the help of the
truncated SVD and obtain the low-dimensional shifts
$\underline{\bm{\Delta}}^{k}(t,\mu)=(\underline{\Delta}_{1}^{k}(t,\mu),\dots,\underline{{\Delta}}^{k}_{n_{\Delta}^{k}}(t,\mu))$.
Once we obtain the $\underline{\bm{\Delta}}^{k}(t,\mu)$ we construct the
training data by stacking all $A^{k}$ and $\underline{\bm{\Delta}}^{k}(t,\mu)$
to construct an
$\hat{A}\in\mathbb{R}^{\left(\sum_{k}r^{k}+\sum_{k}n_{\Delta}^{k}\right)\times
N}$ matrix as shown:

$\hat{A}=\left[\begin{array}[]{ccccc}\bm{a}^{1}(t_{1},\mu_{1})&\bm{a}^{1}(t_{2},\mu_{1})&\ldots&\ldots&\bm{a}^{1}(t_{N_{t}},\mu_{N_{p}})\\\
\vdots&\vdots&&&\vdots\\\
\bm{a}^{f}(t_{1},\mu_{1})&\bm{a}^{f}(t_{2},\mu_{1})&\ldots&\ldots&\bm{a}^{f}(t_{N_{t}},\mu_{N_{p}})\\\
\underline{\bm{\Delta}}^{1}(t_{1},\mu_{1})&\underline{\bm{\Delta}}^{1}(t_{2},\mu_{1})&\ldots&\ldots&\underline{\bm{\Delta}}^{1}(t_{N_{t}},\mu_{N_{p}})\\\
\vdots&\vdots&&&\vdots\\\
\underline{\bm{\Delta}}^{f}(t_{1},\mu_{1})&\underline{\bm{\Delta}}^{f}(t_{2},\mu_{1})&\ldots&\ldots&\underline{\bm{\Delta}}^{f}(t_{N_{t}},\mu_{N_{p}})\\\
\end{array}\right]\,.$ (20)

As for training the deep learning model, we consider $\hat{A}$ and the entries
of $P$ matrix from Equation 9 as the training data. After the successful
training of the network, the time amplitudes and the shifts are predicted for
$(t,\mu)$

$\displaystyle\bm{\hat{a}}(t,\mu)=\left[\begin{array}[]{c}\bm{a}^{1}(t,\mu)\\\
\vdots\\\ \bm{a}^{f}(t,\mu)\\\ \underline{\bm{\Delta}}^{1}(t,\mu)\\\ \vdots\\\
\underline{\bm{\Delta}}^{f}(t,\mu)\end{array}\right]\approx\mathbf{N_{\mathrm{sPOD}}}(t,\mu)\,.$
(27)

By reconstructing the shifts using Equation 18 or eq. 19 depending on the
problem at hand the state $\bm{\tilde{q}^{\prime}}$ for the unseen parameters
can be reconstructed as

$\bm{\tilde{q}^{\prime}}(t,\mu)=\sum^{f}_{k=1}T^{\bm{\Delta}^{k}(t,\mu)}\cdot\left(U^{k}_{r^{k}}\bm{a}^{k}(t,\mu)\right)\,.$
(28)

For later comparison studies, we define the degrees of freedom
($n_{\mathrm{dof}}$) as:

$n_{\mathrm{dof}}=\begin{cases}r,&\text{for POD-NN}\\\
\sum_{k}r^{k}+\sum_{k}n_{\Delta}^{k},&\text{for sPOD-NN}\end{cases}$ (29)

####  2.2.2 Interpolation method with sPOD

For comparison we introduce another approach based on sPOD which we call
sPOD-I (Interpolation) which instead of using the DL techniques uses an
interpolation-based method in the online phase. The sPOD-I method first
extracts the $P$ and $\hat{A}$ matrices as explained in Equation 9 and
Equation 20 respectively. Subsequently, it then uses the scipy
virtanen_scipy_2020  function scipy.interpolate.griddata() for interpolation.
This function triangulates the input domain orourke_computational_1998  and
performs barycentric interpolation on each triangle to construct an
interpolant. For an interpolation point $(t,\mu)$ lying inside a triangle the
interpolated value is given as

$\bm{\hat{a}}(t,\mu)=\sum^{3}_{i=1}\alpha_{i}\bm{\hat{a}}(t_{i},\mu_{i}),$
(30)

where $\alpha_{i}$ are barycentric coordinates with $\sum\alpha_{i}=1$. Once
the time amplitudes and the shifts are obtained, the final state can be
reconstructed as shown in Equation 28.

####  2.2.3 Network architecture

We employ a deep FNN for the predictions. We use PyTorch paszke_pytorch_2019
for constructing the neural networks. For the architecture, we have an input
layer, three hidden layers, and an output layer. The parameters of the
network, the number of inputs and outputs which necessarily are the neurons
for each layer are shown in Table 1. In the table, we have $p$ as the input
for the input layer of the network where $p=n_{\mu}+1$ as described in
Equation 9. For our examples, we only consider a single parameter $\mu$ thus
$p=2$. The number of outputs $u$ in the output layer is problem dependent and
has been pointed out in the respective numerical examples.

Table 1: Deep feed-forward neural network parameters Layer | Num. of Inputs | Num. of Outputs | Activation function
---|---|---|---
Input layer | $p$ | 25 | ELU
$1^{st}$ hidden layer | 25 | 50 | ELU
$2^{nd}$ hidden layer | 50 | 75 | ELU
$3^{rd}$ hidden layer | 75 | 50 | LeakyReLU
Output layer | 50 | $u$ | -

We basically rely on the L1 loss function, namely MAE (Mean Absolute Error):

$\mathcal{L}_{\mathrm{MAE}}=\begin{cases}\frac{1}{n_{\mathrm{dof}}}\left\|\bm{a}(t,\mu)-{\textbf{N}_{\mathrm{POD}}}(t,\mu)\right\|_{1},&\text{for
POD-NN}\\\
\frac{1}{n_{\mathrm{dof}}}\left\|\bm{\hat{a}}(t,\mu)-{\textbf{N}_{\mathrm{sPOD}}}(t,\mu)\right\|_{1},&\text{for
sPOD-NN},\end{cases}$ (31)

where $n_{\mathrm{dof}}$ is defined in eq. 29. For improving the prediction
accuracy we use data scaling. In this work, we use Min-Max scaling to scale
our training data to the interval $[0,1]$. Let us consider $\hat{A}$ from
Equation 20 and $P$ from Equation 9 which serve as the training data for our
network. The scaling for the matrix $P$ is given as:

$P^{\mathrm{scaled}}_{ij}=\frac{P_{ij}-\underset{j=1,\ldots,N_{s}}{\min}(P_{ij})}{\underset{j=1,\ldots,N_{s}}{\max}(P_{ij})-\underset{j=1,\ldots,N_{s}}{\min}(P_{ij})},\quad\text{for}\quad
i=1,\ldots,n_{\mu}+1$ (32)

whereas, for the matrix $\hat{A}$ where the time amplitudes and the shifts are
stacked together, the scaling is performed separately for both quantities as
shown here:

$\displaystyle\hat{A}^{\mathrm{scaled}}=\left[\begin{array}[]{c}\frac{\bm{A}_{ij}-\underset{i,j}{\min}(\bm{A}_{ij})}{\underset{i,j}{\max}(\bm{A}_{ij})-\underset{i,j}{\min}(\bm{A}_{ij})}\\\
\vspace{0.05cm}\hfil\\\
\frac{\underline{\bm{\Delta}}_{ij}-\underset{i,j}{\min}(\underline{\bm{\Delta}}_{ij})}{\underset{i,j}{\max}(\underline{\bm{\Delta}}_{ij})-\underset{i,j}{\min}(\underline{\bm{\Delta}}_{ij})}\end{array}\right]$
(36)

The minimum and maximum values for all calculated quantities are stored and
then used to rescale the prediction results back to their original scale.
Neural network basics are explained in Appendix A in more detail.

###  2.3 Errors

Before we look at the numerical results it is crucial to define the errors
incurred by the aforementioned approaches. Once the network is trained and the
time amplitudes and the shifts are predicted, the errors are as follows:
POD-NN:

$E^{\mathrm{POD}}=\frac{\left\lVert
Q-\tilde{Q}\right\rVert_{\mathrm{F}}}{\left\lVert
Q\right\rVert_{\mathrm{F}}},\quad\text{where}\quad\tilde{Q}=U_{r}\Sigma_{r}(V_{r})^{\top}\
$ (37) $E^{\mathrm{POD-NN}}_{\mathrm{tot}}=\frac{\left\lVert
Q-U_{r}A_{\mathrm{NN}}\right\rVert_{\mathrm{F}}}{\left\lVert
Q\right\rVert_{\mathrm{F}}}$ (38)

The error $E^{\mathrm{POD-NN}}_{\mathrm{tot}}$ is bounded by the POD
truncation error $E^{\mathrm{POD}}$ and the neural network prediction error.
sPOD-NN:

$E^{\mathrm{sPOD}}=\frac{\left\lVert
Q-\tilde{Q}\right\rVert_{\mathrm{F}}}{\left\lVert
Q\right\rVert_{\mathrm{F}}},\quad\text{where}\quad\tilde{Q}\approx\sum^{f}_{k=1}T^{\bm{\Delta}^{k}}Q^{k},$
(39) $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}=\frac{\left\lVert
Q-\sum_{k=1}^{f}T^{\bm{\Delta}^{k}_{\mathrm{NN}}}\left(U^{k}_{r^{k}}A^{k}_{\mathrm{NN}}\right)\right\rVert_{\mathrm{F}}}{\left\lVert
Q\right\rVert_{\mathrm{F}}}$ (40)

The final error $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}$ is bounded by the sPOD
truncation error $E^{\mathrm{sPOD}}$, the neural network prediction error and
the error which arises from the shift operations.
sPOD-I:

$E^{\mathrm{sPOD-I}}_{\mathrm{tot}}=\frac{\left\lVert
Q-\sum_{k=1}^{f}T^{\bm{\Delta}^{k}_{\mathrm{I}}}\left(U^{k}_{r^{k}}A^{k}_{\mathrm{I}}\right)\right\rVert_{\mathrm{F}}}{\left\lVert
Q\right\rVert_{\mathrm{F}}}$ (41)

The final error $E^{\mathrm{sPOD-I}}_{\mathrm{tot}}$ is bounded by the sPOD
truncation error, the interpolation error, and the error due to shift
operations.

We also look at the error calculated at every time instance fresca_pod-dl-
rom_2022  to gain a comprehensive understanding of the performance of the
method. Although this can be computed for all the three methods explained, we
only focus on the sPOD-NN method for our analysis. For $(t_{j},\mu)$ with
Equation 28 we have

$E^{\mathrm{sPOD-
NN}}_{j}=\frac{\left\lvert\bm{q}(t_{j},\mu)-\bm{\tilde{q}^{\prime}}(t_{j},\mu)\right\rvert}{\sqrt{\frac{1}{N_{t}}\sum^{N_{t}}_{j=1}\left\lVert\bm{q}(t_{j},\mu)\right\rVert^{2}_{\mathrm{F}}}}\,$
(42)

which basically gives us the error at each grid point for a particular
$(t_{j},\mu)$.

##  3 Numerical results

In order to assess the effectiveness of the sPOD-NN method, we evaluate its
numerical performance 111The source code is available at
https://github.com/MOR-transport/sPOD-NN-paper by the run time and the errors
outlined in Section 2. For a fair comparison, we also compare it with the POD-
NN and sPOD-I methods, based on the aforementioned errors. The techniques are
evaluated on 1D and 2D wildland fire models. For sPOD we use the algorithm 8
in krah_non-linear_nodate . The tuning parameters can be found in Appendix C.

For benchmarking, the methods were also checked on synthetically generated
test data for which the results can be found in Appendix D. The importance of
the test case can be inferred from the left plot in Figure 25 of Appendix D
where we see that $E^{\mathrm{sPOD-I}}_{\mathrm{tot}}$ decays with the sPOD
basis reconstruction error $E^{\mathrm{sPOD}}$ with an increasing
$n_{\mathrm{dof}}$. This implies that if the online phase is substantially
accurate then the final error will be dominated only by the basis
reconstruction error. This is certainly true for the sPOD-I but not so for the
other methods for the given test case. This is by the design of the problem,
where the time amplitudes belong to a polynomial space and the method sPOD-I
interpolates using polynomials.

###  3.1 Wildland fire model

We use the model shown in mandel_wildland_2008 . Considering a domain
$\Omega\in\mathbb{R}^{d}$ with $d=1,2$ and a finite time horizon
$\mathbb{T}:=[0,t_{f}]$ where $0<t_{f}<\infty$, we compute the temperature $T$
and the fuel supply mass fraction $S$ satisfying the coupled nonlinear PDE

$\displaystyle\partial_{t}T$ $\displaystyle=\nabla\cdot(k\nabla
T)-v\cdot\nabla T+\alpha(Sr(T,\mu,T_{a})-\gamma(T-T_{a})),$ (43)
$\displaystyle\partial_{t}S$ $\displaystyle=-S\gamma_{S}r(T,\mu,T_{a}),$

where $\gamma_{S}r(T,\mu,T_{a})$ is called the reaction rate constant which is
given by the modified Arrhenius law:

$r(T,\mu,T_{a}):=\begin{cases}\mathrm{exp}(-\frac{\mu}{T-T_{a}}),&\text{if
$T>T_{a}$},\\\ 0,&\text{otherwise}\,.\end{cases}$ (44)

Variables and their fixed values used in the model are shown in Table 2.

Table 2: Wildland fire model variables and values for the fixed coefficients involved Name | Symbol | Unit | Value
---|---|---|---
temperature | $T$ | $\mathrm{K}$ | -
supply mass fraction | $S$ | - | -
thermal diffusivity | $k$ | $\mathrm{m^{2}/s}$ | 0.2136
Arrhenius coefficient111Tuning factor | $\mu$ | $\mathrm{K}$ | 558.49
pre-exponential factor | $\gamma_{S}$ | $\mathrm{1/s}$ | 0.1625
wind speed | $v$ | $\mathrm{m/s}$ | 0
temperature rise per second | $\alpha$ | $\mathrm{K/s}$ | 187.93
scaled heat transfer coefficient | $\gamma$ | $\mathrm{1/K}$ | $4.8372\times 10^{-5}$
ambient temperature | $T_{a}$ | $\mathrm{K}$ | 300
|  |  |

For notational convenience, we work with the relative temperature $T-T_{a}$
instead of the temperature which results in $T_{a}=0$ for the rest of the
study.

####  3.1.1 1D model

Here, we will motivate the application of the proposed method to perform non-
linear model order reduction of a one-dimensional wildland fire model. As an
initial step, we solve the model equations Equation 43 and Equation 44 on a
one-dimensional strip of length $l=1000\mathrm{m}$ in a computational domain
$\Omega=[0,1000]$ for $t_{f}=1400\mathrm{s}$ yielding the time domain
$\mathbb{T}=[0,1400]$. The velocity of wind $v=0$. The time domain is chosen
such that the traveling fronts never reach the boundaries. We consider the
initial condition to be:

$T^{0}(x):=1200\mathrm{exp}\left(-\frac{(x-l/2)^{2}}{200m^{2}}\right)\quad\text{and}\quad
S^{0}\equiv 1.$ (45)

The value of Arrhenius coefficient $\mu=558.49\mathrm{K}$ is used for our
simulation. A point to note is that we use $\mu$ as our parameter of choice
for tuning the model. A lower $\mu$ results in a faster spread of the fire and
a higher $\mu$ results in a slowly expanding fire. Periodic boundary
conditions are used for simplicity. The spatial domain is decomposed into
$N_{x}=3000$ grid points making $h_{x}=1/3\mathrm{m}$ and discretized with
finite difference method with $6^{th}$ order central finite difference
stencil. We solve the model for $6000$ time steps. Time integration is done
with the standard RK4 method for which $dt$ is calculated with the help of CFL
criteria:

$dt=\mathrm{cfl}\cdot\frac{\sqrt{h_{x}^{2}+h_{y}^{2}}}{c}$ (46)

where we prescribe $\mathrm{cfl}=0.7$, $h_{y}=0$ and speed of sound $c=1$.

Subsequently, the obtained temperature and the supply mass fraction profiles
are shown in Figure 1. In the temperature plot we observe, a localized
ignition region from which flame fronts emanate and travel in opposite
directions. This is a perfect example of traveling wave-type phenomena. The
supply mass fraction varies in the range $[0,1]$ with the region left and
right of the two traveling fronts respectively being still unburnt while the
central region is completely burnt.

Figure 1: Full order model : Temperature(T) and the supply mass fraction(S)
profiles computed for $\mu=558.49\mathrm{K}$

For constructing non-linear parametric ROMs we solve Equation 43 and Equation
44 for
$\mu\in\\{540\mathrm{K},550\mathrm{K},560\mathrm{K},570\mathrm{K},580\mathrm{K}\\}$
and sample naively every $10^{th}$ snapshot in time to construct a parametric
snapshot matrix $Q$. More advanced sampling strategies can be found in
karcher_adaptive_2022 . With $N_{t}=6000/10=600$ time instances, $N_{p}=5$
parameter instances of $\mu$ and $M=N_{x}=3000$ the constructed matrix
$Q\in\mathbb{R}^{M\times N_{t}N_{p}}=\mathbb{R}^{3000\times 3000}$. Once the
parameter snapshot matrix is assembled we perform dimensionality reduction by
sPOD and obtain a low-dimensional representation of the data as shown in
Figure 2. The number of co-moving frames is decided by the user and for this
example it is chosen as $3$ (one for each of the left and right traveling
fronts and one for the initial ignition regime). The shifts for the frames
required for sPOD are calculated by tracking the gradients of the field
variable as explained in black_efficient_2021 . We note here that we only show
the results for the temperature field for simplicity, however, the entire
procedure can be replicated for the supply mass fraction as well, the results
of which could be found in Appendix B. The motivation for using sPOD is
substantiated by the basis reconstruction error results shown in Table 3. We
observe that as the number of modes increases the drop in $E^{\mathrm{sPOD}}$
is significant in comparison to $E^{\mathrm{POD}}$. For eg. just for
$16+10+16=42$ modes we could get to
$E^{\mathrm{sPOD}}\sim\mathcal{O}(10^{-4})$ whereas
$E^{\mathrm{POD}}\sim\mathcal{O}(10^{-2})$.

Figure 2: sPOD decomposition of the temperature field. The second, third, and fourth plots show the co-moving frames into which the field is decomposed. The transformation operators $T^{\bm{\Delta}^{k}}$ transform the stationary decomposed frames into the co-moving frames which when combined gives the reconstructed snapshot shown in the first plot  $Q=T^{\bm{\Delta}^{1}}Q^{1}+T^{\bm{\Delta}^{2}}Q^{2}+T^{\bm{\Delta}^{3}}Q^{3}$ Table 3: Offline and online error study w.r.t number of modes for 1D wildland fire model Modes111We also have $r=r_{1}+r_{2}+r_{3}+2$ where $r_{1}$ and $r_{3}$ are the ranks of the co-moving frames and $r_{2}$ is the rank of the stationary frame and $2$ accounts for two additional degrees of freedom controlling the shifts. For the comparison, we, therefore, use $r$ modes for POD. | Offline errors | Online errors
---|---|---
$r_{1}+r_{2}+r_{3}$ | $E^{\mathrm{sPOD}}$ | $E^{\mathrm{POD}}$ | $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}$ | $E^{\mathrm{sPOD-I}}_{\mathrm{tot}}$ | $E^{\mathrm{POD-NN}}_{\mathrm{tot}}$
1 + 1 + 1 | $3.85\times 10^{-1}$ | $6.90\times 10^{-1}$ | $3.93\times 10^{-1}$ | $2.36\times 10^{-1}$ | $6.92\times 10^{-1}$
2 + 2 + 2 | $1.97\times 10^{-1}$ | $5.03\times 10^{-1}$ | $2.34\times 10^{-1}$ | $2.31\times 10^{-1}$ | $5.03\times 10^{-1}$
4 + 2 + 4 | $6.13\times 10^{-2}$ | $3.57\times 10^{-1}$ | $1.53\times 10^{-1}$ | $7.83\times 10^{-2}$ | $3.55\times 10^{-1}$
6 + 2 + 6 | $2.49\times 10^{-2}$ | $2.69\times 10^{-1}$ | $1.00\times 10^{-1}$ | $4.05\times 10^{-2}$ | $2.67\times 10^{-1}$
8 + 4 + 8 | $1.26\times 10^{-2}$ | $1.86\times 10^{-1}$ | $7.25\times 10^{-2}$ | $3.83\times 10^{-2}$ | $1.85\times 10^{-1}$
10 + 4 + 10 | $5.51\times 10^{-3}$ | $1.48\times 10^{-1}$ | $4.09\times 10^{-2}$ | $3.67\times 10^{-2}$ | $1.48\times 10^{-1}$
14 + 8 + 14 | $1.20\times 10^{-3}$ | $8.07\times 10^{-2}$ | $3.99\times 10^{-2}$ | $3.64\times 10^{-2}$ | $8.67\times 10^{-2}$
16 + 10 + 16 | $6.94\times 10^{-4}$ | $6.06\times 10^{-2}$ | $3.78\times 10^{-2}$ | $3.64\times 10^{-2}$ | $7.85\times 10^{-2}$
|  |  |  |  |

As a result of the dimensionality reduction we obtain the time amplitude
matrix $A^{k}$ from $Q^{k}$ as shown in Equation 16 and we already have the
shifts calculated for all $k$. Following the procedure outlined in Section 2
for training the neural network model we assemble the time amplitude matrix
$\hat{A}\in\mathbb{R}^{42\times 3000}$ and the parameter matrix
$P\in\mathbb{R}^{2\times 3000}$. We only consider $2$ shifts in the
aforementioned assembly of the matrix $\hat{A}$, as the shift for the second
frame $\underline{\bm{\Delta}}^{2}=0$. The network is trained for
$N_{\mathrm{epochs}}=100000$ with batch size $N_{b}=100$. Early stopping
criteria are imposed to prevent overfitting which stops the training if the
validation loss does not decrease for 3000 consecutive epochs. For testing, we
choose $\mu=[558.49\mathrm{K}]$ and use the trained model to predict the time
amplitudes and the shifts which are shown in Figure 3. For a more quantitative
picture, we look at the error estimates shown in Figure 4.

Figure 3: Plot (a) and (b) show the first modes for the frame 1 and 2
respectively as frame 3 is identical to its counterpart frame 1. (c) and (d)
show the shifts for frame 1 and frame 3 respectively. (e) and (f) show the
first two modes corresponding to the POD-NN approach. Figure 4: The first plot
shows the error decay for an increasing $n_{\mathrm{dof}}$. The errors labeled
$E_{\mathrm{tot}}$ show the final reconstruction error after prediction or
interpolation. The second plot shows the trend of the relative error
$E^{\mathrm{sPOD-NN}}_{j}$ over time.

In the first plot, we see that as the $n_{\mathrm{dof}}$ increases all three
errors go down as expected. However, for both sPOD-NN and sPOD-I the
$E_{\mathrm{tot}}$ starts to stagnate after a certain $n_{\mathrm{dof}}$. We
observe that when we increase the $n_{\mathrm{dof}}$ the amount of new
information added to the training data is not substantial. On the other hand,
the parameter matrix $P$ remains the same. Thus in turn the approximation
capability of the network remains almost the same even while making the
prediction task more and more difficult. This also seems to be the case for
sPOD-I. For POD-NN this effect can not be fully captured here although, slight
stagnation could be seen towards the end of the curve. This effect is also
mentioned in Fresca et.al (2021) fresca_pod-dl-rom_2022 . sPOD-NN is able to
reach to an $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}\sim 0.037$. We also study the
error $E^{\mathrm{sPOD-NN}}_{j}$ which is shown in the right plot of Figure 4
where we see $\max\sim\mathcal{O}(10^{-2})$ whereas mean
$<\mathcal{O}(10^{-3})$.

The full reconstructed results are plotted in Figure 5. We select two random
time instances: one near the ignition and the other nearing the end of the
time domain and the prediction results are shown for all the methods at these
instances. We could observe that the POD-NN although being able to have
$E^{\mathrm{POD-NN}}_{\mathrm{tot}}\sim 0.078$, has spurious oscillations. The
$E^{\mathrm{POD}}\sim 0.06$ thus we infer that the prediction accuracy for
POD-NN is limited by the basis reconstruction error itself, which cannot be
resolved unless more modes are added for the study. As for sPOD-NN and sPOD-I,
we see minor oscillations near the sharp edges of the curve.

Figure 5: This figure captures the prediction accuracy of the studied methods.
We study the cross-sectional views of the temperature profile at $t=35s$ and
$1190s$.

####  3.1.2 2D model (without wind)

The model equations for the wildland fire can also be extended to a 2D case.
We solve the equations Equation 43 and Equation 44 on a two-dimensional square
of side length $l=500$ in a computational domain $\Omega=[0,500]\times[0,500]$
for $t_{f}=1000s$, yielding the time domain $\mathbb{T}=[0,1000]$ such that
the traveling fronts do not reach the boundaries. The velocity of wind
$v=(0,0)$ results in no change in the topology of the fronts. The initial
condition is given to be:

$T^{0}(x,y):=1200\mathrm{exp}\left(-\left(\frac{(x-l/2)^{2}}{200m^{2}}+\frac{(y-l/2)^{2}}{200m^{2}}\right)\right)\quad\text{and}\quad
S^{0}(x,y)\equiv 1.$ (47)

Periodic boundary conditions are employed. The spatial domain is split into
$N_{x}=500,N_{y}=500$ grid points making $h_{x}=1,h_{y}=1$ and discretized
with finite difference method with $6^{th}$ order central finite difference
stencil. We solve the system for $1000$ time steps. Time integration is done
with standard RK4 scheme and the $dt$ is computed by CFL condition shown in
Equation 46 where we prescribe $\mathrm{cfl}=0.7$ and $c=1$. Similar to the 1D
case we consider $\mu=558.49\mathrm{K}$ which also acts as a tuning parameter
in the model reduction procedure. The temperature and the supply mass fraction
profiles are shown in Figure 6.

Figure 6: Full order Model: Temperature(T) and supply mass fraction(S) profile
at $t=1000s$.

For constructing non-linear parametric ROMs we follow the same steps as
mentioned in the 1D case, solving for
$\mu\in\\{540\mathrm{K},550\mathrm{K},560\mathrm{K},570\mathrm{K},580\mathrm{K}\\}$
and sampling every $10^{th}$ snapshot in time to construct a parametric
snapshot matrix $Q$. With $N_{t}=1000/10=100$ time instances, $N_{p}=5$
parameter instances of $\mu$ and $M=N_{x}*N_{y}=250000$ the constructed matrix
$Q\in\mathbb{R}^{M\times N_{t}N_{p}}=\mathbb{R}^{250000\times 500}$. We then
carry out the dimensionality reduction through sPOD. We, however, need the
shifts prior to running the sPOD algorithm. To this end, we observe that the
2D problem is set up in a way that the developed flame front is perfectly
circular throughout the domain. Because of this radial symmetry, we can
transform the data from a Cartesian to a polar coordinate system which
simplifies the problem significantly. For performing model reduction in polar
coordinate system we have
$\underline{\bm{\Delta}}^{k}=(\underline{\Delta}^{k}_{\mathrm{R}},\underline{\Delta}^{k}_{\theta})$
as shifts.

Let us consider the Figure 7, the first plot shows the 1D cross-sectional
views of the temperature profile at two different time instances. $t=1000s$ is
the final instance that can be considered as the reference position. The
$\Delta$ is given as:

$\Delta=\underbrace{x(t=300s)}_{\mathrm{R}}-x(t=1000s)$

There is no dependency of the shifts on the $\theta$ coordinate. The $\Delta$
only depends upon $\mathrm{R}$ and $t$. The other two plots show the
temperature profiles in the polar coordinate system for the two-time instances
$t=300s$ and $t=1000s$.

Figure 7: The left plot shows the cross-sectional view of the temperature
profile for two different time instances. The other two plots on the right
show the temperature profiles at the mentioned time instances in the polar
coordinate system.

With sPOD, we now aim at decomposing the polar temperature field into two
frames: the first one captures the traveling fronts, and the second one
captures the initial ignition regime. For this we consider,
$\underline{\Delta}^{1}_{\mathrm{R}}=\Delta,\>\underline{\Delta}^{1}_{\theta}=0$
and $\underline{\bm{\Delta}}^{2}=(0,0)$. The sPOD decomposition of the
temperature is shown in Figure 8 in Cartesian coordinates.

Figure 8: This figure shows the sPOD decomposition of the temperature field at $t=100s$. The first plot shows the sPOD reconstruction whereas we observe the decomposition of the field into frame 1 consisting of the traveling wave part and frame 2 containing the initial ignition.  Table 4: Offline and online error study w.r.t number of modes for 2D wildland fire model(without wind) Modes111We also have $r=r_{1}+r_{2}+1$, where $r_{1}$ is the rank of the moving frame, $r_{2}$ is the rank of the stationary frame, and $1$ accounts for one additional degree of freedom controlling the shift. For comparison, we therefore use $r$ modes for the POD. | Offline errors | Online errors
---|---|---
$r_{1}+r_{2}$ | $E^{\mathrm{sPOD}}$ | $E^{\mathrm{POD}}$ | $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}$ | $E^{\mathrm{sPOD-I}}_{\mathrm{tot}}$ | $E^{\mathrm{POD-NN}}_{\mathrm{tot}}$
1 + 1 | $3.67\times 10^{-1}$ | $6.20\times 10^{-1}$ | $2.02\times 10^{-1}$ | $2.05\times 10^{-1}$ | $6.17\times 10^{-1}$
5 + 2 | $2.18\times 10^{-2}$ | $3.51\times 10^{-1}$ | $3.20\times 10^{-2}$ | $3.23\times 10^{-2}$ | $3.46\times 10^{-1}$
7 + 2 | $1.13\times 10^{-2}$ | $2.89\times 10^{-1}$ | $3.22\times 10^{-2}$ | $2.99\times 10^{-2}$ | $2.84\times 10^{-1}$
10 + 3 | $5.25\times 10^{-3}$ | $2.03\times 10^{-1}$ | $3.79\times 10^{-2}$ | $2.89\times 10^{-2}$ | $2.01\times 10^{-1}$
16 + 4 | $1.88\times 10^{-3}$ | $1.18\times 10^{-1}$ | $2.35\times 10^{-2}$ | $2.85\times 10^{-2}$ | $1.25\times 10^{-1}$
19 + 5 | $1.33\times 10^{-3}$ | $8.91\times 10^{-2}$ | $1.77\times 10^{-2}$ | $2.85\times 10^{-2}$ | $9.85\times 10^{-2}$
|  |  |  |  |

The basis reconstruction error results are shown in Table 4. We observe that
as we keep on adding more modes for the reconstruction the drop in
$E^{\mathrm{sPOD}}$ is significant in comparison to $E^{\mathrm{POD}}$. For
eg. sPOD converges for $19+5=24$ modes to
$E^{\mathrm{sPOD}}\sim\mathcal{O}(10^{-3})$ whereas
$E^{\mathrm{POD}}\sim\mathcal{O}(10^{-2})$.

As a result of the sPOD, we obtain the time amplitudes for both the frames and
we subsequently construct the time amplitude matrix $\hat{A}$ and the
parameter matrix $P$ for training as shown in Section 2. For this problem
$\hat{A}\in\mathbb{R}^{25\times 500}$ and $P\in\mathbb{R}^{2\times 500}$. The
network is trained for $N_{\mathrm{epochs}}=200000$ with batch size
$N_{b}=50$. Early stopping criteria are imposed to prevent overfitting, which
stops the training if the validation loss does not decrease for 4000
consecutive epochs. For testing, we choose $\mu=[558.49\mathrm{K}]$ to predict
the time amplitudes and the shifts which are shown in Figure 9. Quantitative
error estimates are shown in Figure 10 where in the first plot we observe the
similar stagnation behavior that was observed in the 1D wildland fire model
case. The $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}\sim 0.017$. In the second plot,
we see the error over the whole time interval $E^{\mathrm{sPOD-NN}}_{j}$ where
we see $\max<\mathcal{O}(10^{-2})$ and the mean $<\mathcal{O}(10^{-4})$.

Figure 9: Plot (a) and (b) show the time amplitude predictions of the first
mode for frame 1 and frame 2 respectively. (c) and (d) show the prediction
results for the first two modes corresponding to the POD-NN approach. (e) show
the shift prediction corresponding to the sPOD-NN approach. Figure 10: The
first plot shows the error decay for an increasing $n_{\mathrm{dof}}$. The
errors labeled $E_{\mathrm{tot}}$ show the final reconstruction error after
prediction or interpolation. The second plot shows the trend of the relative
error $E^{\mathrm{sPOD-NN}}_{j}$ over time.

Full reconstruction results are shown in Figure 11. Here we see the cross-
sectional view of the temperature profile at two different time instances. We
observe that the POD-NN has spurious oscillations. The $E^{\mathrm{POD}}\sim
0.089$ and $E^{\mathrm{POD-NN}}_{\mathrm{tot}}\sim 0.098$ thus we infer that
the prediction accuracy for POD-NN is once again limited by the basis
reconstruction error itself. As for sPOD-NN and sPOD-I, we see very minor
oscillations near the sharp edges of the curve.

Figure 11: The top two plots show the temperature profiles at $t=100s$ and
$t=1000s$. The bottom two plots show the cross-sectional view of their top
counterparts with the prediction results for sPOD-NN and POD-NN procedure.

####  3.1.3 2D model (with wind)

We consider the same model equations as in the previous case albeit with a
minor modification in Equation 43 where we prescribe $v=(0.2\>\mathrm{m/s},0)$
which is the wind velocity applied only in the $x$ direction. The two-
dimensional square domain is kept the same as the previous case and the
spatial domain is split into $N_{x}=500,N_{y}=500$ grid points making
$h_{x}=1,h_{y}=1$ and discretized with finite difference method with $6^{th}$
order central finite difference stencil. We solve the system for $500$ time
steps for $t_{f}=500s$. The temperature and the supply mass fraction profiles
are shown in Figure 12.

Figure 12: Full order model: Temperature(T) and supply mass fraction(S)
profile at $t=500s$.

Non-linear parametric ROM is constructed by solving for
$\mu\in\\{540\mathrm{K},550\mathrm{K},560\mathrm{K},570\mathrm{K},580\mathrm{K}\\}$
and sampling every $5^{th}$ snapshot in time to construct a parametric
snapshot matrix $Q\in\mathbb{R}^{250000\times 500}$. The sPOD is subsequently
applied for dimensionality reduction. On similar lines to the previous case,
we convert our data to a polar coordinate system for simplicity. However, this
time due to the changing topology of the flame front the calculation of shifts
is no longer trivial. Similar to the previous case we have
$\underline{\bm{\Delta}}^{k}=(\underline{\Delta}^{k}_{\mathrm{R}},\underline{\Delta}^{k}_{\theta})$
and $\underline{\bm{\Delta}}^{2}=(0,0)$. However,
$\underline{\Delta}^{1}_{\mathrm{R}}$ is now dependent on
$(t,\theta,\mathrm{R})$ for which we look at Figure 13 for detail.

Figure 13: The plots show the supply mass fractions at two different time
instances in the polar coordinate system.

We see the edges of the fronts in the left plot for two different time
instances from which we calculate the shifts for different values of $\theta$
as shown. With sPOD, we now aim at decomposing the polar temperature field
into two frames: the first one captures the traveling fronts, and the second
one captures the initial ignition regime. Subsequently we have
$\underline{\Delta}^{1}_{\mathrm{R}}\in\mathbb{R}^{250000\times
N},\>\underline{\Delta}^{1}_{\theta}=0$ where in
$\underline{\Delta}^{1}_{\mathrm{R}}$ each column gives the amount of shift to
be applied to all the grid points of the domain at instances of $N$. The sPOD
decomposition of the temperature is shown in Figure 14 in Cartesian
coordinates.

Figure 14: This figure shows the sPOD decomposition of the temperature field at $t=50s$. The first plot shows the sPOD reconstruction and the respective decomposition into the traveling part and the stationary initial ignition part are shown in the second and the third plot. Table 5: Offline and online error study w.r.t number of modes for 2D wildland fire model (with wind) Modes111We also have $r=r_{1}+r_{2}+4$, where $r_{1}$ is the rank of the moving frame, $r_{2}$ is the rank of the stationary frame, and $4$ accounts for the added number of amplitudes extracted from the shift matrix $\underline{\bm{\Delta}}^{1}_{\mathrm{R}}$. For comparison, we therefore use $r$ modes for the POD. | Offline errors | Online errors
---|---|---
$r_{1}+r_{2}$ | $E^{\mathrm{sPOD}}$ | $E^{\mathrm{POD}}$ | $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}$ | $E^{\mathrm{sPOD-I}}_{\mathrm{tot}}$ | $E^{\mathrm{POD-NN}}_{\mathrm{tot}}$
1 + 1 | $2.63\times 10^{-1}$ | $4.01\times 10^{-1}$ | $2.04\times 10^{-1}$ | $1.99\times 10^{-1}$ | $3.99\times 10^{-1}$
4 + 2 | $9.88\times 10^{-2}$ | $2.51\times 10^{-1}$ | $1.01\times 10^{-1}$ | $1.00\times 10^{-1}$ | $2.49\times 10^{-1}$
8 + 3 | $4.73\times 10^{-2}$ | $1.55\times 10^{-1}$ | $5.10\times 10^{-2}$ | $4.83\times 10^{-2}$ | $1.54\times 10^{-1}$
13 + 6 | $2.46\times 10^{-2}$ | $8.45\times 10^{-2}$ | $3.39\times 10^{-2}$ | $3.04\times 10^{-2}$ | $8.19\times 10^{-2}$
17 + 9 | $1.97\times 10^{-2}$ | $5.75\times 10^{-2}$ | $2.81\times 10^{-2}$ | $2.73\times 10^{-2}$ | $7.16\times 10^{-2}$
|  |  |  |  |

The basis reconstruction error results are shown in Table 5. We could reach to
$E^{\mathrm{sPOD}}\sim 0.01$ with $17+9=26$ modes whereas
$E^{\mathrm{POD}}\sim 0.06$.

sPOD gives us the time amplitudes for both the frames and we proceed to
construct the time amplitude matrix $\hat{A}$ for training as explained in
Section 2. The shift truncation dimension $n_{\Delta}^{1}=4$ is chosen
suitably keeping in mind the relative reconstruction error of the shift matrix
shown in Figure 15. This is one particular case of extracting low-dimensional
structure from the high-dimensional shifts which is desired to proceed with
the neural network training. The theory behind this is explained in detail in
Section 2.

Figure 15: Relative reconstruction error with respect to the number of shift
truncation modes

For this problem $\hat{A}\in\mathbb{R}^{30\times 500}$ and
$P\in\mathbb{R}^{2\times 500}$. The network parameters remain the same as the
previous case with $N_{\mathrm{epochs}}=200000$ and batch size $N_{b}=50$.
Early stopping criteria are imposed to prevent overfitting after 4000
consecutive epochs. For the test parameter, $\mu=[558.49\mathrm{K}]$ the time
amplitude and the shift predictions are in line with the previous case. We
also show the results of a parameter sweep study for the temperature for 16
different test samples where we sample the values for $\mu\in[540,580]$. The
result is shown in Figure 18. Error estimates are shown in Figure 16. In the
first plot we observe the $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}\sim 0.028$. In
the second plot, we see the error over the whole time interval
$E^{\mathrm{sPOD-NN}}_{j}$ where we see $\max<\mathcal{O}(10^{-2})$.

Figure 16: The first plot shows the error decay for an increasing
$n_{\mathrm{dof}}$. The second plot shows the trend of the relative error
$E^{\mathrm{sPOD-NN}}_{j}$ over time.

Full reconstruction results are shown in Figure 17. The cross-sectional views
of the temperature profile at two different time instances are shown. We
observe that the POD-NN is corrupted by oscillations but for sPOD-NN and
sPOD-I we do not see any oscillations.

Figure 17: The top two plots show the temperature profiles at $t=50s$ and
$t=500s$. The bottom two plots show the cross-sectional views of the profile
with the prediction results for sPOD-NN, POD-NN, and sPOD-I procedures. Figure
18: Parameter sweep study for 1D and 2D temperature values

###  3.2 Timing study

Here we show the computational time analysis of the proposed methods. All the
tests for computing the timing are run on Macbook Air M1(2020) with an 8-core
CPU and 16GB of RAM. We refer to Figure 19 where we see that for both 1D and
2D models, both sPOD-NN and sPOD-I methods are more than 100 times faster than
the FOM for the multi-query scenario. However, this speedup becomes very large
(more than 1000 times) for POD-NN. The reason for this could be explained by
looking at the sub-steps of sPOD-NN and POD-NN. For the 1D model, the sPOD-NN
timing consists of the time for evaluation of the neural network followed by
that of transforming the stationary frames into co-moving frames which
involves matrix-matrix multiplications for every frame and at last the time
for adding the individual frames and producing the final snapshot. However,
for POD-NN the timing only consists of evaluating the neural network and
reconstructing the final snapshot with just one matrix-matrix multiplication
operation. This difference translates into the 2D model as well except we have
one more sub-step for sPOD-NN where the snapshot data has to be converted from
a Cartesian coordinate system to a polar one for analysis and back again. This
makes the difference in speedups more pronounced. The speedup shown in Figure
19 is defined as:

$\text{Speedup}=\frac{t_{\mathrm{FOM}}}{t_{\mathrm{ROM}}}$ (48)

where $t_{\mathrm{ROM}}$ is the time taken for the methods proposed in the
paper. The $t_{\mathrm{FOM}}=5.2s$ for 1D model and $63.2s$ for the 2D model.
Both the test cases for the 2D model described in the paper have similar
computational time consumption. The comparison is made with the same
$n_{\mathrm{dof}}$ across all methods for a specific model. With the results
shown, we substantiate our claims that not only is the proposed method
accurate but is also extremely fast for the online phase.

Figure 19: Computational time analysis for the online phase. For 1D model
$n_{\mathrm{dof}}=44$, for 2D model (without wind) $n_{\mathrm{dof}}=25$ and
for 2D model (with wind) $n_{\mathrm{dof}}=30$. The $t_{\mathrm{FOM}}=5.2s$
for 1D model and $t_{\mathrm{FOM}}=63.2s$ for the 2D model.

Just for the sake of completeness, the time study for the offline phase is
shown in Table 6.

Table 6: Offline phase timings | Basis reconstruction time | Network training time
---|---|---
Model | sPOD | POD | sPOD-NN | POD-NN
1D |  $\sim 1.5$ h |  $54$ s |  $\sim 32$ min |  $\sim 32$ min
2D (without wind) |  $\sim 33$ h |  $34$ s |  $\sim 9$ min |  $\sim 9$ min
2D (with wind) |  $\sim 12$ min |  $34$ s |  $\sim 9$ min |  $\sim 9$ min
|  |  |  |

##  4 Conclusion

In this paper, we proposed a non-intrusive model reduction technique sPOD-NN
for parametric transport-dominated systems. In particular, the technique uses
sPOD for constructing the reduced basis and then subsequently extracts the
time amplitudes for training a neural network for the offline phase. In the
online phase, the trained model is then used to predict the time amplitudes at
unseen parameter values. The core concept builds upon the ability of sPOD,
over the conventional methods like POD in constructing an optimal reduced
basis for transport-dominated problems. For benchmarking we tested sPOD-NN
against two other methods: sPOD-I and POD-NN. We assessed the computational
performance and the prediction accuracy of the methods on non-linear time-
dependent parametrized PDE systems: 1D and 2D wildland fire models. Through
the numerical results shown in Section 3, we saw that sPOD-NN yielded accurate
numerical approximations for the time amplitudes and the shifts and in turn
the final snapshot. Along with the accurate predictions the proposed method
provides substantial speedups in the online phase compared to the FOM. We also
saw that the sPOD-NN outperforms POD-NN in all the examples presented owing to
the ability of sPOD to better construct the optimal reduced basis which will
almost always happen in transport-dominated systems. Also, based on the
results we saw that the scope of improvement for sPOD-NN is huge compared to
POD-NN as the total error is dominated by the neural network prediction error
for sPOD-NN. Given more number of training samples the network prediction
error could be further reduced.

For the scope of this paper, we considered models with either no wind presence
or a constant unidirectional wind. Although for such test cases, the presented
methods would suffice, they will have to be modified to handle more realistic,
non-uniform, and shape-changing fronts krah_front_2022 . Also, the wildland
fire model considered here although is enough to capture the overall dynamics
of the process broadly, more complex and detailed models could be considered
as test cases for the proposed methods. As a future step, our aim is to study
a more challenging test case where the front profiles change in a more
complicated way as time progresses. In such cases, it becomes even more
difficult to accurately compute the shifts from the snapshot data making it a
perfect candidate to run our proposed methods on and substantiate their
generality.

## Code and Data Availability

In order to facilitate reproducibility and transparency of the research
presented in this paper, the source code used for the experiments and analyses
is made publicly available. Interested parties can access the code at the
following GitHub repository:

https://github.com/MOR-transport/sPOD-NN-paper

We encourage researchers to make use of the code and to extend it for their
own research purposes. In addition, the trained network-data is available upon
request, and interested parties can contact the corresponding author for
access.

## Author Contribution Statement (CRediT)

In the following, we declare the authors’ contributions to this work.

Shubhaditya Burela: |  methodology, implementation of neural networks, computations, and visualizations, writing original draft, reviewing & editing
---|---
Philipp Krah: |  initial concept, methodology, implementation of sPOD algorithm, reviewing & editing
Julius Reiss: |  supervision, reviewing & editing

## References

  *   * (1) DaCosta, G., Derigo, D., Liberta, G., Durrant, T., San-miguel-ayanz, J.: European wildfire danger and vulnerability in a changing climate: towards integrating risk dimensions (2020). https://doi.org/10.2760/46951
  * (2) Kondylatos, S., Prapas, I., Ronco, M., Papoutsis, I., Camps‐Valls, G., Piles, M., Fernández‐Torres, M., Carvalhais, N.: Wildfire danger prediction and understanding with deep learning. Geophysical Research Letters 49(17) (2022). https://doi.org/10.1029/2022GL099368
  * (3) Sayad, Y.O., Mousannif, H., Al Moatassime, H.: Predictive modeling of wildfires: A new dataset and machine learning approach. Fire Safety Journal 104, 130–146 (2019). https://doi.org/10.1016/j.firesaf.2019.01.006
  * (4) Vilar, L., Herrera, S., Tafur-Garcia, E., Yebra, M., Martinez-Vega, J., Echavarria, P., Martin, M.P.: Modelling wildfire occurrence at regional scale from land use/cover and climate change scenarios. Environmental Modelling & Software 145, 105200 (2021). https://doi.org/10.1016/j.envsoft.2021.105200
  * (5) Valero, M.M., Jofre, L., Torres, R.: Multifidelity prediction in wildfire spread simulation: Modeling, uncertainty quantification and sensitivity analysis. Environmental Modelling & Software 141, 105050 (2021). https://doi.org/10.1016/j.envsoft.2021.105050
  * (6) Mandel, J., Bennethum, L.S., Beezley, J.D., Coen, J.L., Douglas, C.C., Kim, M., Vodacek, A.: A wildland fire model with data assimilation. Mathematics and Computers in Simulation 79(3), 584–606 (2008). https://doi.org/10.1016/j.matcom.2008.03.015
  * (7) Benner, P., Ohlberger, M., Cohen, A., Willcox, K.: Model Reduction and Approximation. Computational Science & Engineering. Society for Industrial and Applied Mathematics, Philadelphia, PA (2017). https://doi.org/10.1137/1.9781611974829
  * (8) Hesthaven, J.S., Rozza, G., Stamm, B.: Certified Reduced Basis Methods for Parametrized Partial Differential Equations. SpringerBriefs in Mathematics. Springer, Cham (2016). https://doi.org/10.1007/978-3-319-22470-1
  * (9) Quarteroni, A., Manzoni, A., Negri, F.: Reduced Basis Methods for Partial Differential Equations vol. 92. Springer, Cham (2015). https://doi.org/10.1007/978-3-319-15431-2
  * (10) Greif, C., Urban, K.: Decay of the Kolmogorov n-width for wave problems. Applied Mathematics Letters 96, 216–222 (2019). https://doi.org/10.1016/j.aml.2019.05.013
  * (11) Unger, B., Gugercin, S.: Kolmogorov n-widths for linear dynamical systems. Advances in Computational Mathematics 45(5), 2273–2286 (2019). https://doi.org/10.1007/s10444-019-09701-0
  * (12) Krah, P., Büchholz, S., Häringer, M., Reiss, J.: Front Transport Reduction for Complex Moving Fronts. Accepted for publication in J. Sci. Comput. (2022). https://doi.org/10.48550/arXiv.2202.08208,
  * (13) Huang, C., Duraisamy, K., Merkle, C.: Challenges in Reduced Order Modeling of Reacting Flows. In: 2018 Joint Propulsion Conference. AIAA Propulsion and Energy Forum. American Institute of Aeronautics and Astronautics, Cincinnati, Ohio (2018). https://doi.org/10.2514/6.2018-4675
  * (14) Ohlberger, M., Rave, S.: Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing. Comptes Rendus Mathematique 351(23), 901–906 (2013). https://doi.org/10.1016/j.crma.2013.10.028
  * (15) Rim, D., Moe, S., LeVeque, R.J.: Transport Reversal for Model Reduction of Hyperbolic Partial Differential Equations. SIAM/ASA Journal on Uncertainty Quantification 6(1), 118–150 (2018). https://doi.org/10.1137/17M1113679
  * (16) Cagniart, N., Maday, Y., Stamm, B.: Model Order Reduction for Problems with Large Convection Effects. In: Contributions to Partial Differential Equations and Applications. Computational Methods in Applied Sciences, pp. 131–150. Springer, Cham (2019). https://doi.org/10.1007/978-3-319-78325-3_10
  * (17) Welper, G.: Transformed Snapshot Interpolation with High Resolution Transforms. SIAM Journal on Scientific Computing 42(4), 2037–2061 (2020). https://doi.org/10.1137/19M126356X
  * (18) Black, F., Schulze, P., Unger, B.: Projection-based model reduction with dynamically transformed modes. ESAIM: Mathematical Modelling and Numerical Analysis 54(6), 2011–2043 (2020). https://doi.org/10.1051/m2an/2020046
  * (19) Nonino, M., Ballarin, F., Rozza, G., Maday, Y.: Overcoming slowly decaying Kolmogorov n-width by transport maps: application to model order reduction of fluid dynamics and fluid–structure interaction problems. arXiv (2019). https://doi.org/10.48550/arXiv.1911.06598
  * (20) Krah, P., Sroka, M., Reiss, J.: Model Order Reduction of Combustion Processes with Complex Front Dynamics. In: Numerical Mathematics and Advanced Applications ENUMATH 2019\. Lecture Notes in Computational Science and Engineering, pp. 803–811. Springer, Cham (2021). https://doi.org/10.1007/978-3-030-55874-1_79
  * (21) Black, F., Schulze, P., Unger, B.: Efficient Wildland Fire Simulation via Nonlinear Model Order Reduction. Fluids 6(8), 280 (2021). https://doi.org/10.3390/fluids6080280
  * (22) Reiss, J., Schulze, P., Sesterhenn, J., Mehrmann, V.: The Shifted Proper Orthogonal Decomposition: A Mode Decomposition for Multiple Transport Phenomena. SIAM Journal on Scientific Computing 40(3), 1322–1344 (2018). https://doi.org/10.1137/17M1140571
  * (23) Lee, K., Carlberg, K.T.: Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. Journal of Computational Physics 404, 108973 (2020). https://doi.org/10.1016/j.jcp.2019.108973
  * (24) Reiss, J.: Optimization-Based Modal Decomposition for Systems with Multiple Transports. SIAM Journal on Scientific Computing 43(3), 2079–2101 (2021). https://doi.org/10.1137/20M1322005
  * (25) Krah, P.L.: Non-linear reduced order modeling for transport dominated fluid systems [doctoral dissertation] (2023). https://doi.org/10.14279/depositonce-16974
  * (26) Mendible, A., Brunton, S.L., Aravkin, A.Y., Lowrie, W., Kutz, J.N.: Dimensionality reduction and reduced-order modeling for traveling wave physics. Theoretical and Computational Fluid Dynamics 34, 385–400 (2020). https://doi.org/10.1007/s00162-020-00529-9
  * (27) Papapicco, D., Demo, N., Girfoglio, M., Stabile, G., Rozza, G.: The Neural Network shifted-proper orthogonal decomposition: A machine learning approach for non-linear reduction of hyperbolic equations. Computer Methods in Applied Mechanics and Engineering 392, 114687 (2022). https://doi.org/%****␣sn-article.bbl␣Line␣475␣****10.1016/j.cma.2022.114687
  * (28) Sarna, N., Grundel, S.: Hyper-reduction for parametrized transport dominated problems via online-adaptive reduced meshes. arXiv (2021). https://doi.org/10.48550/arXiv.2003.06362
  * (29) Carlberg, K., Farhat, C., Cortial, J., Amsallem, D.: The GNAT method for nonlinear model reduction: Effective implementation and application to computational fluid dynamics and turbulent flows. Journal of Computational Physics 242, 623–647 (2013). https://doi.org/10.1016/j.jcp.2013.02.028
  * (30) Barrault, M., Maday, Y., Nguyen, N.C., Patera, A.T.: An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations. Comptes Rendus Mathematique 339(9), 667–672 (2004). https://doi.org/10.1016/j.crma.2004.08.006
  * (31) Jain, S., Tiso, P.: Hyper-Reduction Over Nonlinear Manifolds for Large Nonlinear Mechanical Systems. Journal of Computational and Nonlinear Dynamics 14, 081008–081008 (2017). https://doi.org/10.1115/1.4043450
  * (32) Chaturantabut, S., Sorensen, D.C.: Nonlinear Model Reduction via Discrete Empirical Interpolation. SIAM Journal on Scientific Computing 32(5), 2737–2764 (2010). https://doi.org/10.1137/090766498
  * (33) Peherstorfer, B.: Model Reduction for Transport-Dominated Problems via Online Adaptive Bases and Adaptive Sampling. SIAM Journal on Scientific Computing 42(5), 2803–2836 (2020). https://doi.org/10.1137/19M1257275
  * (34) Koch, O., Lubich, C.: Dynamical Low‐Rank Approximation. SIAM Journal on Matrix Analysis and Applications 29(2), 434–454 (2007). https://doi.org/10.1137/050639703. Publisher: Society for Industrial and Applied Mathematics
  * (35) Etter, P.A., Carlberg, K.T.: Online adaptive basis refinement and compression for reduced-order models via vector-space sieving. Computer Methods in Applied Mechanics and Engineering 364, 112931 (2020)
  * (36) Dihlmann, M.A., Drohmann, M., Haasdonk, B.: Model reduction of parametrized evolution problems using the reduced basis method with adaptive time partitioning. Proc. of ADMOS, 2011:64 (2011)
  * (37) Zimmermann, R., Peherstorfer, B., Willcox, K.: Geometric Subspace Updates with Applications to Online Adaptive Nonlinear Model Reduction. SIAM Journal on Matrix Analysis and Applications 39(1), 234–261 (2018). https://doi.org/%****␣sn-article.bbl␣Line␣625␣****10.1137/17M1123286
  * (38) Peherstorfer, B., Willcox, K.: Online Adaptive Model Reduction for Nonlinear Systems via Low-Rank Updates. SIAM Journal on Scientific Computing 37(4), 2123–2150 (2015). https://doi.org/10.1137/140989169
  * (39) Koellermeier, J., Krah, P., Kusch, J.: Macro-micro decomposition for consistent and conservative model order reduction of hyperbolic shallow water moment equations: A study using POD-Galerkin and dynamical low rank approximation. arXiv (2023). https://doi.org/10.48550/arXiv.2302.01391
  * (40) Gonzalez, F.J., Balajewicz, M.: Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems. arXiv (2018). https://doi.org/10.48550/arXiv.1808.01346
  * (41) Kim, Y., Choi, Y., Widemann, D., Zohdi, T.: A fast and accurate physics-informed neural network reduced order model with shallow masked autoencoder. Journal of Computational Physics 451, 110841 (2022). https://doi.org/10.1016/j.jcp.2021.110841
  * (42) Salvador, M., Dedè, L., Manzoni, A.: Non intrusive reduced order modeling of parametrized PDEs by kernel POD and neural networks. Computers & Mathematics with Applications 104, 1–13 (2021). https://doi.org/10.1016/j.camwa.2021.11.001
  * (43) San, O., Maulik, R.: Neural network closures for nonlinear model order reduction. Advances in Computational Mathematics 44 (2018). https://doi.org/10.1007/s10444-018-9590-z
  * (44) Hesthaven, J.S., Ubbiali, S.: Non-intrusive reduced order modeling of nonlinear problems using neural networks. Journal of Computational Physics 363, 55–78 (2018). https://doi.org/10.1016/j.jcp.2018.02.037
  * (45) Wang, Q., Hesthaven, J.S., Ray, D.: Non-intrusive reduced order modeling of unsteady flows using artificial neural networks with application to a combustion problem. Journal of Computational Physics 384, 289–307 (2019). https://doi.org/10.1016/j.jcp.2019.01.031
  * (46) Bērziņš, A., Helmig, J., Key, F., Elgeti, S.: Standardized Non-Intrusive Reduced Order Modeling Using Different Regression Models With Application to Complex Flow Problems. arXiv (2021). https://doi.org/10.48550/arXiv.2006.13706
  * (47) Kast, M., Guo, M., Hesthaven, J.S.: A non-intrusive multifidelity method for the reduced order modeling of nonlinear problems. Computer Methods in Applied Mechanics and Engineering 364, 112947 (2020). https://doi.org/10.1016/j.cma.2020.112947
  * (48) Guo, M., Hesthaven, J.S.: Data-driven reduced order modeling for time-dependent problems. Computer Methods in Applied Mechanics and Engineering 345, 75–99 (2019). https://doi.org/10.1016/j.cma.2018.10.029
  * (49) Bhattacharya, K., Hosseini, B., Kovachki, N.B., Stuart, A.M.: Model Reduction And Neural Networks For Parametric PDEs. The SMAI Journal of computational mathematics 7, 121–157 (2021). https://doi.org/10.5802/smai-jcm.74
  * (50) Mohan, A.T., Gaitonde, D.V.: A Deep Learning based Approach to Reduced Order Modeling for Turbulent Flow Control using LSTM Neural Networks. arXiv (2018). https://doi.org/10.48550/arXiv.1804.09269
  * (51) J. Nagoor Kani, Elsheikh, A.H.: Reduced-Order Modeling of Subsurface Multi-phase Flow Models Using Deep Residual Recurrent Neural Networks. Transport in Porous Media 126(3), 713–741 (2019). https://doi.org/10.1007/s11242-018-1170-7
  * (52) Fresca, S., Dede’, L., Manzoni, A.: A Comprehensive Deep Learning-Based Approach to Reduced Order Modeling of Nonlinear Time-Dependent Parametrized PDEs. Journal of Scientific Computing 87(2), 61 (2021). https://doi.org/10.1007/s10915-021-01462-7
  * (53) Fresca, S., Manzoni, A.: POD-DL-ROM: Enhancing deep learning-based reduced order models for nonlinear parametrized PDEs by proper orthogonal decomposition. Computer Methods in Applied Mechanics and Engineering 388, 114181 (2022). https://doi.org/10.1016/j.cma.2021.114181
  * (54) Kovárnová, A., Krah, P., Reiss, J., Isoz, M.: Shifted Proper Orthogonal Decomposition and Artificial Neural Networks for Time-Continuous Reduced Order Models of Transport-Dominated Systems (2022). https://doi.org/10.14311/TPFM.2022.016
  * (55) Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, ., Feng, Y., Moore, E.W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P.: SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods 17(3), 261–272 (2020). https://doi.org/10.1038/s41592-019-0686-2
  * (56) O’Rourke, J.: Computational Geometry In C, 2nd edn. Cambridge University Press, Cambridge (1998). https://doi.org/10.1017/CBO9780511804120
  * (57) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., Chintala, S.: PyTorch: An Imperative Style, High-Performance Deep Learning Library. In: Advances in Neural Information Processing Systems, vol. 32. Curran Associates, Inc., Vancouver, Canada (2019)
  * (58) Karcher, N., Franz, T.: Adaptive sampling strategies for reduced-order modeling. CEAS Aeronautical Journal 13(2), 487–502 (2022). https://doi.org/10.1007/s13272-022-00574-6
  * (59) Jagtap, A.D., Karniadakis, G.E.: How important are activation functions in regression and classification? A survey, performance comparison, and future directions. arXiv (2022). https://doi.org/10.48550/arXiv.2209.02681
  * (60) Nwankpa, C., Ijomah, W., Gachagan, A., Marshall, S.: Activation Functions: Comparison of trends in Practice and Research for Deep Learning. arXiv (2018). https://doi.org/10.48550/arXiv.1811.03378

## Appendix

##  Appendix A Neural network foundation

In recent years, non-intrusive model reduction methods that utilize neural
networks have gained considerable attention. Despite advancements in the type
and architecture of these networks, the core concept remains unchanged. In
this section, we describe the specific neural network architecture used in the
numerical examples presented in this work.

For the architecture, we use a deep FNN (Figure 20). For a clearer picture,
consider a feed-forward network with $L$ hidden layers where
$l\in\\{1,\ldots,L\\}$ is the index of the hidden layer, for such a setup the
forward propagation is described as:

$\displaystyle\begin{split}z^{(l+1)}_{i}&=\bm{w}^{(l+1)}_{i}\bm{y}^{(l)}+b^{(l+1)}_{i},\\\
y^{(l+1)}_{i}&=f(z^{(l+1)}_{i})\end{split}$ (49)

where $\bm{y}^{(l)}$ is the vector of outputs from the layer $l$,
$\bm{z}^{(l)}$ is the vector of inputs into the layer $l$, $W^{(l)}$ is the
weight matrix in the layer $l$ along with a bias vector for the layer $l$,
$\bm{b}^{(l)}$. The activation function is given by $f(\cdot)$. In the
backpropagation step, the training of the network is usually carried out with
a gradient descent algorithm minimizing a certain loss function $J(W,\bm{b})$.
As we deal with a regression problem in our work we basically rely on two
types of loss functions, namely MAE (Mean Absolute Error) often referred to as
L1 loss, and MSE (Mean Squared Error):

$MAE=\frac{1}{n}\sum^{n}_{i=1}\left\lvert y_{i}-\hat{y}_{i}\right\rvert,\qquad
MSE=\frac{1}{n}\sum^{n}_{i=1}(y_{i}-\hat{y}_{i})^{2}$ (50)

where $y_{i}$ and $\hat{y}_{i}$ are the input and the target vectors
respectively. The training data in our problems is heavily skewed and contains
critical outliers that impact the physical significance of the problem. Using
Mean Squared Error (MSE) as a loss function can penalize these outliers and
produce subpar prediction results. To address this, we use Mean Absolute Error
(MAE) as the preferred loss function due to its robustness in the presence of
outliers. A stochastic gradient descent algorithm with momentum term is
utilized for minimizing the loss function and updating the parameters $W$ and
$\bm{b}$ in every iteration as shown:

$\displaystyle\begin{split}&W^{(l)}_{ij}\leftarrow
W^{(l)}_{ij}-\alpha\frac{\partial}{\partial W^{(l)}_{ij}}J(W,\bm{b}),\\\
&b^{(l)}_{i}\leftarrow b^{(l)}_{i}-\alpha\frac{\partial}{\partial
b^{(l)}_{i}}J(W,\bm{b})\end{split}$ (51)

$\alpha$ is the learning rate and can be tuned to get a better decay of the
loss function. $W^{(l)}_{ij}$ is an element of the weight matrix associated
with the connection between a node $j$ in layer $l$ and a node $i$ in layer
$l+1$ and the term $b^{(l)}_{i}$ is an element of the bias vector associated
with the node $i$ in layer $l+1$.

Figure 20: Feed-forward neural network

The Table 1 outlines the activation functions used. The first three layers use
Exponential Linear Unit (ELU) activation, and the next layer uses Leaky
Rectified Linear Unit (LeakyReLU). It is common in the literature to not use
any activation function for the output layer in regression problems, which we
also follow in this case. The detailed comparison of the different activation
functions can be found in jagtap_how_2022 , nwankpa_activation_2018 .

##  Appendix B Results for supply mass fraction

Here we present the results for the supply mass fraction for all the models
presented in the text. Table 7, Table 8 and Table 9 show the offline and
online errors for 1D, 2D(without wind) and 2D(with wind) models. Figure 21
shows the parameter sweep study for supply mass fraction where we test the
model on the samples of $\mu\in[540,580]$

Table 7: Offline and online error study w.r.t number of modes for 1D wildland fire model (supply mass fraction) | Offline errors | Online errors
---|---|---
Modes | $E^{\mathrm{sPOD}}$ | $E^{\mathrm{POD}}$ | $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}$ | $E^{\mathrm{sPOD-I}}_{\mathrm{tot}}$ | $E^{\mathrm{POD-NN}}_{\mathrm{tot}}$
1 + 1 + 1 | $8.95\times 10^{-2}$ | $1.53\times 10^{-1}$ | $6.27\times 10^{-2}$ | $6.17\times 10^{-2}$ | $1.53\times 10^{-1}$
4 + 2 + 4 | $2.15\times 10^{-2}$ | $6.76\times 10^{-2}$ | $2.19\times 10^{-2}$ | $1.79\times 10^{-2}$ | $6.70\times 10^{-2}$
8 + 2 + 8 | $3.88\times 10^{-3}$ | $4.05\times 10^{-2}$ | $8.91\times 10^{-3}$ | $8.62\times 10^{-3}$ | $4.06\times 10^{-2}$
12 + 4 + 12 | $9.00\times 10^{-4}$ | $2.44\times 10^{-2}$ | $7.35\times 10^{-3}$ | $7.81\times 10^{-3}$ | $2.78\times 10^{-2}$
|  |  |  |  |
Table 8: Offline and online error study w.r.t number of modes for 2D (without wind) wildland fire model (supply mass fraction) | Offline errors | Online errors
---|---|---
Modes | $E^{\mathrm{sPOD}}$ | $E^{\mathrm{POD}}$ | $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}$ | $E^{\mathrm{sPOD-I}}_{\mathrm{tot}}$ | $E^{\mathrm{POD-NN}}_{\mathrm{tot}}$
1 + 1 | $9.60\times 10^{-2}$ | $1.24\times 10^{-1}$ | $7.71\times 10^{-2}$ | $7.79\times 10^{-2}$ | $1.24\times 10^{-1}$
4 + 1 | $4.20\times 10^{-2}$ | $7.59\times 10^{-2}$ | $4.17\times 10^{-2}$ | $4.14\times 10^{-2}$ | $7.51\times 10^{-2}$
7 + 3 | $7.45\times 10^{-3}$ | $4.57\times 10^{-2}$ | $1.04\times 10^{-2}$ | $9.06\times 10^{-3}$ | $4.57\times 10^{-2}$
11 + 4 | $5.39\times 10^{-3}$ | $3.10\times 10^{-2}$ | $8.73\times 10^{-3}$ | $7.77\times 10^{-3}$ | $3.15\times 10^{-2}$
|  |  |  |  |
Table 9: Offline and online error study w.r.t number of modes for 2D (with wind) wildland fire model (supply mass fraction) | Offline errors | Online errors
---|---|---
Modes | $E^{\mathrm{sPOD}}$ | $E^{\mathrm{POD}}$ | $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}$ | $E^{\mathrm{sPOD-I}}_{\mathrm{tot}}$ | $E^{\mathrm{POD-NN}}_{\mathrm{tot}}$
1 + 1 | $7.43\times 10^{-2}$ | $3.67\times 10^{-2}$ | $7.15\times 10^{-2}$ | $7.13\times 10^{-2}$ | $3.53\times 10^{-2}$
2 + 2 | $2.95\times 10^{-2}$ | $2.92\times 10^{-2}$ | $2.09\times 10^{-2}$ | $2.08\times 10^{-2}$ | $2.73\times 10^{-2}$
4 + 3 | $1.71\times 10^{-2}$ | $2.23\times 10^{-2}$ | $1.82\times 10^{-2}$ | $1.82\times 10^{-2}$ | $2.17\times 10^{-2}$
|  |  |  |  |
Figure 21: Parameter sweep study for 1D and 2D supply mass fraction values

##  Appendix C sPOD parameters

As suggested by black_efficient_2021  we use the sPOD for decomposing the
snapshot data of the wildland fire model. To determine the co-moving fields
$\\{Q^{k}\\}_{k=1,\dots,f}$ shown in Equation 4 and an error term
$E\in\mathbb{R}^{M\times N}$ to capture noise, we solve the following
constraint optimization problem:

$\displaystyle\min_{Q^{k},E}\sum_{k=1}^{f}\left\lVert
Q^{k}\right\rVert_{*}+\varsigma\left\lVert E\right\rVert_{1}\quad\text{s.t.
}\quad Q=\sum_{k=1}^{f}T^{\Delta^{k}}(Q^{k})+E\,,$ (52)

using algorithm 8 in krah_non-linear_nodate . Here, $\|E\|_{1}=\sum_{ij}\lvert
E_{ij}\rvert$ is not the usual matrix 1-norm, but the vector 1-norm of a long
vector $E\in\mathbb{R}^{M\times N}$, $\|A\|_{*}=\sum_{i}\sigma_{i}(A)$ is the
Schatten one-norm and $\varsigma\geq 0,\eta>0$ are tuning parameters. For the
1D and 2D wildland fire model, we state $\varsigma,\eta$ in Table 10. It is
important to note that once the computation of ${Q^{k}}_{k=1,\dots,f}$ and
$E\in\mathbb{R}^{M\times N}$ is complete, $E\in\mathbb{R}^{M\times N}$ is
disregarded as it solely captures the data noise.

Table 10: sPOD parameters for 1D and 2D wildland fire models Model | $\eta$ | $\varsigma$ | $N_{\mathrm{iter}}$
---|---|---|---
1D | $5.47\times 10^{-6}$ | $1.85\times 10^{-1}$ | $2000$
2D (without wind) | $3.92\times 10^{-7}$ | $4.47\times 10^{-2}$ | $2000$
2D (with wind) | $1.18\times 10^{-5}$ | $4.47\times 10^{2}$ | $7$

##  Appendix D Synthetic test case

As a proof of concept for the proposed method, we consider a 1D test case of
traveling waves $q^{1}$ and $q^{2}$. We assume that both waves move along the
path:

$\underline{\bm{\Delta}}^{1}(t,\mu)=\mu
t,\quad\quad\text{and}\quad\quad\underline{\bm{\Delta}}^{2}(t,\mu)=-\mu t$
(53)

The path is parameterized by the parameter $\mu$. We have
$\mu\in(\mu_{1},\ldots,\mu_{N_{p}})\in\mathcal{P}\subset\mathbb{R}^{N_{p}}$
which means that for different values of the parameter $\mu$ the paths of the
traveling waves will change. We then assume a function $q(x,t,\mu)$ which is
constructed by the superposition of $q^{1}$ and $q^{2}$:

$q(x,t,\mu)=q^{1}(x+\underline{\bm{\Delta}}^{1}(t,\mu),t,\mu)+q^{2}(x+\underline{\bm{\Delta}}^{2}(t,\mu),t,\mu),$
(54)

where $(x,t,\mu)\in[-L/2,L/2[\times[0,T[\times\mathcal{P}$, with $M$ number of
grid points in the spatial domain, $N_{t}$ time steps and $N_{p}$ instances of
the parameter $\mu$. To be as close as possible to a realistic setting, we
assume that $q(x,t,\mu)$ is the solution of a high dimensional ODE resulting
from a discretized PDE. We also assume that we only have access to the
snapshot matrix $Q_{i,j+N_{t}*p}=q(x_{i},t_{j},\mu_{p})$ and the paths
$\underline{\bm{\Delta}}^{k}(t_{j})$ for $N_{p}$ parameters and where $k$ is
the number of individual waves being superposed to construct the solution. The
individual waves for a single parameter instance $\mu_{p}$ are defined as:

$q^{k}(x,t,\mu_{p}):=\sum^{n_{s}-1}_{n=0}\mu_{p}\left(1+e^{-2nt}\right)\cos\left(-2\pi
t/T(n+1)\right)h_{n}(x^{k}/w),\\\ \langle h_{n},h_{m}\rangle=\delta_{n,m}$
(55)

where $w=0.015L$ and $x^{k}=x\pm 0.1L$ for $k=1,2$ respectively. Both $q^{1}$
and $q^{2}$ are defined in terms of dyadic pairs as shown in Equation 55 as it
enables us to tune the singular value spectra ($n_{s}$ being the number of
singular values) in each frame. Furthermore, we choose Gauss-Hermite
polynomials $h_{n}$

$h_{n}(x)=\frac{(-1)^{n}}{\sqrt{2^{n}n!\sqrt{\pi}}}e^{x^{2}/2}\frac{d^{n}}{dx^{n}}e^{-x^{2}}=\frac{1}{\sqrt{2^{n}n!\sqrt{\pi}}}H_{n}(x)e^{-\frac{1}{2}x^{2}}$
(56)

as they nicely mimic strongly localized wave structures (since $w\ll L$). To
build the data, we consider $M=500,\>N_{t}=500$ for $L=1$ and $T=1$, with
$\mu\in[0.1,0.15,0.2,0.25,0.3]$, $k=2$, $n_{s}=8$. The superposed solution $q$
and the individual stationary waves $q^{1}$ and $q^{2}$ for all the five
different $\mu$ are shown in the three rows of Figure 22 respectively. Note
that the pictures have the same structure
$Q_{i,j+N_{t}*p}=q(x_{i},t_{j},\mu_{p})$ as the snapshot matrices
$Q,Q^{1},Q^{2}$.

Figure 22: The second and the third row of the plot shows $q^{1}$ and $q^{2}$
respectively which are the individual waves and the top plot shows $q$ which
is the shifted superposed solution computed at $\mu=[0.1,0.15,0.2,0.25,0.3]$
respectively.

However, for further analysis, we assume that we have access to only $Q$ for
every $\mu$ so for $M=500$, $N_{t}=500$ and $N_{p}=5$ the collection of
snapshots results in a matrix of size $Q\in\mathbb{R}^{M\times
N}=\mathbb{R}^{500\times 2500}$. In the first step, we apply sPOD on this
matrix $Q$ and obtain a low-dimensional representation. We let the algorithm
run for optimal decomposition such that the
$E^{\mathrm{sPOD}}\sim\mathcal{O}(10^{-4})$. We use the ansatz
$Q\approx\tilde{Q}=T^{\bm{\Delta}^{1}}Q^{1}+T^{\bm{\Delta}^{2}}Q^{2}$. The
result of the decomposition is shown for $\mu_{0}$ in Figure 23.

Figure 23: sPOD decomposition on $Q$ corresponding to $\mu_{0}=0.1$. The first
two images are the decomposed stationary frames and the last image is the
reconstructed $Q$.

The low-dimensional description for each of the sPOD decomposed frames $Q^{k}$
for $\mu_{0}$ is given as:

$Q^{k}=q^{k}(x,t,\vec{\mu})\approx\sum^{r_{k}}_{i=1}a^{k}_{i}(t,\vec{\mu})\phi^{k}_{i}(x),\quad
k=1,2$ (57)

The time amplitude matrix $A^{k}$ can be extracted as shown in Equation 16 and
also we already have access to the shifts for all $k$. Thus for training the
neural network model, we assemble the time amplitude matrix $\hat{A}$ and the
parameter matrix $P$ as shown in Equation 20 and Equation 9. For our problem
$\hat{A}\in\mathbb{R}^{18\times 2500}$ and $P\in\mathbb{R}^{2\times 2500}$. We
run the training loop for $N_{\mathrm{epochs}}=150000$ with batch size
$N_{b}=50$. For testing, we choose $\mu=[0.23]$ and use the trained model to
predict the time amplitudes and the shifts which are shown in Figure 24.

Figure 24: (a) and (b) show the prediction results for the first mode of both
the frames with sPOD-NN and sPOD-I. (c) is the prediction result for the first
mode with the POD-NN approach. (d) and (e) show the prediction results for the
shifts for both frames.

We assess the prediction accuracy of the proposed methods through Figure 25.
In the first plot, we see that with an increasing number of modes the
$E^{\mathrm{sPOD}}$ decreases rapidly compared to $E^{\mathrm{POD}}$ which
substantiates our claim about using sPOD for basis construction. Here, it
needs to be pointed out that $E_{\mathrm{tot}}$ is bounded from below by the
basis reconstruction error behavior thus we see that $E^{\mathrm{POD-
NN}}_{\mathrm{tot}}$ is always higher than $E^{\mathrm{POD}}$. Same is true
for $E^{\mathrm{sPOD-NN}}_{\mathrm{tot}}$. It is also noted that the
prediction accuracy of the sPOD-NN is an order of magnitude better than POD-
NN. Surprisingly enough for this particular setting, the neural network is not
able to learn the POD modes at all. We see that with an increasing
$n_{\mathrm{dof}}$ the $E^{\mathrm{POD-NN}}_{\mathrm{tot}}$ remains very near
to 1. The performance for the POD-NN could improve if we add more parameter
samples to the training data. However, in this test case, the sPOD-I performs
the best. It predicts nearly at the limit of $E^{\mathrm{sPOD}}$, without
incurring any extra error in the interpolation step. This can be attributed to
the fact that the basis functions for the aforementioned problem are
polynomials by design and the sPOD-I uses polynomials for interpolation.

We can get more insight into the prediction accuracy of sPOD-NN by looking at
the second plot in Figure 25. It shows the maximum, minimum, and mean of the
relative error $E(t,\mu)$ for all the time instances. We see that the maximum
error recorded throughout the time frame is around $\sim 0.06$ however, the
mean error is always $\sim\mathcal{O}(10^{-3})$.

Figure 25: The first plot shows the error decay for an increasing
$n_{\mathrm{dof}}$. $E^{\mathrm{sPOD}}$, $E^{\mathrm{POD}}$ show the basis
reconstruction error for sPOD and POD respectively, and the errors labeled
$E_{\mathrm{tot}}$ show the final reconstruction error after prediction or
interpolation. The second plot shows the trend of the relative error
$E^{\mathrm{sPOD-NN}}_{j}$ over time.

We now look at the final reconstructed data in Figure 26. We observe that the
profile with $Q^{\mathrm{POD-NN}}$ is completely distorted, which makes sense
when we look at the errors shown in Figure 25. However, there is only a minor
visible difference in the profiles of sPOD-NN and sPOD-I compared to the
original.

Figure 26: Reconstructed snapshot comparison for all the methods