Data-driven reduced order modeling of poroelasticity of heterogeneous media based on a discontinuous Galerkin approximation

Abstract

A simulation tool capable of speeding up the calculation for linear poroelasticity problems in heterogeneous porous media is of large practical interest for engineers, in particular, to effectively perform sensitivity analyses, uncertainty quantification, optimization, or control operations on the fluid pressure and bulk deformation fields. Towards this goal, we present here a non-intrusive model reduction framework using proper orthogonal decomposition (POD) and neural networks based on the usual offline-online paradigm. As the conductivity of porous media can be highly heterogeneous and span several orders of magnitude, we utilize the interior penalty discontinuous Galerkin (DG) method as a full order solver to handle discontinuity and ensure local mass conservation during the offline stage. We then use POD as a data compression tool and compare the nested POD technique, in which time and uncertain parameter domains are compressed consecutively, to the classical POD method in which all domains are compressed simultaneously. The neural networks are finally trained to map the set of uncertain parameters, which could correspond to material properties, boundary conditions, or geometric characteristics, to the collection of coefficients calculated from an $L^2$ projection over the reduced basis. We then perform a non-intrusive evaluation of the neural networks to obtain coefficients corresponding to new values of the uncertain parameters during the online stage. We show that our framework provides reasonable approximations of the DG solution, but it is significantly faster. Moreover, the reduced order framework can capture sharp discontinuities of both displacement and pressure fields resulting from the heterogeneity in the media conductivity, which is generally challenging for intrusive reduced order methods. The sources of error are presented, showing that the nested POD technique is computationally advantageous and still provides comparable accuracy to the classical POD method. We also explore the effect of different choices of the hyperparameters of the neural network on the framework performance.

Appendices

Appendix A: Finite element discretization and solution

We use the discontinuous Galerkin solver from Kadeethum et al. (2019), Kadeethum et al. (2019), Lee et al. (2019), Kadeethum et al. (2020), Kadeethum et al. (2021), and we briefly revisit the discretization in this section. We begin by introducing the necessary notation. Let ${\mathcal {T}}_h$ be a shape-regular triangulation obtained by a partition of $\varOmega$ into d-simplices (segments in $d=1$, triangles in $d=2$, tetrahedra in $d=3$). For each cell $T \in {\mathcal {T}}_h$, we denote by $h_{T}$ the diameter of T, and we set $h=\max _{T \in {\mathcal {T}}_h} h_{T}$ and $h_{l}=\min _{T \in {\mathcal {T}}_h} h_{T}$. We further denote by ${\mathcal {E}}_h$ the set of all facets (i.e., $d - 1$ dimensional entities connected to at least a $T \in {\mathcal {T}}_h$) and by ${\mathcal {E}}_h^{I}$ and ${\mathcal {E}}_h^{\partial }$ the collection of all interior and boundary facets, respectively. The boundary set ${\mathcal {E}}_h^{\partial }$ is decomposed into two disjoint subsets associated with the Dirichlet boundary facets, and the Neumann boundary facets for each of Eqs. (7) and (12). In particular, ${\mathcal {E}}_{h}^{D,u}$ and ${\mathcal {E}}_{h}^{N,u}$ correspond to the facets on $\partial \varOmega _u$ and $\partial \varOmega _{tr}$, respectively, for Eq. (7). On the other hand, for Eq. (12), ${\mathcal {E}}_{h}^{D,m}$ and ${\mathcal {E}}_{h}^{N,m}$ conform to $\partial \varOmega _p$ and $\partial \varOmega _{q}$, respectively.

We also define

$$\begin{aligned} e = \partial T^{+}\cap \partial T^{-}, \ \ e \in {\mathcal {E}}_h^I, \end{aligned}$$

where $T^{+}$ and $T^{-}$ are the two neighboring elements to e. We denote by $h_e$ the characteristic length of e calculated as

$$\begin{aligned} h_{e} :=\frac{{\text {meas}}\left( T^{+}\right) +{\text {meas}}\left( T^{-}\right) }{2 {\text {meas}}(e)}, \end{aligned}$$

(29)

depending on the argument, meas($\cdot$) represents the measure of a cell or of a facet.

Let ${\mathbf {n}}^{+}$ and ${\mathbf {n}}^{-}$ be the outward unit normal vectors to $\partial T^+$ and $\partial T^-$, respectively. For any given scalar function $\zeta : {\mathcal {T}}_h \rightarrow {\mathbb {R}}$ and vector function $\varvec{\tau }: {\mathcal {T}}_h \rightarrow {\mathbb {R}}^d$, we denote by $\zeta ^{\pm }$ and $\varvec{\tau }^{\pm }$ the restrictions of $\zeta$ and $\varvec{\tau }$ to $T^\pm$, respectively. Subsequently, we define the weighted average operator as

$$\begin{aligned} \{\zeta \}_{\delta e}=\delta _{e} \zeta ^{+}+\left( 1-\delta _{e}\right) \zeta ^{-}, \ \text { on } e \in {\mathcal {E}}_h^I, \end{aligned}$$

(30)

and

$$\begin{aligned} \{\varvec{\tau }\}_{\delta e}=\delta _{e} \varvec{\tau }^{+}+\left( 1-\delta _{e}\right) \varvec{\tau }^{-}, \ \text { on } e \in {\mathcal {E}}_h^I, \end{aligned}$$

(31)

where $\delta _{e}$ is calculated by Ern et al. (2009), Ern and Stephansen (2008):

$$\begin{aligned} \delta _{e} :=\frac{{k}^{-}_e}{{k}^{+}_e+{k}^{-}_e}. \end{aligned}$$

(32)

Here,

$$\begin{aligned} {k}^{+}_e :=\left( {\mathbf {n}}^{+}\right) ^{\intercal } \cdot \varvec{k}^{+} {\mathbf {n}}^{+}, \ \text { and } {k}^{-}_e :=\left( {\mathbf {n}}^{-}\right) ^{\intercal } \cdot \varvec{k}^{-} {\mathbf {n}}^{-}, \end{aligned}$$

(33)

where ${k_e}$ is a harmonic average of $k^{+}_e$ and ${k}^{-}_e$ which reads

$$\begin{aligned} {k_{e}}:= \frac{2{k}^{+}_e {k}^{-}_e}{{k}^{+}_e+{k}^{-}_e}, \end{aligned}$$

(34)

and $\varvec{k}$ is defined as in Eq. (11).

The jump across an interior edge will be defined as

$$\begin{aligned} \llbracket \zeta \rrbracket = \zeta ^+{\mathbf {n}}^++\zeta ^-{\mathbf {n}}^- \quad \text{ and } \quad \llbracket \tau \rrbracket = \varvec{\tau }^+\cdot {\mathbf {n}}^+ + \varvec{\tau }^-\cdot {\mathbf {n}}^- \quad \text{ on } e\in {\mathcal {E}}_h^I. \end{aligned}$$

Finally, for $e \in {\mathcal {E}}^{\partial }_h$, we set $\{\zeta \}_{\delta _e} := \zeta$ and $\{{\varvec{\tau }}\}_{\delta _e} := \varvec{\tau }$ for what concerns the definition of the weighted average operator, and $\llbracket \zeta \rrbracket := \zeta {\mathbf {n}}$ and $\llbracket \varvec{\tau } \rrbracket := \varvec{\tau } \cdot {\mathbf {n}}$ as definition of the jump operator.

1.1 Temporal discretization

We adapt the Biot’s system solver from Kadeethum et al. (2021), Kadeethum et al. (2020). The time domain, ${\mathbb {T}} = \left( 0, \mathrm {T}\right]$, is partitioned into $N^t$ open intervals such that, $0=: t^{0}<t^{1}<\cdots <t^{ N^t} := \mathrm {T}$. The length of the interval, $\varDelta t^n$, is defined as $\varDelta t^n=t^{n}-t^{n-1}$ where n represents the current time step. $\varDelta t^0$ is an initial $\varDelta t$, which is defined as $t^{1}-t^{0}$, while the other time steps, $\varDelta t^n$, are calculated as follows

$$\begin{aligned} \varDelta t^n := {\left\{ \begin{array}{ll} \varDelta t_{mult}\times \varDelta t^{n-1} &{} \text {if} \ \varDelta t^n \le \varDelta t_{max} \ \\ \varDelta t_{max} &{} \text {if} \ \varDelta t^n > \varDelta t_{max}, \end{array}\right. } \end{aligned}$$

(35)

where $\varDelta t_{mult}$ is a positive constant multiplier, and $\varDelta t_{max}$ is the maximum allowable time step. Then, let $\varphi (\cdot , t)$ be a scalar function and $\varphi ^{n}$ be its approximation at time $t^n$, i.e. $\varphi ^{n} \approx \varphi \left( t^{n}\right)$. We employ the following backward differentiation formula for time discretization of all primary variables (Ibrahim et al. 2007; Akinfenwa et al. 2013; Lee et al. 2018; Kadeethum et al. 2020)

$$\begin{aligned} \mathrm {BDF}_{1}\left( \varphi ^{n}\right) := \frac{1}{\varDelta t^n}\left( \varphi ^{n}-\varphi ^{n-1}\right) . \end{aligned}$$

(36)

1.2 Full discretization

Following (Liu et al. 2009; Kadeethum et al. 2021, 2020, 2019), in this study, the displacement field is approximated by the classical continuous Galerkin method (CG) method, and the fluid pressure field is discretized by discontinuous Galerkin (DG) method to ensure local mass conservation and provide a better flux approximation (Kadeethum et al. 2021, 2020).

We begin with defining the finite element space for the continuous Galerkin (CG) method for a vector-valued function

$$\begin{aligned} {\mathcal {U}}_{h}^{\mathrm {CG}_{k}}\left( {\mathcal {T}}_{h}\right) :=\left\{ \varvec{\psi }_{\varvec{u}} \in {\mathbb {C}}^{0}(\varOmega {; {\mathbb {R}}^d}) :\left. \varvec{\psi }_{\varvec{u}}\right| _{T} \in {\mathbb {Q}}_{k}(T{; {\mathbb {R}}^d}), \forall T \in {\mathcal {T}}_{h}\right\} , \end{aligned}$$

(37)

where k indicates the order of polynomial that can be approximated in this space, ${\mathbb {C}}^0(\varOmega {; {\mathbb {R}}^d})$ denotes the space of vector-valued piece-wise continuous polynomials, ${\mathbb {Q}}_{k}(T{; {\mathbb {R}}^d})$ is the space of polynomials of degree at most k over each element T, and ${\mathbb {R}}$ is a set of real numbers. We will denote in the following by $N_h^u$ the dimension of the space ${\mathcal {U}}_{h}^{\mathrm {CG}_{k}}\left( {\mathcal {T}}_{h}\right)$, i.e. the number of degrees of freedom for the displacement approximation.

Next, the DG space for scalar-valued functions is defined as

$$\begin{aligned} {\mathcal {P}}_{h}^{\mathrm {DG}_{k}}\left( {\mathcal {T}}_{h}\right) :=\left\{ \psi _p \in L^{2}(\varOmega ) :\left. \psi _p\right| _{T} \in {\mathbb {Q}}_{k}(T), \forall T \in {\mathcal {T}}_{h}\right\} , \end{aligned}$$

(38)

where $L^{2}(\varOmega )$ is the space of square integrable functions. This non conforming finite element space allows us to consider discontinuous coefficients and preserves the local mass conservation. We will denote in the following by $N_h^p$ the dimension of ${\mathcal {P}}_{h}^{\mathrm {DG}_{k}}\left( {\mathcal {T}}_{h}\right)$.

We seek the approximated displacement ($\varvec{u}_{h}$) and pressure (${p}_{h}$) solutions by discretizing the linear momentum balance equation Eq. (7) employing the above CG finite element spaces for $\varvec{u}_{h}$ and the DG spaces for ${p}_{h}$. The fully discretized linear momentum balance equation Eq. (7) can be defined using the following forms

$$\begin{aligned} {\mathcal {A}}_u\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \varvec{\psi }_{\varvec{u}} \right) = {\mathcal {L}}_u\left( \varvec{\psi }_{\varvec{u}} \right) , \quad \forall \varvec{\psi }_{\varvec{u}} \in {\mathcal {U}}_{h}^{\mathrm {CG}_{2}}\left( {\mathcal {T}}_{h}\right) , \end{aligned}$$

(39)

at each time step $t^n$, where

$$\begin{aligned} {\mathcal {A}}_u\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \varvec{\psi }_{\varvec{u}} \right) = \sum _{T \in {\mathcal {T}}_{h}} \int _{T} \varvec{\sigma }^{\prime }\left( \varvec{u}_{h}\right) : \nabla ^{s} \varvec{\psi }_{\varvec{u}} \, d V - \sum _{T \in {\mathcal {T}}_{h}} \int _{T} \alpha p_{h} {\mathbf {I}} : \nabla ^{s} \varvec{\psi }_{\varvec{u}} \, d V,\nonumber \\ \end{aligned}$$

(40)

and

$$\begin{aligned} {\mathcal {L}}_u\left( \varvec{\psi }_{\varvec{u}} \right) =\sum _{T \in {\mathcal {T}}_{h}} \int _{T} \varvec{f} \varvec{\psi }_{\varvec{u}} \, d V+\sum _{e \in {\mathcal {E}}_{h}^{N}} \int _{e} \varvec{t_D} \varvec{\psi }_{\varvec{u}} \, d S. \end{aligned}$$

(41)

Here, $\nabla ^{s}$ is a symmetric gradient operator. We then discretize Eq. (12) as

$$\begin{aligned} {\mathcal {A}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p \right) = {\mathcal {L}}_p\left( \psi _p \right) , \quad \forall \psi _p \in {\mathcal {P}}_{h}^{\mathrm {DG}_{1}}\left( {\mathcal {T}}_{h}\right) , \end{aligned}$$

(42)

for each time step $t^n$, where

$$\begin{aligned} {\mathcal {A}}_p\left( (\varvec{u}_{h}^{n}, p_{h}^{n}), \psi _p \right)= & {} \sum _{T \in {\mathcal {T}}_{h}} \int _{T} \frac{\alpha }{K} \mathrm {BDF}_1(\sigma _{v}) \psi _p \, d V\nonumber \\&+ \sum _{T \in {\mathcal {T}}_{h}} \int _{T} \left( \frac{1}{M}+\frac{\alpha ^{2}}{K}\right) \mathrm {BDF}_{1}\left( p_{h}^{n} \right) \psi _p \, d V \nonumber \\&+ \sum _{T \in {\mathcal {T}}_{h}} \int _{T} \varvec{\kappa } \nabla p_{h}^n \cdot \nabla \psi _{p} \, d V - \sum _{e \in {\mathcal {E}}_h^{I} \cup {\mathcal {E}}_{h}^{D}} \int _{e}\left\{ \varvec{\kappa }\nabla p_{h}^n\right\} _{\delta _{e}} \cdot \llbracket \psi _p \rrbracket \, d S \nonumber \\&- \sum _{e \in {\mathcal {E}}_h^{I} \cup {\mathcal {E}}_{h}^{D}} \int _{e}\left\{ \varvec{\kappa } \nabla \psi _{p}\right\} _{\delta _{e}} \cdot \llbracket p_h^n \rrbracket \, d S \nonumber \\&+ \sum _{e \in {\mathcal {E}}_h^{I} \cup {\mathcal {E}}_{h}^{D}} \int _{e} \frac{\beta }{h_{e}} {\varvec{\kappa }}_{{e}} \llbracket p_h^n \rrbracket \cdot \llbracket \psi _p \rrbracket \, d S, \end{aligned}$$

(43)

and

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_p\left( \psi _p \right)&:= \sum _{T \in {\mathcal {T}}_{h}} \int _{T} g \psi _{p} \, d V+ \sum _{e \in {\mathcal {E}}_{h}^{N}} \int _{e} q_{D} \psi _{p} \, d S \\&\qquad -\sum _{e \in {\mathcal {E}}_{h}^{D}} \int _{e} \varvec{\kappa } \nabla \psi _{p} \cdot p_{D} {\mathbf {n}} \, d S + \sum _{e \in {\mathcal {E}}_{h}^{D}} \int _{e} \frac{\beta }{h_{e}} {\varvec{\kappa }}_{{e}} \llbracket \psi _p \rrbracket \cdot p_D {\mathbf {n}} \, d S. \end{aligned} \end{aligned}$$

(44)

More details regarding block structure and solver algorithm could be found in Kadeethum et al. (2020, 2021).

Appendix B: Intrusive reduced order model by Galerkin projection

The starting point of the intrusive methodology based on a Galerkin projection is the end of Sect. 3.3. Let $\left\{ {\mathbf {w}}_{1}, \cdots , {\mathbf {w}}_{\mathrm {N}}\right\}$ denote the basis functions spanning ${\mathcal {U}}_{\mathrm {N}}$. Similarly, let $\left\{ q_{1}, \cdots , q_{\mathrm {N}}\right\}$ denote the basis functions spanning ${\mathcal {P}}_\mathrm {N}$.

The offline phase of the intrusive method does not require to carry out either the $L^2$ projection step as described in Sect. 3.4, or training the ANN as described in 3.5. Then, during the online phase (which effectively replaces the non-intrusive choice in Sect. 3.6), one carries out a time stepping as follows: given ${\widehat{\theta }}_k^u(t^0, \varvec{\mu })$, ${\widehat{\theta }}_k^p(t^0, \varvec{\mu })$, for every time step $t^n$ find reduced coefficients ${\widehat{\theta }}_k^u(t^n, \varvec{\mu })$, ${\widehat{\theta }}_k^p(t^n, \varvec{\mu })$ such that the reconstructed solutions

$$\begin{aligned} \begin{aligned} \widehat{\varvec{u}}_{h}^n := \widehat{\varvec{u}}_{h}\left( \cdot ; t^n, \varvec{\mu }\right) = \sum _{k=1}^{\mathrm {N}} {\widehat{\theta }}_k^u(t^n, \varvec{\mu }) {\mathbf {w}}_{k},\\ {\widehat{p}}_{h}^n := {\widehat{p}}_{h}\left( \cdot ; t^n, \varvec{\mu }\right) = \sum _{k=1}^{\mathrm {N}} {\widehat{\theta }}_k^p(t^n, \varvec{\mu }) q_{k}, \end{aligned} \end{aligned}$$

are solution to the following Galerkin method

$$\begin{aligned} {\left\{ \begin{array}{ll} {\mathcal {A}}_u\left( (\widehat{\varvec{u}}_{h}^{n}, {\widehat{p}}_{h}^{n}), {\mathbf {w}}_k \right) = {\mathcal {L}}_u\left( {\mathbf {w}}_k \right) , \quad \forall k = 1, \hdots , N,\\ {\mathcal {A}}_p\left( (\widehat{\varvec{u}}_{h}^{n}, {\widehat{p}}_{h}^{n}), q_k \right) = {\mathcal {L}}_p\left( q_k \right) , \quad \forall k = 1, \hdots , N. \end{array}\right. } \end{aligned}$$

More details on efficient assembly of such system by precomputation of time and parameter independent tensors can be found in reduced basis textbooks, see, e.g., Hesthaven et al. (2016).

Appendix C: Comparison with global basis for all unknowns

Coupled HM processes which we have considered throughout this manuscript, are problems characterized by two unknowns, namely bulk displacement ($\varvec{u}$) and fluid pressure (p). The standpoint we have taken in this manuscript is to compress the snapshots, compute reduced bases, and build neural networks for each variable separately (see Sect. 3). In the following, we refer to this method as the partitioned method. However, an alternative is to consider the two unknowns $(\varvec{u}, p)$ as a global (or monolithic) solution and derive global reduced bases from the compression of $\varvec{u}$ and p fields together. One advantage of doing that is that one only needs to train one neural network to approximate two variables. In the following, we will refer to this alternative method as the monolithic approach.

Fig. 20

Sensitivity analysis–errors of reconstruction solutions using 1000 testing $\varvec{\mu }$: a mean squared error (MSE) of displacement field ($\varvec{u}$) with $\mathrm {M} = 400$, $\mathrm {N} = 10$, $\mathrm {N}_{\mathrm{hl}} = 3$, and $\mathrm {N}_{\mathrm{nn}} = 7$, b mean squared error (MSE) of fluid pressure field (p) with $\mathrm {M} = 400$, $\mathrm {N} = 10$, $\mathrm {N}_{\mathrm{hl}} = 3$, and $\mathrm {N}_{\mathrm{nn}} = 7$. We fix $\mathrm {N_{int}} = 10$. The red squares represent outliers, and the box plot covers the interval from the 25th percentile to 75th percentile, highlighting the mean (50th percentile) with an orange line. We note that this figure presents the results of monolithic approach applied to Example 4. Please refer to Figure 19 for the results of partitioned method

We present the results of the monolithic method in Fig. 20. The problem setting including geometry, material properties, and boundary conditions are as in Example 4. Moreover, the ROM parameters are applied according to model 1: $\mathrm {M} = 400$, $\mathrm {N}_{\mathrm{int}} = 10$, $\mathrm {N} = 10$, $\mathrm {N}_{\mathrm{hl}} = 3$, and $\mathrm {N}_{\mathrm{nn}} = 7$, in Sect. 4.2.5. From these results, we clearly observe that the error of $\varvec{u}$ field in Fig. 20a is significantly higher than that of Fig. 19a. However, the results of p field between Figs. 20b and 19c are approximately similar.

This observation is as expected because the magnitude of $\varvec{u}$ and p are significantly different (six to seven orders of magnitude apart). The monolithic method compresses and builds a single set of reduced basis and neural networks for both variables. Due to the different orders of magnitude, the monolithic approach would only try to capture the p field, neglecting the $\varvec{u}$ field. Therefore, in practical applications, one would require considering the range of each state variable before training the monolithic method. To counteract this effect, in the loss function of the ANN, one could scale the contributions coming from the two-state variables trying to take into account the different orders of magnitude. However, such tuning may be cumbersome for problems with a large number of state variables. Instead, the partitioned approach results in improved performance for our HM problems and does not require any scaling.