# Progressive reduced order modeling: empowering data-driven modeling with
selective knowledge transfer

Teeratorn Kadeethum  Sandia National Laboratories, Albuquerque, NM, 87185, USA
Daniel O’Malley  Los Alamos National Laboratory, Los Alamos, NM, 87545, USA
Youngsoo Choi  Lawrence Livermore National Laboratory, Livermore, CA, 94550,
USA  Hari S. Viswanathan  Los Alamos National Laboratory, Los Alamos, NM,
87545, USA  Hongkyu Yoon  Sandia National Laboratories, Albuquerque, NM,
87185, USA<EMAIL_ADDRESS>
###### Abstract

Data-driven modeling can suffer from a constant demand for data, leading to
reduced accuracy and impractical for engineering applications due to the high
cost and scarcity of information. To address this challenge, we propose a
progressive reduced order modeling framework that minimizes data cravings and
enhances data-driven modeling’s practicality. Our approach selectively
transfers knowledge from previously trained models through gates, similar to
how humans selectively use valuable knowledge while ignoring unuseful
information. By filtering relevant information from previous models, we can
create a surrogate model with minimal turnaround time and a smaller training
set that can still achieve high accuracy. We have tested our framework in
several cases, including transport in porous media, gravity-driven flow, and
finite deformation in hyperelastic materials. Our results illustrate that
retaining information from previous models and utilizing a valuable portion of
that knowledge can significantly improve the accuracy of the current model. We
have demonstrated the importance of progressive knowledge transfer and its
impact on model accuracy with reduced training samples. For instance, our
framework with four parent models outperforms the no-parent counterpart
trained on data nine times larger. Our research unlocks data-driven modeling’s
potential for practical engineering applications by mitigating the data
scarcity issue. Our proposed framework is a significant step toward more
efficient and cost-effective data-driven modeling, fostering advancements
across various fields.

## Introduction

In our quest to sustain global economics and address population growth, we
face the pressing need to produce energy sustainably, extract groundwater
responsibly, and store undesirable substances underground [1, 2]. However,
understanding the complex multiphysics phenomena that govern subsurface
physics and characterize underground structures poses significant challenges.
These structures are often anisotropic, heterogeneous, and exhibit
discontinuous material properties that span multiple orders of magnitude [3,
4]. Traditional approaches, such as lab experiments, are expensive, time-
consuming, and limited in their ability to simulate complex initial and
boundary conditions [5]. On the other hand, field trials can be risky,
compromising operational safety and potentially leading to irreversible damage
to assets. Furthermore, a lack of understanding of these complex behaviors may
prevent optimizing development plans [6]. High-fidelity numerical models like
finite difference, finite volume, or finite element approaches present an
alternative to overcome these limitations. These models can effectively
simulate coupled processes across different scales, accommodating complex
geometries, boundary conditions, and heterogeneous materials [7]. However,
their use often entails a high computational cost, mainly when dealing with
highly nonlinear and high-resolution problems [8]. For instance, in subsurface
physics applications such as geothermal or geologic carbon storage,
considerable uncertainties exist in parameter space, including rock properties
and reservoir management [9, 10]. To comprehensively understand the impact of
these uncertainties on operational risks like induced seismicity or
groundwater contamination, conducting over 10,000 Monte Carlo simulations may
be necessary [11, 12]. Each simulation can be time-consuming, taking hours or
even days to complete, making this procedure impractical.

Machine learning has revolutionized the field of data-driven modeling,
offering a faster and more accessible approach to tackle a wide array of
complex physical problems [13, 14, 15]. By leveraging this technology, we can
now expedite crucial workflows like large-scale inverse modeling,
optimizations, and uncertain quantification processes [16, 17, 18],
significantly boosting our understanding and decision-making capabilities.
However, data-driven modeling faces significant hurdles, including the
constant hunger for data. Acquiring or generating high-quality data can be
time-consuming, costly, and resource-intensive, mainly when dealing with high-
dimensional datasets. Furthermore, as the number of parameters increases, the
demand for training samples grows exponentially, leading to intractable
problems. It’s important to highlight another significant limitation of data-
driven modeling: it does not consider physical principles, potentially leading
to predictions that lack a physical basis. One viable approach to address this
concern is introducing a residual evaluation mechanism based on a system of
partial differential equations, as demonstrated in the concept of physics-
informed neural networks [19]. However, for the purpose of this study, our
primary focus will be on tackling this specific challenge.

One can utilize transfer learning to overcome the data-hungry challenge and
make data-driven machine learning modeling more practical and efficient.
Transfer learning is a powerful approach that leverages previously trained
models to enhance the current model’s performance [20]. There are various
techniques for transfer learning, two of which are commonly used: fine-tuning
and adding fresh layers. In fine-tuning, we take the entire set of trainable
parameters from a pre-trained model and use them as an initialization for the
current model. A relatively small learning rate is employed during the
training process to allow the model to adjust and adapt to the new task while
retaining the knowledge gained from the previous training. This approach is
practical when the pre-trained model’s features are applicable to the current
task.

Alternatively, in adding fresh layers, we utilize the pre-trained trainable
parameters and introduce new layers to the model. These fresh layers can be
added to the pre-trained model or replaced with specific layers. The pre-
trained layers are typically frozen, preventing their weights and biases from
being updated during training, while only the newly introduced layers are
trained. This method is handy when the pre-trained model captures generic
features that can be used as a starting point for the current task. While
transfer learning techniques offer significant benefits, it is essential to
acknowledge their limitations. Two notable drawbacks are (1) current
procedures only leverage knowledge from one model and apply it to the next
(i.e., one-to-one transfer) [21]. Imagine you have many (e.g., 10’s-100’s)
trained models that might be useful but can only utilize one; that would be a
considerable loss of opportunities and (2) the phenomenon known as
catastrophic forgetting, where the newly trained model completely loses the
knowledge of the previous task and cannot be applied without retraining [22].
This limitation highlights the need for progressive learning approaches that
enable knowledge transfer while preserving previously acquired information.
The survey by Zhuang et al. [21] provides a valuable resource for a
comprehensive understanding of transfer learning techniques, including their
advantages and disadvantages.

Progressive neural networks (PNN) have been developed to address the challenge
of many-to-one transfer and mitigate the issue of catastrophic forgetting,
leveraging prior knowledge from previously trained models through their
lateral connections [22]. PNN is often associated with lifelong or continual
learning, where a model continues learning from new data and tasks without
forgetting previously learned information. Imagine a scenario where we aim to
build a data-driven model that emulates high-fidelity models for single-phase
flow problems. Traditionally, we would create a training set, train our proxy
model, optimize hyperparameters, and utilize it for our desired applications
(e.g., [23, 24, 25]). However, this approach becomes inefficient when we
tackle multi-phase problems the next day, as we repeat the entire process,
disregarding the valuable knowledge gained from previously trained models.

From our perspective, the solution lies in progressive reduced order models
(built upon the PNN concept), which aim to capitalize on prior knowledge and
make the most of previously developed models, see Fig. 1. Compared with
ongoing research on transfer learning, this framework progressively
accumulates common features of physics of interest and expands them instead of
fine-tuning (including adding fresh layers) [20, 21] or enforcing loss
functions [26] resulting in immunization to catastrophic forgetting. Besides,
it can transfer learning from many trained models instead of one-to-one
transfer. The main ingredient of this framework lies within _gates_ used to
control information flow from previous generations (e.g., grandparents and
parents).

These gates are pivotal in controlling information flow while concurrently
training a child model. The total trainable parameters consist of a
combination of trainable parameters from the child model and trainable
parameters from all the gates associated with the child model, as indicated by
the red square in Fig. 1. Each gate corresponds to a connection between the
child and a preceding model (grandparents or parents in this illustration).
During the training process, if the information from a prior model proves
beneficial (i.e., it reduces the evaluation of a loss function), the
corresponding gate increases the flow of information and vice versa. In this
context, we employ a linear layer as a gate, meaning that the trainable
parameters of each gate are adjusted in response to evaluating the loss
function. It’s worth noting that one can enhance the complexity of these gates
by employing techniques such as recurrent units [27], attention mechanisms
[28], or adaptive algorithms [29].

Figure 1: A schematic figure of multi-level transfer learning. $N$ pre-trained
data-driven ML models for $N$ different problems can accelerate the training
and improve the accuracy of a new data-driven ML model for the new problems
through multi-level transfer learning techniques. E and D represent the
Encoder and Decoder, respectively. The red square represents combined
trainable parameters for the child model.

Since this framework, Fig. 1 is non-intrusive or data-driven in nature (i.e.,
it does not require a knowledge of underlying physics, governing equations, or
modifications inside high-fidelity software), it enables our proposed platform
to learn from physics simulators that can be treated as black-box. Throughout
this study, we will use the Barlow Twins reduced order model (BT-ROM) [24] to
showcase our progressive reduced order modeling framework (p-BT-ROM
hereafter). We use a series of physics-based problems ranging from transport
in porous media and gravity-driven flow to a finite deformation of a
hyperelastic material since we want to emphasize that transferring knowledge
among different physics (fluid to solid mechanics, in this case) is possible.
Besides, we employ a range of topologies from 2- to 3-Dimensional domains to
further emphasize the possibilities of transferring knowledge across different
topologies.

We aim to demonstrate that our framework effectively mitigates the data
scarcity challenge, often associated with approaches like few-shot learning
[30], within data-driven methodologies. Instead of starting anew for each
task, we systematically build upon existing models, harnessing the knowledge
and insights acquired from prior models to inform and enhance subsequent
iterations. This strategy conserves time and resources and significantly
improves the overall efficiency of data-driven machine learning. This
transformative approach expedites research and development in subsurface
physics, facilitating safer and more cost-effective decision-making for energy
production, groundwater extraction, and underground storage.

## Results

We introduce the proposed progressive learning framework that demonstrates
enhanced accuracy compared to models lacking prior knowledge (without a
parent) and achieves comparable or superior accuracy to models trained on
larger datasets (more training samples but without a parent). In this section,
we assess the performance of the progressive Barlow Twins reduced order model
(p-BT-ROM) by evaluating its results on a range of physics problems discussed
in Sections Data generation and Supplementary sec. 2.

### Data generation

We will be showcasing four physics problems, which are detailed in Table 1.
The first problem, referred to as Problem #1, involves the study of transport
in porous media with varying velocity fields. The parameter space for this
problem encompasses fluid velocity, and our primary focus is on the
concentration. Moving on to Problem #2, it revolves around transport scenarios
with varying diffusivity coefficients. In this case, the parameter space
pertains to the diffusivity coefficient, and we continue to analyze the
concentration as the quantity of interest. Problem #3 explores gravity-driven
flow in porous media, with the Rayleigh number ($\mathrm{Ra}$) as a parameter.
In this instance, the Rayleigh number defines the parameter space, and our
quantity of interest is the fluid temperature. This problem is adapted from
the research conducted by Kadeethum et al. [31]. Finally, Problem #4 addresses
the finite deformation of a hyperelastic material. The parameter space for
this problem varies according to dimensionality. For the 2-Dimensional case,
it involves the displacement at a boundary, while for the 3-Dimensional case,
it includes the displacement within the source term. In both cases, our
primary interest lies in the magnitude of solid displacement. This problem is
adapted from the work of Kadeethum et al. [32]. For comprehensive details
regarding the problem settings, quantities of interest, parameter space, and
the methodology for generating our training, validation, and testing datasets,
please refer to Section Supplementary sec. 4.

Table 1: Summary of physical problems.  Problem | Dimension | Physics | State variables | Quantity of interest | Parameters | Section
---|---|---|---|---|---|---
#1 | 2 |  | transport in
---
porous media
concentration | concentration |  | fluid
---
velocity
Supplementary sec. 2.1
#2 | 2 |  | transport in
---
porous media
concentration | concentration |  | diffusivity
---
coefficient
Supplementary sec. 2.2
#3 | 2 |  gravity-driven flow | pressure | temperature |  Rayleigh number | Supplementary sec. 2.3
velocity
temperature
#4 | 2 | solid deformation | displacement |  magnitude of displacement |  | boundary
---
condition
Supplementary sec. 2.4
3 | source term

### Investigation strategy

Throughout this Result section in the main text, we have designed the
numerical experiments to address two primary questions. The following is an
outline of the experimental setup and structure:

  1. 1.

First Question: Knowledge Transfer between Physics with Fixed Topology

     * •

Objective: Determine if the progressive learning framework can transfer
knowledge between different physics problems while maintaining a fixed
topology (i.e., same mesh and connectivity).

     * •

Parent Problems (Fluid): Problems #1, #2, and #3

     * •

Child Problem (Solid): Problem #4

     * •

Implication: Successful knowledge transfer implies the framework’s ability to
share knowledge across different physics domains.

  2. 2.

Second Question: Knowledge Sharing between Different Topologies

     * •

Objective: Investigate whether the progressive learning framework can share
knowledge between different topologies.

     * •

Parent Problems (2-Dimensional Topology): Problems #1, #2, #3, and #4

     * •

Child Problem (3-Dimensional Topology): Problem #4

     * •

Implication: Successful knowledge sharing indicates the framework’s potential
to transfer knowledge across different dimensions and topologies.

To ensure a systematic approach, we divide the investigation into two stages:

  * •

Stage 1: Similar Topology (Section Similar Topology and Supplementary sec.
5.5)

Concentration on problems with different physics but the same topology. Focus
on evaluating the effects of the number of training samples and the number of
parents on the performance of the Barlow Twins reduced order model (BT-ROM)
[24] and p-BT-ROM.

  * •

Stage 2: Different Topologies (Section Different Topologies)

An investigation involving problems with different physics and different
topologies. Emphasis on analyzing the impact of the number of training samples
and the number of parents on the performance of BT-ROM and p-BT-ROM. Detailed
descriptions of each topology are provided in Supplementary sec. 1. By
following this systematic structure, we aim to provide clear insights into the
progressive learning framework’s capabilities regarding knowledge transfer and
sharing across various physics domains, topologies, and dimensions.

In addition to our research findings in the main text, we offer comprehensive
numerical examples in supplementary sections to enhance the understanding of
our work. Specifically, Section Supplementary sec. 5.1 demonstrates the
conventional approach of transfer learning, where pre-trained models’ weights
and biases are utilized to initialize a new model (referred to as "init.")
[20, 21]. This approach serves as a benchmark for comparison with our
progressive learning framework, known as p-BT-ROM, presented in Section
Supplementary sec. 5.2. Furthermore, we explore the potential benefits of
combining initialization with p-BT-ROM, denoted as "p-BT-ROM with init."
Through our analysis in Section Supplementary sec. 5.4, we conclude that the
preferred approach going forward is to employ p-BT-ROM with init. To
supplement our findings in Section Similar topology, we provide additional
results of p-BT-ROM with init. for Problems #1, #2, #3, and #4 with fixed
topology in Section Supplementary sec. 5.5. This supplementary section sheds
light on the performance of our framework when applied to different problem
instances. Additionally, we have thoroughly investigated the impact of varying
the number of parents within our framework to assess its effectiveness,
yielding valuable insights.

### Similar topology

We will illustrate that progressive learning enables us to transfer knowledge
to problems with different underlying physics (e.g., fluid to solid mechanics)
and parameter spaces (e.g., boundary conditions or material properties). We
present a sample of a reconstruction loss (Supplementary eq. 14) result of the
validation set (as a function of the number of epochs and parents) for Problem
#4 in a 2-Dimensional setting (the detail of this geometry is shown in
Supplementary fig. 1a) in Fig. 2. We note that more comprehensive results can
be found in Supplementary sec. 5.3, but here, we want to introduce our
framework’s benefit briefly for the training phase of the (p-)BT-ROM (see
Supplementary fig. 2 \- third step). To reiterate, our approach involves a
combination of p-BT-ROM, utilizing the weights and biases of the trained
models as an initialization. We provide in-depth results and rationale for
adopting this approach in Sections Supplementary sec. 5.1, Supplementary sec.
5.2, Supplementary sec. 5.3, and Supplementary sec. 5.4. From this figure, the
more parents we provide the p-BT-ROM model, the more accurate it becomes. We
note that with 0 parents, the p-BT-ROM is the original BT-ROM [24]. These
results imply that the p-BT-ROM could transfer knowledge from fluid problems
(Problems #1 to #3) to a solid problem (Problem #4). We note that as we add
more parents to the framework, our trainable parameters inevitably become
larger (see Supplementary sec. 6). These additional parameters lead to a
heavier computational time and memory load, but the cost remains far below
running the physical model.

Figure 2: Illustration that our framework can transfer knowledge from fluid
mechanics (Problems #1-#3) to solid mechanics (Problem #4) problems.
Reconstruction loss ((Supplementary eq. 14)) of the validation set for Problem
#4 (Supplementary sec. 2.4) with 2-Dimensional domain (Supplementary fig. 1a).
This loss reflects the training phase of the (p-)BT-ROM (see Supplementary
fig. 2 \- third step). We note that for Problem #4, we show the results in the
deformed configuration.

To summarize our testing process, we train our model with a specific number of
epochs (100 in this case, however, one should always perform trial and error
for any new physics-based problems or new neural network architecture.) with
our _training set_ (see Supplementary fig. 2 \- third step). Subsequently, we
pick a set of weights and biases that deliver the lowest loss for a
_validation set_ (still Supplementary fig. 2 \- third step). Finally, we use a
decoder of p-BT-ROM with the set of weights and biases in combination with
radial basis function (RBF) interpolator for the evaluation using a _test set_
(see Supplementary fig. 2 \- fourth and fifth steps). We present the test
sets’ results in Supplementary sec. 5.5. In a nutshell, the results of the
test samples are similar to the behavior we observed in the validation cases
(Fig. 2). To elaborate, as we increase the number of parents, our p-BT-ROM
gains more accuracy. As expected, the BT-ROM (0 parent) has the worst
accuracy. However, even the model with many parents still could not achieve
the same level of accuracy as the model that trained with a very large
training set with 0 parent (BT-ROM). This behavior implies that there are
benefits of adding training samples that progressive neural networks might not
be able to compensate for or substitute. On the other hand, they can produce
more accurate results for fixed-sized datasets (supplemented with data from
the simpler parent problems).

### Different topologies

In this study, we demonstrate the effectiveness of the progressive approach in
transferring knowledge to problems with significantly different topologies.
Specifically, we investigate the transfer from 2-dimensional to 3-dimensional
domains or between 2-dimensional domains with varying numbers of holes (1 to 3
holes). We assess the validation reconstruction loss ((Supplementary eq. 14))
as a function of the number of epochs and parents for Problem #4 in a
3-dimensional context (see Supplementary fig. 1d), shown in Fig. 3. Consistent
with our earlier observations, we find that increasing the number of parents
provided to the p-BT-ROM model improves its accuracy. Notably, there is a
minimal discrepancy between the model’s performance with 3 parents and that
with 4 parents. These results indicate the feasibility of transferring
knowledge across different topologies and dimensions, ultimately enhancing the
precision of our data-driven framework. Additionally, as discussed in the
previous section on similar topology, incorporating more parents into the
framework expands the trainable parameter space, increasing training costs.
For a comprehensive understanding of the number of parameters involved, please
refer to Supplementary sec. 6.

Figure 3: Illustration that our framework can transfer knowledge from
2-Dimensional to 3-Dimensional problems. Reconstruction loss ((Supplementary
eq. 14)) of the validation set for Problem #4 (Supplementary sec. 2.4) with
3-Dimensional domain (Supplementary fig. 1d). This loss reflects the training
phase of the (p-)BT-ROM (see Supplementary fig. 2 \- third step). We note that
for Problem #4, we show the results in the deformed configuration.

From Section Similar topology, we have illustrated that our p-BT-ROM can
improve the accuracy of BT-ROM by adding more parents to the model. However,
our investigation is limited to one mesh; see Supplementary fig. 1a, for two
physics problems, Problems #3 to #4. Here, we will investigate the p-BT-ROM’s
performance using different meshes and topologies, see Supplementary figs.
1b-c for three physics problems, Problems #2 to #4 (test model performance on
different physics and topologies). The schematic of p-BT-ROM specifies
child$-$parent(s) relationship used throughout this section, and its
subsections are shown in Fig. 4. In short, we use the p-BT-ROM with the small
dataset from Supplementary sec. 5.5.2. To elaborate, for our 4 parents, we use
Problem #1 (Supplementary sec. 2.1) with a training set of
$\mathrm{M}N_{t}=2405$ ($\mathrm{M}=5$), Problem #2 (Supplementary sec. 2.2)
with a training set of $\mathrm{M}N_{t}=2405$ ($\mathrm{M}=5$), Problem #3
(Supplementary sec. 2.3) with a training set of $\mathrm{M}N_{t}=3623$
($\mathrm{M}=5$), and Problem #4 (Supplementary sec. 2.4) with a training set
of $\mathrm{M}=100$. This is representative of many realistic scenarios where
the model for which we need a ROM can only produce a small training set with
$M=100$ samples, but larger datasets are available for simpler problems. We
want to emphasize that we can combine different types of physics problems
(i.e., fluid and solid mechanics problems) and different characteristics
(i.e., transient and steady-state problems) into our pool of parents.

Figure 4: Different topologies - schematic of p-BT-ROM specifies
child$-$parent(s) relationship. See Supplementary fig. 1 for more information
on the geometries.

#### Transport in porous media - Problem #2

Here, we focus on a physics problem of transport in porous media - Problem #2
(Supplementary sec. 2.2) with mesh shown in Supplementary fig. 1b. A sample of
snapshots is presented in Supplementary fig. 4a. We have one primary variable,
saturation. We use a small training set of $\mathrm{M}N_{t}=1443$
($\mathrm{M}=3$), large training set of $\mathrm{M}N_{t}=9620$
($\mathrm{M}=20$), and a test set of $\mathrm{M_{test}}N_{t}=4810$
($\mathrm{M_{test}}=10$). We present the mean squared error (MSE),
Supplementary eq. 19, and mean absolute error (MAE), Supplementary eq. 20,
results in Figs. 5a-b. From these plots, the p-BT-ROM with 4 parents with the
small training set performs much better than the BT-ROM (0 parent) with the
small training set. Besides, the p-BT-ROM with 4 parents’ performance is very
similar to that of the BT-ROM (0 parent) with the large training set while
using eight times fewer training samples. We present the MAE results as a
function of parameter $\bm{\mu}$ in Figs. 15a-b. The pattern of these three
figures is similar (i.e., the shape), but the magnitude of the MAE is
different. This observation implies that the p-BT-ROM only alters the
magnitude of the error, not its pattern.

Figure 5: Different topologies - (a) mean squared error (MSE) and (b) mean
absolute error (MAE) results for Problem #2 (Supplementary sec. 2.2) using the
topology shown in Supplementary fig. 1b (for $\mathrm{M}N_{t}=1443$
($\mathrm{M}=3$) and $\mathrm{M}N_{t}=9620$ ($\mathrm{M}=20$)), (c) mean
squared error (MSE) and (d) mean absolute error (MAE) results for Problem #3
(Supplementary sec. 2.3) using the topology shown in Supplementary fig. 1c
(for $\mathrm{M}N_{t}=2706$ ($\mathrm{M}=5$) and $\mathrm{M}N_{t}=10832$
($\mathrm{M}=20$)), and (e) mean squared error (MSE) and (f) mean absolute
error (MAE) results for Problem #4 (Supplementary sec. 2.4) using the topology
shown in Supplementary fig. 1d. We present the small training set in red for 0
parent and in green for 4 parents. Solid black lines represent the large
training set. Dashed black lines show accuracy gain by including knowledge
from parents: p-BT-ROM.

#### gravity-driven in porous media - Problem #3

Next, we focus on a physics problem of gravity-driven in porous media -
Problem #3 (Supplementary sec. 2.3) with mesh shown in Supplementary fig. 1c.
A sample snapshot is illustrated in Supplementary fig. 4b. Here, we have three
primary variables, fluid pressure, velocity, and temperature. We only focus on
the fluid temperature in this case. We use a small training set of
$\mathrm{M}N_{t}=2706$ ($\mathrm{M}=5$), a large training set of
$\mathrm{M}N_{t}=10832$ ($\mathrm{M}=20$), and a test set of
$\mathrm{M_{test}}N_{t}=5346$ ($\mathrm{M_{test}}=10$). We present the MSE and
MAE results in Figs. 5c-d. Again, the p-BT-ROM with 4 parents trained with the
small training set delivers much better results than those of the BT-ROM (0
parent) trained with the small training set. Moreover, the p-BT-ROM with 4
parents trained with the small training set outperforms the BT-ROM (0 parent)
trained with the large training set (approximately five times larger). Similar
to the previous section, the pattern of MAE results as a function of parameter
$\bm{\mu}$ are also similar between the p-BT-ROM (4 parents) and BT-ROM (0
parent), see Supplementary figs. 15c-d.

#### Hyperelasticity problem - Problem #4

Lastly, we move on to a physics problem of finite deformation of hyperelastic
material - Problem #4 (Supplementary sec. 2.4) with mesh shown in
Supplementary fig. 1d, one sample snapshot is shown in Supplementary fig. 4c.
We want to emphasize that not only do we test a setting where all parents are
fluid mechanics problems but the child is solid mechanics problems, but we
also want to illustrate that we can use 2-Dimensional meshes for the parents
but 3-Dimensional mesh for the child. We have only one primary variable,
displacement, and we focus on its magnitude of it. We use a small training set
of $\mathrm{M}=100$, a large training set of $\mathrm{M}=1600$, and a test set
of $\mathrm{M_{test}}=100$. We present the MSE and MAE results in Figs. 5e-f
and the pattern of MAE results as a function of parameter $\bm{\mu}$ in
Supplementary figs. 15e-f. Again, the p-BT-ROM with 4 parents trained with the
small training set delivers much better results than those of the BT-ROM (0
parent) trained with the small and large training sets. The pattern of MAE
results as a function of parameter $\bm{\mu}$ are also not much different
between the p-BT-ROM (4 parents) and BT-ROM (0 parent) - see Supplementary
figs. 15e-f (i.e., only the shape of the distribution but the magnitude of the
MAE is drastically different.).

## Discussion

By conducting a series of numerical experiments, we have successfully
demonstrated the effectiveness of a progressive reduced order model. The core
element of this framework is centered around gates used to regulate the flow
of information from previously trained models. These gates are crucial in
managing information flow while simultaneously training the current model. The
total trainable parameters encompass a combination of the trainable parameters
from the current model and those from all the gates associated with it. Each
gate corresponds to a connection between the current model and a preceding
model, with a dedicated gate for each connection. During the training process,
if information from a prior model proves advantageous, meaning it reduces the
loss function, the corresponding gate increases the flow of information.
Conversely, if the information is detrimental, the gate reduces the flow. In
this context, we employ a linear layer as a gate, which means that the
$\mathbf{W}$ and $\mathbf{b}$ parameters of each gate are adjusted in response
to evaluating the loss function. It’s important to note that the complexity of
these gates can be enhanced by utilizing techniques such as recurrent units
[27], attention mechanisms [28], or adaptive algorithms [29].

The outcomes of our experiments are succinctly presented in Tab.
LABEL:tab:summary_results. These results provide compelling evidence that
leveraging the knowledge gained from previously trained models offers a
promising solution to address the limited availability of training data.
Notably, our framework, which incorporates 4 parent models, surpasses its
counterpart that lacks such parental guidance, even when the latter is trained
with 9 times more data. This evidence highlights the remarkable potential of
utilizing prior knowledge from trained models to enhance performance.
Moreover, our research showcases the feasibility of transferring knowledge
across diverse physics domains, such as fluid mechanics, to solid mechanics
problems. Additionally, we demonstrate the applicability of our approach to
different domain topologies, including transitions from 2-Dimensional to
3-Dimensional domains and variations in the number of holes inside a
2-Dimensional domain. Our findings underscore the considerable benefits of
harnessing previously acquired knowledge within a progressive reduced order
model, enabling improved performance and knowledge transfer across various
domains and topologies.

Table 2: Summary of results for test sets.  setting | num. parent(s) | training samples | $\operatorname{avg.}{\mathrm{MAE}}$ | $\operatorname{std.}{\mathrm{MAE}}$ | remark
---|---|---|---|---|---
similar topology | Problem #3 | 0 |  $\mathrm{M}=40$ ($\mathrm{M}N_{t}=28966$) | 2.41E-03 | 1.27E-03 | large training set
0 |  $\mathrm{M}=5$ ($\mathrm{M}N_{t}=3623$) | 1.61E-02 | 5.83E-03 | small training set
1 |  $\mathrm{M}=5$ ($\mathrm{M}N_{t}=3623$) | 5.87E-03 | 2.75E-03 | from transport to gravity-driven problems
2 |  $\mathrm{M}=5$ ($\mathrm{M}N_{t}=3623$) | 4.82E-03 | 2.36E-03
Problem #4 | 0 | $\mathrm{M}=1600$ | 1.08E-04 | 4.84E-05 | large training set
0 | $\mathrm{M}=100$ | 3.69E-03 | 1.93E-03 | small training set
1 | $\mathrm{M}=100$ | 3.79E-04 | 7.19E-05 | from fluid to solid problems
2 | $\mathrm{M}=100$ | 3.70E-04 | 6.52E-05
3 | $\mathrm{M}=100$ | 3.20E-04 | 4.40E-05
different topologies | Problem #2 | 0 |  $\mathrm{M}=20$ ($\mathrm{M}N_{t}=9620$) | 3.99E-03 | 1.64E-03 | large training set
0 |  $\mathrm{M}=3$ ( $\mathrm{M}N_{t}=1443$) | 2.68E-02 | 7.77E-03 | small training set
4 |  $\mathrm{M}=3$ ( $\mathrm{M}N_{t}=1443$) | 4.01E-03 | 2.48E-03 | from 1 hole to 2 holes topologies
Problem #3 | 0 |  $\mathrm{M}=20$ ($\mathrm{M}N_{t}=10832$) | 9.77E-03 | 5.25E-03 | large training set
0 |  $\mathrm{M}=5$ ($\mathrm{M}N_{t}=2706$) | 4.06E-02 | 1.32E-02 | small training set
4 |  $\mathrm{M}=5$ ($\mathrm{M}N_{t}=2706$) | 3.84E-03 | 2.16E-03 | from 1 hole to 3 holes topologies
Problem #4 | 0 | $\mathrm{M}=1600$ | 2.64E-03 | 1.73E-03 | large training set
0 | $\mathrm{M}=100$ | 6.66E-03 | 3.17E-03 | small training set
4 | $\mathrm{M}=100$ | 8.79E-04 | 3.61E-04 | from 2- to 3-Dimensional topologies

However, there is an additional computational cost associated with this
benefit. As the number of parents increases, the number of trainable
parameters also grows, which requires more computational resources in terms of
time and space complexities (Supplementary sec. 6). It’s important to note
that this additional cost primarily affects the training process (the
prediction or inferring phase remains speedy evaluation). For the specific
case illustrated in Fig. 2 and Tab. 3, when training the model with 0 and 3
parents using NVIDIA Quadro RTX 8000, there is an approximate increase of 10%
in the wall time required. However, the prediction time remains very similar,
with both cases taking less than a second. In contrast, running a single high-
fidelity model for Problem #3 would typically take at least 30 minutes.

It is worth mentioning that the framework we have developed, known as
progressive reduced order model concept, can be applied to various machine-
learning architectures and is not limited to the Barlow Twins reduced order
model (p-BT-ROM) alone. Moreover, our framework is not constrained to any
specific application domain, such as fluid or solid mechanics, and can be
utilized for any physical problem. Additionally, it is adaptable to different
sources of data, including high-fidelity numerical simulations or field
measurements. This flexibility is highly advantageous as the framework can
accumulate knowledge from diverse sources, enhancing its overall performance.

Even though this study illustrates the potential of growing our data-driven
model via the concept of progressive neural networks [22], six limitations
need to be acknowledged and discussed. The first one is more complex
geometries and topologies should be used to test this framework. Our framework
does not have direct information on topologies at this stage. Still, one could
use an additional network such as (graph-)trunk nets [33] to embed the
topologies of each parent and child into the framework. Secondly, as mentioned
in Supplementary sec. 4 we uniformly sample our parameter space $\bm{\mu}$;
subsequently, one might argue that if we perform a smart sampling or adaptive
sampling [34], we might also end up using less training samples similarly to
p-BT-ROM. Hence, one also wants to combine our progressive framework with an
adaptive sampling technique to reduce training samples further (or maintain
the number of training samples but improve the models’ accuracy).

Thirdly, it is worth noting that all the physics problems discussed here,
labeled as Problems #1 to #4, operate on the same physical scale in terms of
computational size. However, future investigations should consider
incorporating information from different scales, such as integrating pore-
scale and field-scale data. To control the flow of information across these
scales, adapting an attention mechanism [28] could prove beneficial.

Fourthly, it is essential to acknowledge that increasing the number of parents
in our framework leads to a higher number of trainable parameters, as
explained in Supplementary sec. 6. This results in a more intricate and
computationally expensive training process. To address this challenge, a
potential solution could involve employing a reinforcement learning technique
to choose parents that are relevant for the current generation selectively.
This approach would help optimize the selection of parents, ensuring that only
the most useful models contribute to training the reduced order model (ROM)
for a specific application, such as $\mathrm{CO_{2}}$ storage.

For instance, in Fig. 6, we have multiple trained models available in our
inventory. In the case of developing a ROM for $\mathrm{CO_{2}}$ storage, a
straightforward approach would be to utilize all the trained models (i.e., six
parents) as parents for the $\mathrm{CO_{2}}$ storage ROM. However,
considering the increased number of trainable parameters with more parents, a
reinforcement learning (or meta-learning) framework can be employed to
determine which parents are truly influential for the $\mathrm{CO_{2}}$
storage ROM. By selectively choosing and utilizing only the relevant ROMs as
parents, such as excluding the ROM for the contact problem if it does not
contribute significantly to the performance of the $\mathrm{CO_{2}}$ storage
ROM, we can optimize the training process and streamline the computational
resources.

Figure 6: A hypothetical situation where we want to develop a ROM of
$\mathrm{CO_{2}}$ storage and have trained ROMs of linear elasticity, contact
problem, plasticity, transport problem, gravity-driven flow, and multiphase
flow in our inventory that we can pick as (grand)parent(s) of
$\mathrm{CO_{2}}$ storage ROM. We use red to represent fluid problems and grey
to show solid problems. Lighter shades are used to represent grandparents,
while darker ones are for parents.

The fourth challenge also involves reaching knowledge saturation. Based on our
observations, simply increasing the number of parents does not consistently
lead to enhanced model performance. To illustrate, as shown in Fig. 3, having
three or four parents yields roughly equivalent results. A potential solution
to this issue could involve employing a meta-learning framework (as well as
reinforcement learning) to select relevant parents.

The fifth challenge revolves around integrating physical information, either
by including an additional loss function or by modifying the network
architecture itself. This entails finding ways to embed fundamental physical
principles into the learning process. Given that our framework strictly relies
on data, it might not yield a genuinely physical solution and could generate
predictions lacking a physical basis. To overcome this issue for each
individual component (each child or parent), one can employ the concept of
physics-informed neural networks [19]. However, when dealing with a
combination of parents and children, each influenced by different sets of
physical phenomena, this framework may not be applicable. In such cases, it
might be beneficial to explore transfer learning techniques using physics-
guided models [35].

Lastly, it is important to highlight that our progressive learning framework
follows a strict hierarchical structure. To elaborate, consider Fig. 6. If we
select ROM: multiphase flow as a parent for the ROM: $\mathrm{CO_{2}}$
storage, it automatically implies that ROM: transport problem becomes a
grandparent. In such a scenario, if we decide to retrain ROM: transport
problem, it necessitates the retraining of ROM: multiphase flow, followed by
ROM: $\mathrm{CO_{2}}$ storage. This strict hierarchical growth can introduce
complexities and dependencies in the training process.

To address this challenge, one possible approach is to utilize a directed
acyclic graph (DAG) that allows for information exchange among different
sources without assuming a strict hierarchical connection. The use of a DAG
can provide more flexibility in incorporating information from multiple models
or sources, eliminating the need for a fixed hierarchical structure. This
approach has been explored in the context of machine learning models, as
demonstrated in the work of Gorodetsky et al. [36]. By adopting a DAG-based
approach, we can mitigate the problem of strict hierarchical dependencies and
enhance the adaptability and scalability of our framework.

## Methodology

We propose the use of progressive learning applied in the context of
progressive ROMs. We note our framework can be applied to any data-driven
ROMs; however, here in our application, each machine learning model
corresponds to the Barlow-Twins reduced order model (BT-ROM), as described in
Supplementary sec. 3.1 and detailed in [24]. In Fig. 7, we present our
progressive BT-ROM (p-BT-ROM) consisting of three columns. Each column
represents a model, with two previously trained models referred to as parent 1
($c_{1}$) and parent 2 ($c_{2}$), and the current model referred to as the
child ($c_{3}$). It is important to note that in practice, the p-BT-ROM can be
extended to a finite number of columns based on the available machine memory.
The relationship between the parents and child in the p-BT-ROM is not limited
to a hierarchical structure, where parent 2 is the child of parent 1. The
parents can also be independent, forming standalone BT-ROMs without lateral
connections. However, the child can still establish lateral connections with
each parent individually, building separate relationships with them.

By employing progressive learning in the p-BT-ROM, we facilitate the transfer
of knowledge and information from the previously trained models to enhance the
performance and learning capability of the current model. This approach allows
us to leverage prior knowledge effectively and overcome the limitations of
one-to-one transfer and catastrophic forgetting.

Figure 7: A schematic of p-BT-ROM. Assuming we have three columns, $c_{1}$,
$c_{2}$, and $c_{3}$ representing parent 1, parent 2, and child, respectively.
Each column represents a BT-ROM; see Supplementary sec. 3.1 for a detailed
description. Each BT-ROM consists of one encoder, one decoder, and one
projector. Assuming each one of them has only three layers, $r_{1}$, $r_{2}$,
and $r_{3}$. According to Eq. (1), for child $c_{3}$, $r_{3}$, we have a
contribution from parent 1 ($c_{1}$) and parent 2 ($c_{2}$) as
$U_{r_{3}}^{c_{3}:c_{1}}h_{r_{2}}^{c_{1}}+U_{r_{3}}^{c_{3}:c_{2}}h_{r_{2}}^{c_{2}}$,
and child $c_{3}$, $r_{2}$, we have a contribution from parent 1 ($c_{1}$) and
parent 2 ($c_{2}$) as
$U_{r_{2}}^{c_{3}:c_{1}}h_{r_{1}}^{c_{1}}+U_{r_{2}}^{c_{3}:c_{2}}h_{r_{1}}^{c_{2}}$.
The reduced manifolds plots are constructed using t-Distributed Stochastic
Neighbor Embedding (t-SNE). Different colors represent each value of
$\bm{\mu}$ value. We calculate the t-SNE plots using the Scikit-Learn package
using its default setting and perplexity of 15.

Each BT-ROM consists of one encoder, one decoder, and one projector, as shown
in Fig. 7. Each of them has different numbers of layers or rows (i.e.,
$r_{1},r_{2},r_{3},\dots,r_{n}$). We use $\mathcal{N}\left(\cdot\right)$ to
represent any nonlinear activation functions (e.g., tanh, Sigmoid, or ReLU).
$\mathcal{L}\left(\cdot\right)$ is a linear layer within each column (i.e.,
take the output from the previous layer within the same column,
$h_{r_{i}-1}^{c_{k}}$, and perform a linear transformation with its
$\mathbf{W}$ and $\mathbf{b}$). ${U}\left(\cdot\right)$ is also a linear layer
that acts as a gate for each lateral connection (i.e., take the output from
the previous layer of the previous columns, $h_{r_{i}-1}^{c_{j}}$ where $j<k$
and $k$ is the current column, and perform a linear transformation with its
$\mathbf{W}$ and $\mathbf{b}$). To this end, our output from the $k$ column is
read as follows

$h_{r_{i}}^{c_{k}}=\mathcal{N}\left(\mathcal{L}_{r_{i}}^{c_{k}}\left(h_{r_{i}-1}^{c_{k}}\right)+\sum_{j<k}\left({U}_{r_{i}}^{{c_{k}}:{c_{j}}}\left(h_{r_{i}-1}^{c_{j}}\right)\right)\right),$
(1)

and for the first column (i.e., $c_{k}=0$)

$h_{i}^{c_{k}=0}=\mathcal{N}\left(\mathcal{L}_{i}^{c_{k}=0}\left(h_{i-1}^{c_{k}=0}\right)\right).$
(2)

From these equations, one can see that if information from the previously
trained models is not useful, ${U}\left(\cdot\right)$ will be zero out its
output during the training process. Besides, this framework can be used to
tackle its previous tasks by simply focusing its output from the focus column
(i.e., immune to catastrophic forgetting). To elaborate, assuming we want to
predict a job that parent 2 is trained for, we only need to ignore the outputs
from parent 1 and child and only consider the output from parent 2. We note
that each component of p-BT-ROM acts independently (i.e., there are lateral
connections within the encoder but no lateral connections among the encoder,
decoder, and projector).

## Acknowledgments

This work is supported by the US Department of Energy Office of Fossil Energy
and Carbon Management project - Science-Informed Machine Learning to
Accelerate Real-Time Decisions in Subsurface Applications (SMART) initiative.
Sandia National Laboratories is a multi-mission laboratory managed and
operated by National Technology & Engineering Solutions of Sandia, LLC
(NTESS), a wholly owned subsidiary of Honeywell International Inc., for the
U.S. Department of Energy’s National Nuclear Security Administration
(DOE/NNSA) under contract DE-NA0003525. This written work is authored by an
employee of NTESS. The employee, not NTESS, owns the right, title, and
interest in and to the written work and is responsible for its contents. Any
subjective views or opinions that might be expressed in the written work do
not necessarily represent the views of the U.S. Government. The publisher
acknowledges that the U.S. Government retains a non-exclusive, paid-up,
irrevocable, world-wide license to publish or reproduce the published form of
this written work or allow others to do so, for U.S. Government purposes. The
DOE will provide public access to results of federally sponsored research in
accordance with the DOE Public Access Plan.

## Author contributions statement

T. Kadeethum: Conceptualization, Formal analysis, Software, Validation,
Writing - original draft, Writing - review & editing. D. O’Malley:
Conceptualization, Formal analysis, Supervision, Validation, Writing - review
& editing. Y. Choi: Conceptualization, Formal analysis, Supervision,
Validation, Writing - review & editing. H.S. Viswanathan: Conceptualization,
Supervision, Writing - review & editing. H. Yoon: Conceptualization, Formal
analysis, Funding acquisition, Supervision, Writing - review & editing.

## Competing interests

The authors declare no competing interests

## Supplementary information

##  Supplementary sec. 1 A collections of geometries used in this manuscript

Here, we summarize all meshes used throughout this manuscript in Supplementary
fig. 1 and Supplementary tab. 1. In short, Supplementary fig. 1a is a
2-Dimensional domain with a size of $2\times 1$ and has a hole inside the
domain. It has 1503 nodes and is used for all physics problems ( Supplementary
sec. 2.1 to Supplementary sec. 2.4). Supplementary fig. 1b is also a
2-Dimensional domain with a size of $2\times 1$ and has two holes inside the
domain. It has 1698 nodes and is used for Problem #2 (Supplementary sec. 2.2).
Supplementary fig. 1c is a 2-Dimensional domain with a size of $2\times 1$ and
has three holes inside the domain. It has 620 nodes and is used for Problem #3
(Supplementary sec. 2.3). Supplementary fig. 1d is a 3-Dimensional domain with
a size of $1\times 1\times 1$ and has no hole inside the domain. It has 1330
nodes and is used for Problem #4 (Supplementary sec. 2.4). We note that the
details of each physics problem can be found in Supplementary sec. 2.

Supplementary fig. 1: Illustrations of geometries and their corresponding mesh used in this manuscript. The details of each mesh can be found in Table 1. Supplementary tab. 1: Summary of main information for each mesh. | Dimension | Size |  | Number of
---
nodes
Used for
| A rectangular
---
with 1 hole
2 | 2 x 1 | 1503 | Problems #1, 2, 3, 4
| A rectangular
---
with 2 holes
2 | 2 x 1 | 1698 | Problem #2
| A rectangular
---
with 3 holes
2 | 2 x 1 | 620 | Problem #3
A cube | 3 | 1 x 1 x 1 | 1330 | Problem #4

##  Supplementary sec. 2 A collections of physical problems used in this
manuscript

For all problems, let $\Omega\subset\mathbb{R}^{d}$ ($d\in\\{1,2,3\\}$) denote
the computational domain and $\partial\Omega$ denote the boundary.
$\partial\Omega$ can be represented by the Dirichlet boundary condition
($\partial\Omega_{D}$) or the Neumann boundary condition
($\partial\Omega_{N}$). For Problems # 1 to 3, the time domain is denoted by
$\mathbb{T}=\left(0,\tau\right]$ with $\tau>0$ (i.e., $\tau$ is the final
time). We note that Problem # 4 is a steady-state problem.

###  Supplementary sec. 2.1 Problem #1: transport problem with velocity as a
parameter

In this problem, we focus on the transport problem, and we have a
concentration $c:\Omega\times\mathbb{T}\rightarrow\mathbb{R}$
($\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$)
as a primary variable. Our governing equations are shown below

$\begin{split}\frac{\partial}{\partial t}(\phi
c)+\nabla\cdot\left(\bm{q}c\right)-\nabla\cdot\left(\bm{D}\nabla
c\right)=q&\mbox{ \> in \> }\Omega\times(0,\mathbb{T}],\\\
\eta(\bm{q},c)\cdot\mathbf{n}={c_{in}}\bm{q}\cdot\mathbf{n}&\mbox{ \> on \>
}\partial\Omega_{{\rm in}}\times(0,\mathbb{T}],\\\ {\bm{D}}\nabla
c\cdot\mathbf{n}=0&\mbox{ \> on \> }\partial\Omega_{{\rm
out}}\times(0,\mathbb{T}],\\\ c={c_{0}}&\text{ \> in \> }\Omega\text{ at
}t=0,\end{split}$ (Supplementary eq. 1)

where $\phi$ is a porosity, $\bm{D}$ is a tensor of diffusivity coefficient,
$q$ represents a source term, $c_{in}$ is the inflow concentration, $c_{0}$ is
the initial concentration, $\mathbf{n}$ is a normal unit vector to each
element surface, the mass flux $\eta(\bm{q},c)$ is defined as

$\eta(\bm{q},c):=\bm{q}c-{\bm{D}}(\phi)\nabla c.$ (Supplementary eq. 2)

$\bm{q}$ is the superficial velocity vector defined by
$\bm{q}=(20.0,\mu_{1}\cos(\pi x))$. To approximate $c_{h}$, we use the first-
order discontinuous Galerkin approximation, and the detail of discretization
and the FOM source codes can be found in [37, 31]. Here, we have one parameter
$\bm{\mu}=(\mu_{1}=[5.0,25.0])$. The rest parameters are fixed as $\phi=0.3$,
$\bm{D}=(0.2,0.0,0.0,0.2)$ (i.e., homogeneous and isotropic), and we enforce
$c_{0}=0.0$ and $c_{in}=1.0$ on the left boundary. We use geometries and
meshes shown in Supplementary fig. 1.

###  Supplementary sec. 2.2 Problem #2: Transport problem with diffusivity
coefficient as a parameter

Here, we also use the same governing equations shown in (Supplementary eq. 1).
For fixed parameters, we use $\phi=0.3$, $\bm{q}=(20.0,\cos(\pi x))$, and we
enforce $c_{0}=0.0$ and $c_{in}=1.0$ on the left boundary. The parameter
$\bm{\mu}=(\bm{D}=(\mu_{1},0.0,0.0,\mu_{1}))$, again, homogeneous and
isotropic, and $\mu_{1}=[0.1,1.0]$. We use geometries and meshes shown in
Figs. 1a, b.

###  Supplementary sec. 2.3 Problem #3: gravity-driven problem with Rayleigh
number as a parameter

For Problem #3, we use a gravity-driven flow in porous media. We follow the
FOM used in [31, 24], and here we briefly discuss the governing equations of
this problem. Essentially, we deal with coupled PDEs of mass conservation and
heat advection-diffusion equations. As presented in [31, 24], we write our
system of equations in dimensionless form, and primary variables in this
setting are $\bm{u}(\cdot,t):\Omega\times\mathbb{T}\to\mathbb{R}^{d}$, which
is a vector-valued Darcy velocity
($\mathrm{d}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{s}$),
$p(\cdot,t):\Omega\times\mathbb{T}\to\mathbb{R}^{d}$, which is a scalar-valued
fluid pressure
($\mathrm{d}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{s}$),
and $T(\cdot,t):\Omega\times\mathbb{T}\to\mathbb{R}^{d}$, which is a scalar-
valued fluid temperature
($\mathrm{d}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{s}$).
Time is denoted as $t$
($\mathrm{d}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{s}$).
The system of our governing equations read

$\begin{split}\bm{u}+\nabla p-\mathbf{y}\operatorname{Ra}T=0,&\text{ \> in \>
}\Omega\times\mathbb{T},\\\ \nabla\cdot\bm{u}=0,&\text{ \> in \>
}\Omega\times\mathbb{T},\\\ p=p_{D}&\text{ \> on \>
}\partial\Omega_{p}\times\mathbb{T},\\\ \bm{u}\cdot\mathbf{n}=q_{D}&\text{ \>
on \>}\partial\Omega_{q}\times\mathbb{T},\\\ p=p_{0}&\text{ \> in \>
}\Omega\text{ at }t=0,\end{split}$ (Supplementary eq. 3)

for the mass balance equation where $\partial\Omega_{p}$ and
$\partial\Omega_{q}$ are the prescribed pressure and flux boundaries,
respectively, and $\mathrm{Ra}$ is the Rayleigh number

$\mathrm{Ra}:=\frac{g\alpha\kappa\Delta T^{*}H}{K}.$ (Supplementary eq. 4)

Subsequently, the heat advection-diffusion equation is described as

$\begin{split}\frac{\partial T}{\partial t}+\bm{u}\cdot\nabla
T-\nabla^{2}T-f_{c}=0,&\mbox{ \> in \> }\Omega\times(0,\mathbb{T}],\\\
T=T_{D}&\mbox{ \> on \> }\partial\Omega_{{T}}\times(0,\mathbb{T}],\\\
(-\bm{u}T+\nabla T)\cdot\mathbf{n}={T_{\rm in}}\bm{u}\cdot\mathbf{n}&\mbox{ \>
on \> }\partial\Omega_{{\rm in}}\times(0,\mathbb{T}],\\\ \nabla
T\cdot\mathbf{n}=0&\mbox{ \> on \> }\partial\Omega_{{\rm
out}}\times(0,\mathbb{T}],\\\ T={T_{0}}&\text{ \> in \> }\Omega\text{ at
}t=0,\end{split}$ (Supplementary eq. 5)

where $\partial\Omega_{{T}}$ is prescribed temperature boundary,
$\partial\Omega_{\rm in}$ and $\partial\Omega_{\rm out}$ denote inflow and
outflow boundaries, respectively, defined as

$\partial\Omega_{\rm
in}:=\\{\bm{X}\in\partial\Omega:\bm{u}\cdot\mathbf{n}<0\\}\quad\mbox{ and
}\quad\partial\Omega_{\rm
out}:=\\{\bm{X}\in\partial\Omega:\bm{u}\cdot\mathbf{n}\geq 0\\}.$
(Supplementary eq. 6)

The detail of discretization, similar to problems #1 and 2 that we use
discontinuous Galerkin approximation, as well as the FOM source codes, could
be found in [37, 31]. We note that as our finite element solver utilizes an
adaptive time-stepping [37, 31], each snapshot may have a different number of
time-steps $N^{t}$. We use geometries and meshes shown in Figs. 1a, c. Here,
we have one parameter $\bm{\mu}=(\mu_{1}=\mathrm{Ra})$, and
$\mathrm{Ra}=[350.0,450.0]$.

###  Supplementary sec. 2.4 Problem #4: Hyperelasticity problem with external
forces as a parameter

Problem #4 focuses on the finite deformation of a hyperelastic material, which
is applicable to many engineering applications, such as biomedical engineering
and material science [38, 39, 40]. The weak form of mechanical equilibrium
equations in the reference configuration, where the domain before deformation
is represented with $\Omega_{0}$ and after deformation with $\Omega_{i}$,
where $i$ represents different states of deformation (i.e.,
$\bm{\delta}\mathbf{u}=\mathbf{u}_{i}-\mathbf{u}_{0}$), is

$\int_{\Omega_{0}}\mathbf{P}:\nabla(\boldsymbol{\delta}\mathbf{u})\,dV-\int_{\Omega_{0}}\mathbf{B}\cdot\boldsymbol{\delta}\mathbf{u}\,dVs-\int_{\partial\Omega_{N}}\mathbf{T}\cdot\boldsymbol{\delta}\mathbf{u}\,dS=0,$
(Supplementary eq. 7)

where $\mathbf{P}$ is the 1st Piola Kirchhoff stress tensor, which is a
function of Young’s modulus ($\mathrm{E}$) and Poisson ratio ($\nu$).
$\mathbf{B}$ and $\mathbf{T}$ are the body and traction forces, respectively,
and $\mathbf{u}$ is a displacement vector (primary variable). To approximate
$\mathbf{u}$ (i.e., $\mathbf{u}_{h}$), following [32], we use a discontinuous
Galerkin approximation of the first order. We use PETSc SNES as a nonlinear
solver, and MUMPS as a linear solver [41] with absolute and relative
tolerances of $1\times 10^{-6}$ and $1\times 10^{-16}$, respectively. We
utilize a backtracking line search with slope descent parameter of $1\times
10^{-4}$, initial step length of $1.0$, and quadratic order of the
approximation. We utilize both 2- and 3-Dimensional domains for this problem.

####  Supplementary sec. 2.4.1 2-Dimensional domain

First, we use geometry and mesh shown in Figs. 1a-c. $\mathrm{E}=100$
$\mathrm{M}\mathrm{P}\mathrm{a}$, and $\nu=0.46$. We set
$\mathbf{B}=(0.0,0.0)$ and $\mathbf{T}=(0.0,0.0)$. We enforce
$\mathbf{u}=(\mu_{1},\mu_{2})$ at the boundary condition on the face of where
$y=1.0$, top surface. Other faces have a roller boundary condition (i.e., no
normal displacement). Our parameter space, here we have two parameters,
$\bm{\mu}=(0.05\mu_{1},0.05\mu_{2})$, and $\mu_{1}=[-1.0,1.0]$ and
$\mu_{2}=[-1.0,1.0]$.

####  Supplementary sec. 2.4.2 3-Dimensional domain

Here, we use a geometry and mesh shown in Supplementary fig. 1d,
$\mathrm{E}=10$ $\mathrm{P}\mathrm{a}$, and $\nu=0.3$. We set
$\mathbf{B}=(0.0,-0.5,0.0)$ and $\mathbf{T}=({\mu_{1}},0.0,0.0)$. We enforce
$\mathbf{u}=(0.0,0.0,0.0)$ at the boundary condition on the face of where
$x=0.0$, and

$\displaystyle\mathbf{u}=$ $\displaystyle(0.0,$ (Supplementary eq. 8)
$\displaystyle\mu_{2}(0.5+(y-0.5)\cos(\pi/3)-(z-0.5)\sin(\pi/3)-y)/2,$
$\displaystyle\mu_{2}(0.5+(y-0.5)\sin(\pi/3)+(z-0.5)\cos(\pi/3)-x))/2)$

at the boundary condition on the face of where $x=1.0$. Other faces have a
roller boundary condition (i.e., no normal displacement). Our parameter space
is similar to the 2-Dimensional case as we have two parameters,
$\bm{\mu}=(\mu_{1},\mu_{2})$, and $\mu_{1}=[0.1,0.9]$ and $\mu_{2}=[0.1,0.9]$.

##  Supplementary sec. 3 A collections of reduced order model used in this
manuscript

###  Supplementary sec. 3.1 Reduced order modeling with Barlow Twins self-
supervised learning

[24] has proposed a unified data-driven ROM, BT-ROM, that (1) bridges the
performance gap between the linear and nonlinear manifold approaches and (2)
can operate on unstructured meshes, which provides flexibility in its
application to standard numerical solvers, on-site measurements, or
experimental data. Even though it has been extended to boosting BT-ROM (BBT-
ROM) and BT-ROM with uncertainty quantification (UQ-BT-ROM) [32, 25],
throughout this study, we will focus only on the BT-ROM variation.

The summary of BT-ROM is shown in Supplementary fig. 2. Our BT-ROM starts with
an initialization of the training set $\bm{\mu}$, our parameter space.
Subsequently, we query a FOM (see Supplementary sec. 2) for each parameter
$\bm{\mu}$ in the training set. We note that the same procedures apply for
validation and testing sets, but we will only discuss the training set for
brevity. The third step entails a data compression stage through training BT-
AE developed in [24].

Supplementary fig. 2: The summary of procedures taken to establish the
proposed BT-AE

We illustrate our algorithm in Al. 1, but for further details and
implementation. Here, we only summarize the main components; please refer to
[24, 32]. The training involves two tasks; the first is the training of BT
(encoder and projector) - the outer loop with $\mathbf{B}_{\mathrm{outer}}$.
The second sub-task is the training of AE (encoder and decoder) - the inner
loop with $\mathbf{B}_{\mathrm{inner}}$. First, we alter our training set by
creating $\bm{X}_{h,A}\left(t,\bm{\mu}\right)$ and
$\bm{X}_{h,B}\left(t,\bm{\mu}\right)$ from $\bm{X}\left(t,\bm{\mu}\right)$)
through the addition of random noise

$\widetilde{\bm{X}_{h,A}}\left(t,\bm{\mu}\right),\widetilde{\bm{X}_{h,B}}\left(t,\bm{\mu}\right)=\bm{X}\left(t,\bm{\mu}\right)+\epsilon\operatorname{SD}\left(\bm{X}\left(t,\bm{\mu}\right)\right)\mathcal{G}\left(0,1\right)$
(Supplementary eq. 9)

where
$\widetilde{\bm{X}_{h,A}}\left(t,\bm{\mu}\right),\widetilde{\bm{X}_{h,B}}\left(t,\bm{\mu}\right)$
are altered $\bm{X}\left(t,\bm{\mu}\right)$ by (Supplementary eq. 9). The
constant $\epsilon$, which is set to $0.1$, determines the noise level, and
$\mathcal{G}\left(0,1\right)$ is a random value that is sampled from the
standard normal distribution with mean and standard deviation of zero and one,
respectively. Subsequently, we pass
$\widetilde{\bm{X}_{h,A}}\left(t,\bm{\mu}\right),\widetilde{\bm{X}_{h,B}}\left(t,\bm{\mu}\right)$
through Gaussian blur operation

${\bm{X}_{h,A}}\left(t,\bm{\mu}\right),{\bm{X}_{h,B}}\left(t,\bm{\mu}\right)=\frac{1}{\sqrt{2\pi\operatorname{SD}\left(\widetilde{\bm{X}_{h,A}}\left(t,\bm{\mu}\right),\widetilde{\bm{X}_{h,B}}\left(t,\bm{\mu}\right)\right)^{2}}}\mathrm{exp}\left({-\frac{{\widetilde{\bm{X}_{h,A}}\left(t,\bm{\mu}\right),\widetilde{\bm{X}_{h,B}}\left(t,\bm{\mu}\right)}^{2}}{2\operatorname{SD}\left(\widetilde{\bm{X}_{h,A}}\left(t,\bm{\mu}\right),\widetilde{\bm{X}_{h,B}}\left(t,\bm{\mu}\right)\right)^{2}}}\right)$
(Supplementary eq. 10)

to obtain $\bm{X}_{h,A}\left(t,\bm{\mu}\right)$ and
$\bm{X}_{h,B}\left(t,\bm{\mu}\right)$. Next, we pass
$\bm{X}_{h,A}\left(t,\bm{\mu}\right)$ and
$\bm{X}_{h,B}\left(t,\bm{\mu}\right)$ to the encoder (it is noted we have only
one encoder) resulting in $\bm{z}_{A}^{\bm{X}}\left(t,\bm{\mu}\right)$ and
$\bm{z}_{B}^{\bm{X}}\left(t,\bm{\mu}\right)$. We then use
$\bm{z}_{A}^{\bm{X}}\left(t,\bm{\mu}\right)$ and
$\bm{z}_{B}^{\bm{X}}\left(t,\bm{\mu}\right)$ as an input to the projector
resulting in the cross-correlation matrix
$\mathbf{C}^{\bm{X}}\left(t,\bm{\mu}\right)$.
$\mathbf{C}^{\bm{X}}\left(t,\bm{\mu}\right)$ is a square matrix with the
dimensionality of the projector’s output. The Barlow Twins loss
$\mathcal{L}_{\mathrm{BT}}^{\bm{X}}$, BT loss, is then calculated using

$\mathcal{L}_{\mathrm{BT}}^{\bm{X}}:=\mathcal{L}_{\mathrm{I}}^{\bm{X}}+\mathcal{L}_{\mathrm{RR}}^{\bm{X}}$
(Supplementary eq. 11)

where

$\mathcal{L}_{\mathrm{I}}^{\bm{X}}:=\sum_{i}\left(1-\mathbf{C}_{ii}^{\bm{X}}\left(t,\bm{\mu}\right)\right)^{2},$
(Supplementary eq. 12)

and

$\mathcal{L}_{\mathrm{RR}}^{\bm{X}}:=\lambda\sum_{i}\sum_{j\neq
i}{\mathbf{C}_{ij}^{\bm{X}}\left(t,\bm{\mu}\right)}^{2}.$ (Supplementary eq.
13)

Here, $\mathbf{C}_{ii}^{\bm{X}}\left(t,\bm{\mu}\right)$ denotes the $i$-th
diagonal entry of $\mathbf{C}^{\bm{X}}\left(t,\bm{\mu}\right)$, $\lambda$ is
set to $5\times 10^{-3}$, and $\mathbf{C}_{ij}^{\bm{X}}$ are off-diagonal
entries of $\mathbf{C}^{\bm{X}}$.

Subsequently, we obtain $\bm{z}^{\bm{X}}\left(t,\bm{\mu}\right)$ by passing
$\bm{X}_{h}\left(t,\bm{\mu}\right)$ to the encoder. The
$\bm{z}^{\bm{X}}\left(t,\bm{\mu}\right)$ is used to reconstruct
$\widehat{\bm{X}}_{h}\left(t,\bm{\mu}\right)$ through the decoder, and we
calculate our data compression loss or AE loss
($\mathcal{L}_{\mathrm{AE}}^{\bm{X}}$) using

${\mathcal{L}_{\mathrm{AE}}^{\bm{X}}}:={\mathrm{MSE}^{\bm{X}}}=\frac{1}{\mathrm{M}N^{t}}\sum_{i=1}^{\mathrm{M}}\sum_{k=0}^{N^{t}}\left|\widehat{\bm{X}}_{h}\left(t^{k},\bm{\mu}^{(i)}\right)-\bm{X}_{h}\left(t^{k},\bm{\mu}^{(i)}\right)\right|^{2}.$
(Supplementary eq. 14)

We use the adaptive moment estimation (ADAM) algorithm [42] to train the
framework. The learning rate ($\eta$) is calculated as [43]

$\eta_{c}=\eta_{\min}+\frac{1}{2}\left(\eta_{\max}-\eta_{\min}\right)\left(1+\cos\left(\frac{\mathrm{step_{c}}}{\mathrm{step_{f}}}\pi\right)\right)$
(Supplementary eq. 15)

where $\eta_{c}$ is a learning rate at step $\mathrm{step_{c}}$, $\eta_{\min}$
is the minimum learning rate, which is set as $1\times 10^{-16}$,
$\eta_{\max}$ is the initial learning rate, which is selected as $1\times
10^{-5}$, $\mathrm{step_{c}}$ is the current step, and $\mathrm{step_{f}}$ is
the final step. To prevent our networks from overfitting behavior, we follow
early stopping and generalized cross-validation criteria [44, 45]. Note that
instead of literally stopping our training cycle, we only save the set of
trained weights and biases from being used in the online phase when the
current validation loss is lower than the lowest validation from all the
previous training cycles.

Algorithm 1 Training autoencoder (AE) with Barlow Twins (BT) self-supervised
learning (BT-ROM)

# Training data $\bm{X}_{h}$ and distorted data $\bm{X}_{h,A}$, $\bm{X}_{h,B}$
are input of encoder
# latent spaces $\bm{z}^{\bm{X}}$, $\bm{z}_{A}^{\bm{X}}$, and
$\bm{z}_{B}^{\bm{X}}$ are output encoder
# latent space $\bm{z}^{\bm{X}}$ is output of decoder
# Approximation of $\bm{X}_{h}$, i.e., $\widehat{\bm{X}_{h}}$ is output of
decoder
# latent spaces $\bm{z}_{A}^{\bm{X}}$ and $\bm{z}_{B}^{\bm{X}}$ are input of
projector
# cross-correlation matrix $\mathbf{C}^{\bm{X}}$ is output of projector

1:Initialize (or load pre-trained models) encoder, decoder, and projector
$\triangleright$ size of latent space $\mathrm{Q}$ has to be specified.

2:Initialize (or load pre-trained optimizers) three optimizers for each of
encoder, decoder, and projector

3:Load training set $\bm{\mu}$ $\triangleright$ the total training data is
$\mathrm{M}N^{t}$

4:Randomly select 5% of $\mathrm{M}N^{t}$ for validation set
$\bm{\mu}_{\mathrm{validation}}$ $\triangleright$ the total training data
becomes 95% of $\mathrm{M}N^{t}$

5:Add random noise $\triangleright$ see (Supplementary eq. 9)

6:Add Gaussian blur $\triangleright$ see (Supplementary eq. 10)

7:From step 5 and 6, we obtain $\bm{X}_{h,A}\left(t,\bm{\mu}\right)$ and
$\bm{X}_{h,B}\left(t,\bm{\mu}\right)$ from $\bm{X}\left(t,\bm{\mu}\right)$

8:for each epoch do

9:  _Outer loop: training BT_ $\triangleright$ Batch size
$\mathbf{B}_{\mathrm{outer}}$

10: for each $\mathbf{B}_{\mathrm{outer}}$ do

11:
$\bm{z}^{\bm{X}}_{A}\left(t,\bm{\mu}\right)=\mathrm{encoder}\left({\bm{X}}_{h,A}\left(t,\bm{\mu}\right)\right)$

12:
$\bm{z}^{\bm{X}}_{B}\left(t,\bm{\mu}\right)=\mathrm{encoder}\left({\bm{X}}_{h,B}\left(t,\bm{\mu}\right)\right)$

13:
$\mathbf{C}^{\bm{X}}\left(t,\bm{\mu}\right)=\mathrm{projector}\left(\bm{z}^{\bm{X}}_{A}\left(t,\bm{\mu}\right),\bm{z}^{\bm{X}}_{B}\left(t,\bm{\mu}\right)\right)$

14: Calculate BT loss $\mathcal{L}_{\mathrm{BT}}^{\bm{X}}$ $\triangleright$
see (Supplementary eq. 11)

15: Back-propagation of BT loss w.r.t. each
$\mathrm{encoder}\left({\mathbf{W}},{\mathbf{b}}\right)$ and
$\mathrm{projector}\left({\mathbf{W}},{\mathbf{b}}\right)$

16: Update $\mathrm{encoder}\left({\mathbf{W}},{\mathbf{b}}\right)$ and
$\mathrm{projector}\left({\mathbf{W}},{\mathbf{b}}\right)$ using BT optimizer

17: Update learning rate $\eta_{c}$ of BT optimizer $\triangleright$ see
(Supplementary eq. 15)

18:  _Inner loop: training AE_ $\triangleright$ Batch size
$\mathbf{B}_{\mathrm{inner}}$

19: for each $\mathbf{B}_{\mathrm{inner}}$ do

20:
$\bm{z}^{\bm{X}}\left(t,\bm{\mu}\right)=\mathrm{encoder}\left({\bm{X}}_{h}\left(t,\bm{\mu}\right)\right)$

21:
$\widehat{{\bm{X}}_{h}}\left(t,\bm{\mu}\right)=\mathrm{decoder}\left(\bm{z}^{\bm{X}}\left(t,\bm{\mu}\right)\right)$

22: Calculate AE loss ${\mathcal{L}_{\mathrm{AE}}^{\bm{X}}}$ (data compression
loss) $\triangleright$ see (Supplementary eq. 14)

23: Back-propagation of AE loss w.r.t. each
$\mathrm{encoder}\left({\mathbf{W}},{\mathbf{b}}\right)$ and
$\mathrm{decoder}\left({\mathbf{W}},{\mathbf{b}}\right)$

24: Update $\mathrm{encoder}\left({\mathbf{W}},{\mathbf{b}}\right)$ and
$\mathrm{decoder}\left({\mathbf{W}},{\mathbf{b}}\right)$ using AE optimizer

25: Update learning rate $\eta_{c}$ of AE optimizer $\triangleright$ see
(Supplementary eq. 15)

26: end for

27: end for

28:end for

Reflecting the third step in Supplementary fig. 2.

We reiterate the autoencoder’s architecture used for BT-ROM in the
Supplementary tab. 2; please refer to [24, 32] for further details on
architecture search and hyper-parameters tuning. To sum up, the BT-ROM model
is composed of one encoder, one decoder, and one projector with their own sets
of weight matrices ($\mathbf{W}$) and biases (${\mathbf{b}}$). Each linear
layer is in the Supplementary tab. 2 is subjected to the LeakyReLU activation
function with a negative slope of 0.2. No batch normalization or dropout
layers are used in this architecture, making it preferable and lightweight.

Supplementary tab. 2: Autoencoder - reflecting the third step in Supplementary fig. 2. (input and output sizes are represented by [$\mathrm{B}$, $\mathrm{DOF}$]. $\mathrm{B}$ is a batch size, and $\bm{z}$ is nonlinear manifolds.) block | input size | output size
---|---|---
$1^{\mathrm{st}}$ linear layer | [$\mathrm{B}$, $\mathrm{DOF}$] | [$\mathrm{B}$, $\mathrm{int(DOF}/2)$]
$2^{\mathrm{nd}}$ linear layer | [$\mathrm{B}$, $\mathrm{int(DOF}/2)$] | [$\mathrm{B}$, $\mathrm{int(DOF}/4)$]
$3^{\mathrm{rd}}$ linear layer | [$\mathrm{B}$, $\mathrm{int(DOF}/4)$] | [$\mathrm{B}$, $\mathrm{int(DOF}/8)$]
$4^{\mathrm{th}}$ linear layer | [$\mathrm{B}$, $\mathrm{int(DOF}/8)$] | [$\mathrm{B}$, $\mathrm{int(DOF}/16)$]
$5^{\mathrm{th}}$ linear layer | [$\mathrm{B}$, $\mathrm{int(DOF}/16)$] | [$\mathrm{B}$, $\mathrm{int(DOF}/32)$]
$1^{\mathrm{st}}$ bottleneck | reshape([$\mathrm{B}$, $\mathrm{int(DOF}/32)$]) | [$\mathrm{B}$, $\bm{z}$ ]
$2^{\mathrm{nd}}$ bottleneck | [$\mathrm{B}$, $\bm{z}$ ] | reshape([$\mathrm{B}$, $\mathrm{int(DOF}/32)$])
$6^{\mathrm{th}}$ linear layer | [$\mathrm{B}$, $\mathrm{int(DOF}/32)$] | [$\mathrm{B}$, $\mathrm{int(DOF}/16)$]
$7^{\mathrm{th}}$ linear layer | [$\mathrm{B}$, $\mathrm{int(DOF}/16)$] | [$\mathrm{B}$, $\mathrm{int(DOF}/8)$]
$8^{\mathrm{th}}$ linear layer | [$\mathrm{B}$, $\mathrm{int(DOF}/8)$] | [$\mathrm{B}$, $\mathrm{int(DOF}/4)$]
$9^{\mathrm{th}}$ linear layer | [$\mathrm{B}$, $\mathrm{int(DOF}/4)$] | [$\mathrm{B}$, $\mathrm{int(DOF}/2)$]
$10^{\mathrm{th}}$ linear layer | [$\mathrm{B}$, $\mathrm{int(DOF}/2)$] | [$\mathrm{B}$, $\mathrm{DOF}$]

For the next step, we map $t$ and $\bm{\mu}$ to its representation in the
linear and nonlinear manifold $\bm{z}^{\bm{X}}\left(t,\bm{\mu}\right)$. Note
that if we deal with steady-state problems, we simply drop the $t$ term. We
follow a procedure proposed by [24] and use linear radial basis function (RBF)
interpolation [46] to map $t$ and $\bm{\mu}$ to
$\bm{z}^{\bm{X}}\left(t,\bm{\mu}\right)$. During the online phase, we utilize
the _trained_ RBF and the _trained_ decoder to approximate
$\widehat{{\bm{X}}}_{h}\left(\cdot;t,\bm{\mu}\right)$ for each inquiry (i.e.,
a value of $\bm{\mu}$) through

$\widehat{\bm{z}}^{\bm{X}}\left(t,\bm{\mu}\right)=\operatorname{RBF}\left(t,\bm{\mu}\right),$
(Supplementary eq. 16)

and, subsequently,

$\widehat{{\bm{X}}}_{h}\left(t,\bm{\mu}\right)=\operatorname{decoder}\left(\widehat{\bm{z}}^{\bm{X}}\left(t,\bm{\mu}\right)\right).$
(Supplementary eq. 17)

##  Supplementary sec. 4 Data generation

Here, we briefly explain how we create our training, validation, and testing
sets. We divide our investigation into two phases. In the first phase, we
focus on different physics but similar topology (Sec. Similar topology in the
main text). Here, we illustrate the effects of the number of training samples
and the number of parents on the BT-ROM’s and p-BT-ROM’s performance. In the
second phase, we use different physics and topologies (Sec. Different
topologies in the main text). For this phase, we only focus on the effect of
the number of parents on the BT-ROM’s and p-BT-ROM’s performance. We have
described each topology in detail in Supplementary sec. 1.

In terms of physics problems utilized throughout this study, we have four main
problems. We have described each problem in detail in Supplementary sec. 2.
The first one (Problem #1) is a transport problem with given different
velocity fields, i.e., $\bm{q}$ is the superficial velocity vector defined by
$\bm{q}=(20.0,\mu_{1}\cos(\pi x))$, and our parameter
$\bm{\mu}=(\mu_{1}=[5.0,25.0])$ Supplementary sec. 2.1). For this problem, we
use a mesh provided in Supplementary fig. 1a, and a sample of this problem is
shown in Supplementary fig. 3a. The second problem (Problem #2) is a transport
problem with given different diffusivity coefficients, i.e.,
$\bm{\mu}=(\bm{D}=(\mu_{1},0.0,0.0,\mu_{1}))$ and $\mu_{1}=[0.1,1.0]$
(Supplementary sec. 2.2). For this problem, we use two meshes presented in
Figs. 1a, b, and their samples are illustrated in Supplementary fig. 3b and
Supplementary fig. 4a, respectively.

The third problem (Problem #3) is a gravity-driven flow in porous media with
the Rayleigh number ($\mathrm{Ra}$) as a parameter
$\bm{\mu}=(\mu_{1}=\mathrm{Ra})$, and $\mathrm{Ra}=[350.0,450.0]$
(Supplementary sec. 2.3). For this problem, we use two meshes presented in
Figs. 1a, c, and their samples are illustrated in Supplementary fig. 3c and
Supplementary fig. 4b, respectively. The fourth problem (Problem #4) is the
finite deformation of a hyperelastic material (Supplementary sec. 2.4). For
this problem, $\bm{\mu}$ is different for the problem in a 2-Dimensional
domain, mesh: Supplementary fig. 1a and sample: Supplementary fig. 3d, and a
3-Dimensional domain, mesh: Supplementary fig. 1d and sample: Supplementary
fig. 4c. For the 2-Dimensional setting, we enforce
$\mathbf{u}=(\mu_{1},\mu_{2})$ at the boundary condition on the face of where
$y=1.0$, top surface, and $\bm{\mu}=(0.05\mu_{1},0.05\mu_{2})$, and
$\mu_{1}=[-1.0,1.0]$ and $\mu_{2}=[-1.0,1.0]$. For the 3-Dimensional setting,
we have traction forces $\mathbf{T}=({\mu_{1}},0.0,0.0)$, and we enforce

$\displaystyle\mathbf{u}=$ $\displaystyle(0.0,$ (Supplementary eq. 18)
$\displaystyle\mu_{2}(0.5+(y-0.5)\cos(\pi/3)-(z-0.5)\sin(\pi/3)-y)/2,$
$\displaystyle\mu_{2}(0.5+(y-0.5)\sin(\pi/3)+(z-0.5)\cos(\pi/3)-x))/2)$

at the boundary condition on the face of where $x=1.0$.
$\bm{\mu}=(\mu_{1},\mu_{2})$, and $\mu_{1}=[0.1,0.9]$ and $\mu_{2}=[0.1,0.9]$.

Supplementary fig. 3: Samples of physics problems simulated using mesh
provided in Supplementary fig. 1a (similar topology): (a) Problem #1 (
Supplementary sec. 2.1), (b) Problem #2 (Supplementary sec. 2.2), (c) Problem
#3 (Supplementary sec. 2.3), and (d) Problem #4 (Supplementary sec. 2.4). We
note that problems #1 to #3 are time-dependent problems while Problem #4 is at
a steady-state. For (a) - (c), we show the last time-step. Supplementary fig.
4: Samples of physics problems simulated using different meshes provided in
Supplementary fig. 1b-d (different topologies): (a) Problem #2 (Supplementary
sec. 2.2) with mesh shown in Supplementary fig. 1b, (b) Problem #3
(Supplementary sec. 2.3) with mesh shown in Supplementary fig. 1c, and (c)
Problem #3 (Supplementary sec. 2.4) with mesh shown in Supplementary fig. 1d.
We note that problems #2 and #3 are time-dependent problems while Problem #4
is at a steady-state. For (a) - (b), we show the last time-step.

We generate training, validation, and test sets to train and evaluate our
model. For the training set, we uniformly sample our parameter space
$\bm{\mu}$ with $\mathrm{M}$ snapshots (i.e., assuming $\bm{\mu}=[\mu_{1}]$,
$\mu_{1}=(0.0,1.0)$, and $\mathrm{M}=3$, our $\mu_{1}=[0.0,0.5,1.0]$). For
Problems #1 to #3, our actual training set would be $\mathrm{M}N_{t}$ (i.e.,
$\mathrm{M}$ snapshots $\times$ number of time-steps for each $\mathrm{M}$).
We design our framework this way to handle data provided from FOM with
adaptive time-stepping, which is essential for an advection-dominated problem
(i.e., to satisfy Courant–Friedrichs–Lewy (CFL) condition, Problems #1 to #3).
The framework can also deliver the solution at any time, including times that
do not exist in the training set. To elaborate, this does not imply that our
model can extrapolate, but it can deliver any timestamps inside the range of
the training set (see [24, 31] for more details on this). For Problem #4, the
actual training set is $\mathrm{M}$ because it is a steady-state problem. We
then randomly select 5% of the training samples to be used as a validation
set. For Problems #1 to #3, the final total training samples are
$0.95\mathrm{M}N^{t}$ while, for Problem #4, the final total training samples
are $0.95\mathrm{M}$ For the testing set, we randomly sample our parameter
space $\bm{\mu}$ with $\mathrm{M_{test}}$ snapshots.

##  Supplementary sec. 5 Complementary information on numerical experiments

We present supplement information on numerical experiments throughout the
following sections. These results are used to support numerical examples in
the main text, Sec. Results. Throughout this section, we use The mean squared
error (MSE)

${\mathrm{MSE}}=\frac{1}{\mathrm{M}N^{t}}\sum_{i=1}^{\mathrm{M}}\sum_{j=0}^{N^{t}}\left|\widehat{\bm{X}}_{h}\left(t^{j},\bm{\mu}^{(i)}\right)-\bm{X}_{h}\left(t^{j},\bm{\mu}^{(i)}\right)\right|^{2},$
(Supplementary eq. 19)

and mean absolute error (MAE)

$\mathrm{MAE}=\frac{1}{\mathrm{M}N^{t}}\sum_{i=1}^{\mathrm{M}}\sum_{j=0}^{N^{t}}|\widehat{\bm{X}}_{h}\left(t^{j},\bm{\mu}^{(i)}\right)-\bm{X}_{h}\left(t^{j},\bm{\mu}^{(i)}\right)|.$
(Supplementary eq. 20)

as our evaluation matrices.

###  Supplementary sec. 5.1 Using trained models for an initialization

First, we use trained models’ weights and biases to initialize a new model. If
we use more than one trained model, we simply take an arithmetic average of
each layer for the new model. For illustration purposes, throughout this
section, we use only one trained model - Problem #1 (Supplementary sec. 2.1) -
for the initialization of the new model. Our results are presented in
Supplementary fig. 5 for the reconstruction loss of the validation set
((Supplementary eq. 14)) and Supplementary fig. 6 for the Barlow Twins loss
((Supplementary eq. 11)) of the validation set. To emphasize this process, we
pick the weights and biases of the trained model that delivers the least
validation loss, shown in a red circle in Figs. 5a and 6a, for the
initialization. In short, from these figures, one can observe that using the
trained model as a starting point, the BT-ROM delivers lower reconstruction
validations loss for Problem #2 to #4 ( Supplementary sec. 2.1 to
Supplementary sec. 2.4). However, this trend is not always true for the Barlow
Twins loss, as using the trained model as a starting point can cause an
adverse effect, especially for Problems #3 and #4.

Supplementary fig. 5: Using a trained model for an initialization: validation
reconstruction loss ((Supplementary eq. 14)) results of (a) physical Problem
#1 (Supplementary sec. 2.1) used as an initialization for the following
physical problems, (b) physical Problem #2 (Supplementary sec. 2.2) for both
large dataset ($\mathrm{M}N_{t}=19240$) and small dataset
($\mathrm{M}N_{t}=2405$), (c) physical Problem #3 (Supplementary sec. 2.3) for
both large dataset ($\mathrm{M}N_{t}=28966$) and small dataset
($\mathrm{M}N_{t}=3623$), and (d) physical Problem #4 (Supplementary sec. 2.4)
for both large dataset ($\mathrm{M}=1600$) and small dataset
($\mathrm{M}=100$). We note that problems #1 to #3 are time-dependent problems
while Problem #4 is at a steady-state. The red circle in (a) represents the
epoch that has the lowest validation loss; and subsequently, be used as an
initialization for the other problems. Supplementary fig. 6: Using a trained
model for an initialization: validation Barlow Twins loss ((Supplementary eq.
11)) results of (a) physical Problem #1 (Supplementary sec. 2.1) used as an
initialization for the following physical problems, (b) physical Problem #2
(Supplementary sec. 2.2) for both large dataset ($\mathrm{M}N_{t}=19240$) and
small dataset ($\mathrm{M}N_{t}=2405$), (c) physical Problem #3 (Supplementary
sec. 2.3) for both large dataset ($\mathrm{M}N_{t}=28966$) and small dataset
($\mathrm{M}N_{t}=3623$), and (d) physical Problem #4 - 2-Dimensional domain
(Supplementary sec. 2.4) for both large dataset ($\mathrm{M}=1600$) and small
dataset ($\mathrm{M}=100$). We note that problems #1 to #3 are time-dependent
problems while Problem #4 is at a steady-state. The red circle in (a)
represents the epoch that has the lowest validation loss; and subsequently, be
used as an initialization for the other problems.

###  Supplementary sec. 5.2 Progressive Barlow Twins reduced order modeling
without initialization

Next, we investigate the effect of using p-BT-ROM (see Sec. Methodology in the
main text and Supplementary sec. 3.1) on the validation losses: the
reconstruction loss of the validation set ((Supplementary eq. 14)) and the
Barlow Twins loss ((Supplementary eq. 11)) of the validation set. We emphasize
that here we do not use the trained models’ (or parents’) weights and biases
to initialize our p-BT-ROM. The results of reconstruction (Supplementary fig.
7) and Barlow Twins (Supplementary fig. 8)losses are in line; so, we will only
discuss the reconstruction loss results. The schematic of p-BT-ROM specifies
child$\sim$parent(s) relationship is shown in Supplementary fig. 7a. To
elaborate, for Supplementary fig. 7b, we test our model using Problem #2
(Supplementary sec. 2.2), and our p-BT-ROM has only one parent, Problem #1.
One can see that p-BT-ROM delivers more accurate results.

For Supplementary fig. 7c, we test our model using Problem #3 (Supplementary
sec. 2.3), and the p-BT-ROM has two parents, Problem #1 and Problem #2. Again,
the validation loss of p-BT-ROM is much lower than that of BT-ROM (no parent).
Lastly, we use Problem #4 (Supplementary sec. 2.4) to test our model, and the
results align with previous observations. Hence, from these experiments, we
observe that using p-BT-ROM can improve the model accuracy substantially.

Supplementary fig. 7: p-BT-ROM without initialization: (a) schematic of p-BT-
ROM specifies child$\sim$parent(s) relationship, validation reconstruction
loss ((Supplementary eq. 14)) results of (b) physical Problem #2
(Supplementary sec. 2.2) for both large dataset ($\mathrm{M}N_{t}=19240$) and
small dataset ($\mathrm{M}N_{t}=2405$), (c) physical Problem #3 (Supplementary
sec. 2.3) for both large dataset ($\mathrm{M}N_{t}=28966$) and small dataset
($\mathrm{M}N_{t}=3623$), and (d) physical Problem #4 - 2-Dimensional domain
(Supplementary sec. 2.4) for both large dataset ($\mathrm{M}=1600$) and small
dataset ($\mathrm{M}=100$). We note that problems #1 to #3 are time-dependent
problems while Problem #4 is at a steady-state. Supplementary fig. 8: p-BT-ROM
without initialization: (a) schematic of p-BT-ROM specifies
child$\sim$parent(s) relationship, validation Barlow Twins loss
((Supplementary eq. 11)) results of (b) physical Problem #2 (Supplementary
sec. 2.2) for both large dataset ($\mathrm{M}N_{t}=19240$) and small dataset
($\mathrm{M}N_{t}=2405$), (c) physical Problem #3 (Supplementary sec. 2.3) for
both large dataset ($\mathrm{M}N_{t}=28966$) and small dataset
($\mathrm{M}N_{t}=3623$), and (d) physical Problem #4 - 2-Dimensional domain
(Supplementary sec. 2.4) for both large dataset ($\mathrm{M}=1600$) and small
dataset ($\mathrm{M}=100$). We note that problems #1 to #3 are time-dependent
problems while Problem #4 is at a steady-state.

###  Supplementary sec. 5.3 Progressive Barlow Twins reduced order modeling
with initialization

As previously illustrated that using a trained model as an initialization
(Supplementary sec. 5.1) as well as a progressive reduced order model
(Supplementary sec. 5.2) generally increases our BT-ROM accuracy. Throughout
this section, we investigate the p-BT-ROM with initialization (p-BT-ROM with
init.) effects on the model’s accuracy (i.e., investigate the effect of using
p-BT-ROM (see Sec. Methodology in the main text and Supplementary sec. 3.1) on
the validation losses: the reconstruction loss ((Supplementary eq. 14)) of the
validation set and the Barlow Twins loss ((Supplementary eq. 11)) of the
validation set). To reiterate, if we use more than one trained model (i.e.,
p-BT-ROM has more than one parent model), we take an arithmetic average of
each layer for all parents for the new model.

We present the reconstruction loss ((Supplementary eq. 14)) of the validation
set and the Barlow Twins loss ((Supplementary eq. 11)) of the validation set
in Figs. 9 and 10, respectively. The schematic of p-BT-ROM specifies
child$\sim$parent(s) relationship is shown in Supplementary fig. 9a. For Figs.
9b and 10b, we focus on Problem #2 (Supplementary sec. 2.2) by using Problem
#1 (Supplementary sec. 2.1) as a parent. We observe that for the
reconstruction loss using p-BT-ROM and p-BT-ROM with init. are quite similar;
however, p-BT-ROM with init. delivers the best Barlow Twins loss result.

For Supplementary fig. 9c and 10c, we focus on Problem #3 (Supplementary sec.
2.3) with Problem #1 and #2 as its parents. Here, the reconstruction and
Barlow Twins loss using p-BT-ROM and p-BT-ROM with init. are not much
different and are best. Lastly, For Supplementary fig. 9d and 10d, we focus on
Problem #4 (Supplementary sec. 2.4) with Problem #1, #2, and #3 as its
parents. For this problem, the p-BT-ROM with init. delivers the best result
for the reconstruction loss, while the p-BT-ROM delivers the best result for
the Barlow Twins loss.

Supplementary fig. 9: p-BT-ROM with initialization: (a) schematic of p-BT-ROM
specifies child$\sim$parent(s) relationship and how to initialize, validation
reconstruction loss ((Supplementary eq. 14)) results of (b) physical Problem
#2 (Supplementary sec. 2.2) for both large dataset ($\mathrm{M}N_{t}=19240$)
and small dataset ($\mathrm{M}N_{t}=2405$), (c) physical Problem #3
(Supplementary sec. 2.3) for both large dataset ($\mathrm{M}N_{t}=28966$) and
small dataset ($\mathrm{M}N_{t}=3623$), and (d) physical Problem #4 -
2-Dimensional domain (Supplementary sec. 2.4) for both large dataset
($\mathrm{M}=1600$) and small dataset ($\mathrm{M}=100$). We note that
problems #1 to #3 are time-dependent problems while Problem #4 is at a steady-
state. Supplementary fig. 10: p-BT-ROM with initialization: (a) schematic of
p-BT-ROM specifies child$\sim$parent(s) relationship and how to initialize,
validation Barlow Twins loss ((Supplementary eq. 11)) results of (b) physical
Problem #2 (Supplementary sec. 2.2) for both large dataset
($\mathrm{M}N_{t}=19240$) and small dataset ($\mathrm{M}N_{t}=2405$), (c)
physical Problem #3 (Supplementary sec. 2.3) for both large dataset
($\mathrm{M}N_{t}=28966$) and small dataset ($\mathrm{M}N_{t}=3623$), and (d)
physical Problem #4 - 2-Dimensional domain (Supplementary sec. 2.4) for both
large dataset ($\mathrm{M}=1600$) and small dataset ($\mathrm{M}=100$). We
note that problems #1 to #3 are time-dependent problems while Problem #4 is at
a steady-state.

###  Supplementary sec. 5.4 Discussion on supplement information on numerical
experiments

From Supplementary sec. 5.1, Supplementary sec. 5.2, and Supplementary sec.
5.3, we can see that by using a trained model as an initialization (init.),
p-BT-ROM, and p-BT-ROM with init. can improve our BT-ROM’s accuracy. However,
p-BT-ROM and p-BT-ROM with init. seem to provide a better result compared to
using init. alone. However, there is no clear winner between these two
approaches. Hence, for all numerical examples shown in Sec. Results in the
main text, we use p-BT-ROM with init. We note that the computation cost of
init. is negligible compared to the training p-BT-ROM itself.

###  Supplementary sec. 5.5 Supplement information of similar topology section
(Sec. Similar topology in the main text)

Here, we investigate p-BT-ROM’s performance using one mesh, see Supplementary
fig. 1a, for physics problems, Problems #1 to #4 (test model performance on
different physics but similar mesh or topology).

####  Supplementary sec. 5.5.1 Gravity-driven in porous media - Problem #3

Throughout this section, we focus on a gravity-driven in porous media problem
- Problem #3 (Supplementary sec. 2.3). One sample of this problem is shown in
Supplementary fig. 3c. We aim to investigate the model’s performance on a test
set ($\mathrm{M_{test}}N_{t}=7260$, $\mathrm{M_{test}}=10$) with different
training set sizes and a number of parents. This problem has three primary
variables, fluid pressure, fluid velocity, and fluid temperature. Since our
BT-ROM is non-intrusive, we only pick the fluid temperature as our quantity of
interest. The MSE ((Supplementary eq. 19)) and MAE ((Supplementary eq. 20))
results are presented in Supplementary fig. 11. The schematic of p-BT-ROM
specifies child$\sim$parent(s) relationship is shown in Supplementary fig.
12a. For the model with 1 parent, we use Problem #1 (Supplementary sec. 2.1)
as a parent with a training set of $\mathrm{M}N_{t}=2405$ ($\mathrm{M}=5$,
i.e., small dataset in Supplementary sec. 5). For the model with 2 parents,
Problems #1 and #2 (Supplementary sec. 2.2) with both have a training set of
$\mathrm{M}N_{t}=2405$ ($\mathrm{M}=5$, i.e., small dataset in Supplementary
sec. 5).

From Supplementary fig. 11, we observe that the model with no parent and has a
small training set ($\mathrm{M}N_{t}=3623$, $\mathrm{M}=5$), shown in red,
performs the worst. As we add more parents, shown in blue and green, the
models gain their capability and deliver better accuracy. However, even the
model with two parents, shown in blue, still could not achieve the same level
of accuracy as the model that trained with a large training set
($\mathrm{M}N_{t}=28966$, $\mathrm{M}=40$), but has zero parent.

Supplementary fig. 11: Similar topology - Problem #3 (Supplementary sec. 2.3):
(a) mean squared error and (b) mean absolute error results for large dataset
($\mathrm{M}N_{t}=28966$, $\mathrm{M}=40$) and small dataset
($\mathrm{M}N_{t}=3623,\mathrm{M}=5$) with a different number of parents. We
note that schematic of p-BT-ROM specifies child$\sim$parent(s) relationship is
shown in Supplementary fig. 12a.

Figs. 12b-c illustrate MAE as as a function of parameter $\bm{\mu}$ for three
different cases; large dataset ($\mathrm{M}N_{t}=28966$, $\mathrm{M}=40$) with
0 parent, small dataset ($\mathrm{M}N_{t}=3623$, $\mathrm{M}=5$) with 0
parent, and small dataset ($\mathrm{M}N_{t}=3623$, $\mathrm{M}=5$) with 2
parents, respectively. One can observe that the pattern of these three figures
is similar (i.e., the shape), but the magnitude of the MAE is different. This
observation implies that the p-BT-ROM only alters the magnitude of the error,
not its pattern.

Supplementary fig. 12: Similar topology - Problem #3 (Supplementary sec. 2.3):
(a) schematic of p-BT-ROM specifies child$\sim$parent(s) relationship, mean
absolute error as a function of parameter $\bm{\mu}$ results for (b) large
dataset ($\mathrm{M}N_{t}=28966,\mathrm{M}=40$) with 0 parent, (c) small
dataset ($\mathrm{M}N_{t}=3623,\mathrm{M}=5$) with 0 parent, and (d) small
dataset ($\mathrm{M}N_{t}=3623,\mathrm{M}=5$) with 2 parents. We note that
problems #1 to #3 are time-dependent problems, and their description can be
found in ( Supplementary sec. 2.1, Supplementary sec. 2.2, and Supplementary
sec. 2.3).

####  Supplementary sec. 5.5.2 Hyperelasticity problem - Problem #4

Next, we move to a physics problem of finite deformation of hyperelastic
material - Problem #4 (Supplementary sec. 2.4). One sample of this problem is
shown in Supplementary fig. 3d. Here, we want to investigate the model’s
performance by using a solid mechanics at a steady-state solution (i.e., all
parents, Problems #1 to #3, are transient fluid mechanics problems). We use a
test set of $\mathrm{M_{test}}=100$. with different training set sizes and a
number of parents. The quantity of interest, here, is a magnitude of
displacement.

The MSE and MAE are shown in Supplementary fig. 13, and the schematic of p-BT-
ROM specifies child$\sim$parent(s) relationship is shown in Supplementary fig.
14a. It is similar to the previous section. For the model with 1 parent, we
use Problem #1 (Supplementary sec. 2.1) as a parent with a training set of
$\mathrm{M}N_{t}=2405$ ($\mathrm{M}=5$), for the model with 2 parents,
Problems #1 and #2 (Supplementary sec. 2.2) with both have a training set of
$\mathrm{M}N_{t}=2405$ ($\mathrm{M}=5$), for the model with 3 parents, we add
Problem #3 (Supplementary sec. 2.3) with a training set of
$\mathrm{M}N_{t}=3623$ ($\mathrm{M}=5$). We note that all parents are trained
using a small dataset (more details can be found in Supplementary sec. 5).

From Supplementary fig. 13, again, we observe that the model with no parent
and has a small training set ($\mathrm{M}=100$), shown in red, performs the
worst, while the model with a large training set ($\mathrm{M}=1600$) delivers
the best accuracy. However, with a small training set, as we add more parents,
shown in blue, green, and purple, the models gain their capability and deliver
better accuracy. However, we speculate that it would be challenging to reach
the same level of accuracy as that of the large dataset, even if we add even
more parents.

Supplementary fig. 13: Similar topology - Problem #4 (Supplementary sec. 2.4):
(a) mean squared error and (b) mean absolute error results for large dataset
($\mathrm{M}=1600$) and small dataset ($\mathrm{M}=100$) with a different
number of parents. We note that schematic of p-BT-ROM specifies
child$\sim$parent(s) relationship is shown in Supplementary fig. 14a.

Similar to the previous section as we use Figs. 14b-c to present the MAE
results as a function of parameter $\bm{\mu}$ for three different cases; large
dataset ($\mathrm{M}=1600$) with 0 parent, small dataset ($\mathrm{M}=100$)
with 0 parent, and small dataset ($\mathrm{M}=100$) with 3 parents,
respectively. Again, one can observe that the pattern of these three figures
is similar (i.e., the shape), but the magnitude of the MAE is different. This
observation implies that the p-BT-ROM only alters the magnitude of the error,
not its pattern.

Supplementary fig. 14: Similar topology - Problem #4 (Supplementary sec. 2.4):
(a) schematic of p-BT-ROM specifies child$\sim$parent(s) relationship, mean
absolute error as a function of parameter $\bm{\mu}$ results for (b) large
dataset ($\mathrm{M}=1600$) with 0 parent, (c) small dataset
($\mathrm{M}=100$) with 0 parent, and (d) small dataset ($\mathrm{M}=100$)
with 3 parents. We note that problems #1 to #3 are time-dependent problems
while Problem #4 is at a steady-state, and their description can be found in (
Supplementary sec. 2.1, Supplementary sec. 2.2, Supplementary sec. 2.3, and
Supplementary sec. 2.4).

##  Supplementary sec. 6 Summary of number of parameters of each model

We present the number of parameters of the (p-)BT-ROM with different numbers
of parents in Supplementary tab. 3 for Supplementary fig. 2 and Supplementary
tab. 4 for Supplementary fig. 3. We observe that, as expected, as we increase
the number of parents, the more parameters the models contain (i.e., more
gates required to control the flow of information). We also note that even
though the number of parameters increases. The training and prediction costs
are not much different. We train and test our model using NVIDIA Quadro RTX
8000.

Supplementary tab. 3: Summary of the number of parameters of each model presented in Fig. 2 in the main text. parents | encoder | decoder | projector
---|---|---|---
0 | 47,744,112 | 47,752,562 | 11,168
1 | 59,647,413 | 95,500,636 | 21,728
2 | 71,550,714 | 143,248,710 | 32,288
3 | 83,454,015 | 190,996,784 | 42,848
Supplementary tab. 4: Summary of the number of parameters of each model presented in Fig. 3 in the main text. parents | encoder | decoder | projector
---|---|---|---
0 | 59,018,187 | 59,022,139 | 13,805
1 | 70,920,019 | 106,770,213 | 25,834
2 | 82,822,340 | 154,518,287 | 37,374
3 | 94,724,940 | 202,266,361 | 48,635
4 | 106,627,716 | 250,014,435 | 59,720

##  Supplementary sec. 7 Complementary information on numerical results for
different topologies section (Sec. Different topologies in the main text)

Supplementary fig. 15 is used to showing the pattern of MAE as a function of
parameter $\bm{\mu}$. We discuss this figure in Sec. Different topologies in
the main text.

Supplementary fig. 15: Different topologies - mean absolute error as a
function of parameter $\bm{\mu}$ results for (a) 0 and (b) 4 parents for
Problem #2 (Supplementary sec. 2.2) using the topology shown in Supplementary
fig. 1b, (c) 0 and (d) 4 parents for Problem #3 (Supplementary sec. 2.3) using
the topology shown in Supplementary fig. 1c, and (e) 0 and (f) 4 parents for
Problem #4 (Supplementary sec. 2.4) using the topology shown in Supplementary
fig. 1d.

## References

  * [1] Masood, E., Tollefson, J., Irwin, A. _et al._ Cop27 climate talks: what succeeded, what failed and what’s next.  _Nature_ 612, 16–17 (2022).
  * [2] Falk, J. _et al._ An urgent need for cop27: confronting converging crises.  _Sustainability Science_ 1–5 (2022).
  * [3] Hu, L. & Chugunova, T.  Multiple-point geostatistics for modeling subsurface heterogeneity: A comprehensive review.  _Water Resources Research_ 44 (2008).
  * [4] Hartmann, A.  Putting the cat in the box: why our models should consider subsurface heterogeneity at all scales.  _Wiley Interdisciplinary Reviews: Water_ 3, 478–486 (2016).
  * [5] Ginn, T. _et al._ Processes in microbial transport in the natural subsurface.  _Advances in Water Resources_ 25, 1017–1042 (2002).
  * [6] McCarthy, J. & McKay, L.  Colloid transport in the subsurface: Past, present, and future challenges.  _Vadose Zone Journal_ 3, 326–337 (2004).
  * [7] Evans, G., Blackledge, J. & Yardley, P.  _Numerical methods for partial differential equations_ (Springer Science & Business Media, 2012).
  * [8] Hesthaven, J., Rozza, G., Stamm, B. _et al._ _Certified reduced basis methods for parametrized partial differential equations_ (Springer, 2016).
  * [9] Chen, F. _et al._ Capacity assessment and cost analysis of geologic storage of hydrogen: A case study in intermountain-west region usa.  _International Journal of Hydrogen Energy_ 48, 9008–9022 (2023).
  * [10] Wen, G. _et al._ Real-time high-resolution co 2 geological storage prediction using nested fourier neural operators.  _Energy & Environmental Science_ 16, 1732–1741 (2023).
  * [11] Lengler, U., De Lucia, M. & Kühn, M.  The impact of heterogeneity on the distribution of co2: Numerical simulation of co2 storage at ketzin.  _International Journal of Greenhouse Gas Control_ 4, 1016–1025 (2010).
  * [12] Cho, Y. & Jun, H.  Estimation and uncertainty analysis of the co2 storage volume in the sleipner field via 4d reversible-jump markov-chain monte carlo.  _Journal of Petroleum Science and Engineering_ 200, 108333 (2021).
  * [13] Choi, Y., Coombs, D. & Anderson, R.  Sns: a solution-based nonlinear subspace method for time-dependent model order reduction.  _SIAM Journal on Scientific Computing_ 42, A1116–A1146 (2020).
  * [14] Kapteyn, M., Knezevic, D., Huynh, D., Tran, M. & Willcox, K.  Data-driven physics-based digital twins via a library of component-based reduced-order models.  _International Journal for Numerical Methods in Engineering_ 123, 2986–3003 (2022).
  * [15] Silva, V., Heaney, C., Li, Y. & Pain, C.  Data assimilation predictive gan (da-predgan) applied to a spatio-temporal compartmental model in epidemiology.  _Journal of Scientific Computing_ 94, 1–31 (2023).
  * [16] Qin, T., Chen, Z., Jakeman, J. & Xiu, D.  Deep learning of parameterized equations with applications to uncertainty quantification.  _International Journal for Uncertainty Quantification_ 11 (2021).
  * [17] Xu, K. & Darve, E.  Physics constrained learning for data-driven inverse modeling from sparse observations.  _Journal of Computational Physics_ 453, 110938 (2022).
  * [18] Pachalieva, A., O’Malley, D., Harp, D. & Viswanathan, H.  Physics-informed machine learning with differentiable programming for heterogeneous underground reservoir pressure management.  _Scientific Reports_ 12, 1–12 (2022).
  * [19] Raissi, M., Perdikaris, P. & Karniadakis, G.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ 378, 686–707 (2019).
  * [20] Weiss, K., Khoshgoftaar, T. & Wang, D.  A survey of transfer learning.  _Journal of Big data_ 3, 1–40 (2016).
  * [21] Zhuang, F. _et al._ A comprehensive survey on transfer learning.  _Proceedings of the IEEE_ 109, 43–76 (2020).
  * [22] Rusu, A. _et al._ Progressive neural networks.  _arXiv preprint arXiv:1606.04671_ (2016).
  * [23] Kadeethum, T. _et al._ A framework for data-driven solution and parameter estimation of pdes using conditional generative adversarial networks.  _Nature Computational Science_ 1, 819–829, DOI: https://doi.org/10.1038/s43588-021-00171-3 (2021).
  * [24] Kadeethum, T. _et al._ Reduced order modeling for flow and transport problems with barlow twins self-supervised learning.  _Scientific Reports_ 12, 1–18 (2022).
  * [25] Kadeethum, T., Jakeman, J., Choi, Y., Bouklas, N. & Yoon, H.  Epistemic uncertainty-aware barlow twins reduced order modeling for nonlinear contact problems.  _under review_ (2022).
  * [26] Goswami, S., Kontolati, K., Shields, M. & Karniadakis, G.  Deep transfer operator learning for partial differential equations under conditional shift.  _Nature Machine Intelligence_ 4, 1155–1164 (2022).
  * [27] Huang, Z., Dang, L., Xie, Y., Ma, W. & Chen, B.  Gate recurrent unit network based on hilbert-schmidt independence criterion for state-of-health estimation.  _arXiv preprint arXiv:2303.09497_ (2023).
  * [28] Vaswani, A. _et al._ Attention is all you need.  _Advances in neural information processing systems_ 30 (2017).
  * [29] Alesiani, F., Yu, S. & Yu, X.  Gated information bottleneck for generalization in sequential environments.  _Knowledge and Information Systems_ 65, 683–705 (2023).
  * [30] Drori, I. _et al._ A neural network solves, explains, and generates university math problems by program synthesis and few-shot learning at human level.  _Proceedings of the National Academy of Sciences_ 119, e2123433119 (2022).
  * [31] Kadeethum, T. _et al._ Non-intrusive reduced order modeling of natural convection in porous media using convolutional autoencoders: comparison with linear subspace techniques.  _Advances in Water Resources_ 104098 (2022).
  * [32] Kadeethum, T. _et al._ Enhancing high-fidelity nonlinear solver with reduced order model.  _Scientific Reports_ 12, 1–15 (2022).
  * [33] Lu, L., Jin, P., Pang, G., Zhang, Z. & Karniadakis, G.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  _Nature Machine Intelligence_ 3, 218–229 (2021).
  * [34] Choi, Y., Boncoraglio, G., Anderson, S., Amsallem, D. & Farhat, F.  Gradient-based constrained optimization using a database of linear reduced-order models.  _Journal of Computational Physics_ 423 (2020).
  * [35] Oommen, V., Shukla, K., Goswami, S., Dingreville, R. & Karniadakis, G.  Learning two-phase microstructure evolution using neural operators and autoencoder architectures.  _arXiv preprint arXiv:2204.07230_ (2022).
  * [36] Gorodetsky, A., Jakeman, J. & Geraci, G.  Mfnets: data efficient all-at-once learning of multifidelity surrogates as directed networks of information sources.  _Computational Mechanics_ 68, 741–758 (2021).
  * [37] Kadeethum, T., Lee, S., Ballarin, F., Choo, J. & Nick, H.  A locally conservative mixed finite element framework for coupled hydro-mechanical-chemical processes in heterogeneous porous media.  _Computers & Geosciences_ 104774 (2021).
  * [38] Kumar, N. & DasGupta, A.  On the contact problem of an inflated spherical hyperelastic membrane.  _International Journal of Non-Linear Mechanics_ 57, 130–139 (2013).
  * [39] Luo, Y., Li, M. & Kang, Z.  Topology optimization of hyperelastic structures with frictionless contact supports.  _International Journal of Solids and Structures_ 81, 373–382 (2016).
  * [40] Nezamabadi, S., Zahrouni, H. & Yvonnet, J.  Solving hyperelastic material problems by asymptotic numerical method.  _Computational mechanics_ 47, 77–92 (2011).
  * [41] Balay, S. _et al._ PETSc Users Manual.  Tech. Rep. ANL-95/11 - Revision 3.10, Argonne National Laboratory (2018).
  * [42] Kingma, D. & Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * [43] Loshchilov, I. & Hutter, F.  Sgdr: Stochastic gradient descent with warm restarts.  _arXiv preprint arXiv:1608.03983_ (2016).
  * [44] Prechelt, L.  Early stopping-but when?  In _Neural Networks: Tricks of the trade_ , 55–69 (Springer, 1998).
  * [45] Prechelt, L.  Automatic early stopping using cross validation: quantifying the criteria.  _Neural Networks_ 11, 761–767 (1998).
  * [46] Wright, G.  _Radial basis function interpolation: numerical and analytical developments_ (University of Colorado at Boulder, 2003).