Adaptive space-time model order reduction with dual-weighted residual (MORe DWR) error control for poroelasticity

In this work, the space-time MORe DWR (Model Order Reduction with Dual-Weighted Residual error estimates) framework is extended and further developed for single-phase flow problems in porous media. Specifically, our problem statement is the Biot system which consists of vector-valued displacements (geomechanics) coupled to a Darcy flow pressure equation. The MORe DWR method introduces a goal-oriented adaptive incremental proper orthogonal decomposition (POD) based-reduced-order model (ROM). The error in the reduced goal functional is estimated during the simulation, and the POD basis is enriched on-the-fly if the estimate exceeds a given threshold. This results in a reduction of the total number of full-order-model solves for the simulation of the porous medium, a robust estimation of the quantity of interest and well-suited reduced bases for the problem at hand. We apply a space-time Galerkin discretization with Taylor-Hood elements in space and a discontinuous Galerkin method with piecewise constant functions in time. The latter is well-known to be similar to the backward Euler scheme. We demonstrate the efficiency of our method on the well-known two-dimensional Mandel benchmark and a three-dimensional footing problem.


Introduction
Porous media problems have long-standing applications in subsurface modeling, groundwater flow, hydraulic fracturing, geothermal energy recovery, and nuclear waste storage.The resulting mathematical models yield coupled systems of partial differential equations (PDEs) from which one of the most well-known is the so-called Biot problem [18][19][20][21]71].Developments over the last two decades have led to more complicated models, due to additional physics that can be incorporated.Nonstationary, nonlinear coupled PDE systems are often obtained.They are also sometimes subject to inequality constraints as in multiphysics phase-field fracture in porous media [75,78].Despite advances in numerical solvers of iterative or multigrid type, high-performance paral-lel computing, and improved hardware performance, the computational cost for solving such multiphysics problems remains high.Specific examples in the field of poroelasticity are [5,9,10,22,42,44,53,54,76,77].
Model order reduction (MOR) and reduced-order modeling (ROM) techniques [14,15,30,37,47,56,57,65,68,80] provide one possible solution to significantly reduce the computational cost.Therein, the problem is split into two phases, an offline phase (solving a costly original high-fidelity model) to build a representative reduced basis and an online phase in which the reduced-order model is solved fast.This splitting comes with an additional cost, which pays off when the original model needs to be solved several (hundreds and more) times as for example in uncertainty quantification, inverse modeling, and optimal control.However, providing a robust and problem-specific reduced basis for nonlinear and general problems is a challenge.
Other model-order reduction approaches are based on solving the Galerkin problem with a separation of variables [2].For nonlinear problems, some robust methods have been proposed, but they require a highly intrusive and specific numerical framework [17], or less intrusive and more flexible approaches, but without a mathematical proof of convergence for nonlinear cases [16].Thus, there is a need for reduced-order numerical strategies which are flexible, general and for which the error accompanying the reduced solution space can be controlled for any quantity of interest.Additionally, a goal-oriented version of the proper generalized decomposition (PGD) has been applied to stationary problems in [43].In contrast to our work, an a priori computed dual solution was used to control the error of the primal reduced-order solution.
Thanks to our own recent advances in space-time modeling and goal-oriented spacetime error control for nonstationary, coupled problems [64], an incremental errorcontrolled proper orthogonal decomposition (POD) based ROM method, i.e.MORe DWR (Model Order Reduction with Dual-Weighted Residual error estimates), was suggested for the heat and elastodynamics equations in [32].Therein, satisfying results for adaptively refining the ROM basis according to some distributed-in-time goal functionals were obtained.
In the current work, the Biot system of poroelasticity, namely vector-valued displacements and a Darcy-type pressure equation, is considered and modeled in a space-time fashion.Specifically, the temporal discretization is based on discontinuous Galerkin (dG) finite elements, while the spatial discretization employs classical continuous Galerkin (cG) finite elements.This serves as a basis for space-time error-controlled adaptivity by applying the dual-weighted residual method (DWR) [8,12] to our ROM concepts.Importantly, the ROM error is computed with respect to a goal functional, which is often motivated from technical quantities of interest arising in physics, engineering and practical applications.The ROM updates are performed in an incremental fashion as recently suggested in [46] and do not require a full FOM (full-order model) run in advance.Thus, the ROM is updated on the fly using FOM solutions estimated only for one temporal element.Here, the singular value decomposition (SVD) for computing the leading components reduces to a truncated version and provides the snapshot matrix.An important point in coupled problems resorts to considering a monolithic version and computing one single SVD or, alternatively, splitting the SVD into its PDE components.In this work, we aim for the latter and split the displacements and pressure into separate SVDs in order to facilitate computing different POD sizes for the individual solution components.These develop-ments result into a final MORe DWR algorithm, which is tested for two-dimensional and three-dimensional configurations.First, the well-known Mandel benchmark is addressed.In a second numerical test, a poroelastic three-dimensional block is considered.Therein, on the one hand, typical cost factors such as speedup, FOM solves and ROM sizes are investigated.On the other hand, standard quantities of a posteriori error control such as effectivity indices and indicator indices are computed as well, demonstrating the performance of our overall proposed concepts.
The outline of this paper is as follows.In Sect. 2 the strong form problem statement is presented first.Then, a space-time formulation is derived.The dual-weighted residual method for goal-oriented error control is exposed, which requires a linear, backward-intime running, adjoint problem.Next, in Sect.

Problem formulation
In this work, we model a porous medium by coupling a vector-valued displacement equation (geomechanics) to Darcy flow in a poro-elastic medium.In the problem description, I := (0, T ) denotes the temporal domain and ⊂ R d a sufficiently smooth spatial domain with spatial dimension d ∈ {2, 3}.

Strong form of poroelasticity and function spaces
The governing equations for poroelasticity read [28,49] with the isotropic stress tensor and the (constrained specific) storage coefficient c ≥ c * > 0, which is assumed strictly positive in this paper, the Biot-Willis constant α ∈ [0, 1], the permeability tensor K , the fluid viscosity ν and the Lamé parameters λ, μ > 0. We notice that the constant c may depend on space, i.e., c(x), and is linked to the compressibility M > 0. When c tends to zero, numerical instabilities may arise (e.g., [58,60]) for which mixed methods or enriched Galerkin or discontinuous Galerkin finite elements in the numerical discretization should be employed [48,50,55,59,69].This coupled system of equations is also known as the Biot system [18][19][20][21]71].A rigorous mathematical analysis of this poroelasticity problem can be found in [67].
Assuming homogeneous Dirichlet boundary conditions for the displacement on D and inhomogeneous Neumann/traction boundary conditions the spatial function spaces are given by All boundary conditions, including pressure conditions as well, are carefully specified in Sect. 4 and the function spaces are redefined w.r.t. the modified boundary conditions, too.

Weak space-time form of the primal problem
We define the space-time function space X(I, V ( )) as with the dual space of V ( ) being denoted as V * ( ) = L(V ( ), R).Furthermore, we will use the notation In this notation, f • g represents the Euclidean inner product if f and g are scalar-or vector-valued and represents the Frobenius inner product if f and g are matrices.
We can now derive the space-time variational formulation for the (primal) poroelasticity problem.By integration by parts, the variational formulation reads: Here, n denotes the normal vector on the Neumann boundary N .
For the time discretization, let Then, in an intermediate step similar to spatial discontinuous Galerkin (dG) method, we can define a broken function space that allows for discontinuities at the temporal grid points: Due to these discontinuities, the limits of a function f at time t m from above and from below are and the jump of the function value of f at time t m is Accounting for the discontinuities in time, we thus need to solve the problem: The bilinear form and right-hand side read and where Remark 2.1 This space-time formulation can also be easily extended to a multirate-intime setting with different timestep sizes for displacement and pressure [63]; a rigorous analysis starting directly with the backward Euler time discretization and mixed spaces for flow and conformal Galerkin for geomechanics is provided in [3].
For the projection of the problem onto the POD vectors, the block structure of the linear system is exploited to build separate bases for displacement and pressure.We can decompose the bilinear form into the four following blocks

dG(0), i.e., backward Euler, time discretization of the primal problem
Using the space-time formulation of the poroelasticity problem (2), a dG(0) time discretization is derived by using discontinuous, piecewise-constant finite elements in time; see e.g., [79].This results into a backward Euler scheme.Then, on each temporal element I m = (t m−1 , t m ), we have temporally constant functions (u m , p m ) ∈ V ( ) such that Using piecewise-constant functions f in time, we can insert the relations and arrive at the time-stepping scheme: Here, k := t m − t m−1 denotes a (constant) timestep size.
Remark 2.2 Instead of starting from the space-time formulation, we can also directly derive the backward Euler time discretization from the strong formulation of poroelasticity (1).The main difference is that now we do not scale the mechanics equation by the timestep size k.This is motivated by the displacement equation being quasi-static.Hence, in traditional time discretizations the time integral is omitted for this equation.For a consistent mathematical description, the previous time-stepping scheme fits better in the space-time setting, but in the actual implementation we can neglect the time step size k in the mechanics equation.We then have: We will use this formulation for our numerical tests since it is equivalent to the dG(0) time discretization.

Weak space-time form of the adjoint problem
The MORe DWR method [32] measures the adjoint sensitivity of the primal reducedorder solution with respect to some quantities of interest.Let a goal functional J : X(T k , V ( )) → R of the form be given, which represents some physical quantity of interest (QoI).The adjoint equations for porous media are derived, following [79].Since the problem is linear, we just need to switch trial and test functions in the bilinear form and replace the primal with the adjoint solution.The forcing function on the right-hand side is being substituted by the goal functional J , which is justified by the Lagrangian formalism, c.f. Sect.3. The space-time formulation of the adjoint problem thus reads: More concretely, the bilinear form reads where The time derivatives are moved from the test function to the adjoint solution by integration by parts in time e.g., on the time-continuous level we use integration by parts for The adjoint problem runs backward in time and starts with a final time condition, which depends on the quantity of interest.In particular, we have i.e., Z(T ) = 0 for J 2 = 0.The integration by parts can be applied to each temporal element, which leads exemplarily to Overall, we can rewrite the left-hand side of the adjoint problem as

dG(0), i.e., backward Euler, time discretization of the adjoint problem
A dG(0) time discretization is again used for the adjoint problem.Then on each temporal element Using piecewise-constant functions f in time, the relations hold and give the adjoint time-stepping scheme:

MORe DWR: model order reduction with dual-weighted residual error estimates
The MORe DWR approach aims at evaluating the quantities of interest using a reducedorder model, while guaranteeing the error due to this approximation.Thus, the difference between the reduced-order-model (ROM) solution controlled by employing a dual-weighted residual method [8,11,12].To this end, we obtain an optimization problem in which the model error between FOM and ROM, both measured in terms of some goal functional J (•), shall be minimized: subject to the constraint that the variational formulation of the poroelasticity problem ( 2) is satisfied by ).For more information on space-time error control, we refer the reader to [64,66,70,79].We focus in this work on the enrichment of the reduced basis depending on the temporal evolution of the quantities of interest.The reduced basis is refined in a goal-oriented way to accurately and efficiently compute the solution over the whole temporal domain.In principle coarsening would also be possible, but is not the objective in this work.For coarsening, we would need to follow the work of Meyer and Matthies [52].

Space-time dual-weighted residual method
The two Lagrange functionals for the constrained optimization problem (5) are defined as with ∈ {FOM, ROM}.The stationary points (U FOM , Z FOM ) and (U ROM , Z ROM ) of the Lagrange functionals L FOM and L ROM need to satisfy the Karush-Kuhn-Tucker first-order optimality conditions.

Primal problem
Firstly, the stationary points are solution to the following primal problems The corresponding solutions denoted U FOM and U ROM are the primal solutions.We observe that the primal solution can be obtained by solving the original problem, i.e., the quasi-static poroelasticity problem, forward in time as presented in Sect.2.3.

Adjoint problem
Secondly, the stationary points must also satisfy the following adjoint (or dual) problems which give the adjoint solutions Z FOM and Z ROM .Note that here we already use that the poroelasticity equations are linear, which ensures that Hence, the adjoint solution is obtained by solving which is backward in time as discussed in Sect.2.5.
Remark 3.1 For linear goal functionals, the right-hand side of the adjoint problem (6) reduces to and the adjoint problem does not depend on the primal solution anymore.

Error estimator
As a first result, we have Proposition 3.1 The computable ROM error estimator is given by Proof The proof follows from Theorem 4.2 from [32].It holds The previous result is assembled on each temporal element I m separately to localize the error in time.More concretely, for poroelasticity the temporally-localized error estimator reads: Proposition 3.2 Given the primal problem (2), the goal functional (3), and the corresponding adjoint problem (6) with the adjoint bilinear form (4), and η from Proposition 3.1, the temporally-localized error estimator is given by for m = 1, . . ., M.Moreover, from (7), we have The superscripts indicating reduced-order-model solutions for both the primal and dual solutions are omitted for a clearer presentation.The (reciprocal) effectivity index [6], which is the ratio between the estimated and the true errors, i.e., is used to measure the quality of our error estimator.We desire I eff ≈ 1 since then the error estimator can reliably predict the reduced-order-modeling error.We refer to [31] for two-sided proofs of discretization errors measured in goal functionals using saturation assumptions.We also observe that the effectivity index I eff is close to one in the numerical tests in Sect. 4. Finally, the quality of the adaptive enrichemnt is measured in terms of the indicator index [62]:

Error estimator-based ROM updates
In this section, we give a compact summary of our novel approach of a goal-oriented incremental reduced-order model from [32] and discuss its extension to porous media settings.
In the MORe DWR method, we marry a reduced-order model with a dual-weightedresidual-based error estimator and an incremental version of the POD algorithm.The goal is then to use the error estimator to identify when the solution behavior is not captured accurately by the reduced basis, such that it is incrementally enriched on-the-fly with new full-order-model snapshots.The approach is particularly cheap because the FOM snapshots are computed only for one temporal element.
In more detail, we apply our findings on error control of Sect.3.1 to a backward Euler reduced-order model of poroelasticity.Further, an incremental basis generation is mandatory for the method to reduce computational operations and thus to be fast, which we realize by means of an incremental SVD.The incremental SVD is presented in Sect.3.3.1.In this context, we also introduce the incremental POD as a trimmed version of the incremental SVD.Subsequently, the overall MORe DWR framework is depicted in Sect.3.4, where all the ingredients are assembled and the final algorithm is presented.
In summary, our novel approach avoids a computationally heavy offline phase and directly solves the reduced model.Our approach can be thought of as a replacement for traditional full-order-model simulations, since we highly reduce the computational cost, while controlling the error between the full-order (FOM) and the reduced-order model (ROM).

Incremental proper orthogonal decomposition
The proper orthogonal decomposition (POD) [7,15,25,27,30,39,45,47,61,68,80] can be employed to project the differential equation onto a lower-dimensional solution manifold.In traditional applications of POD, one first collects full-order solution snapshots and finds the POD basis by computing the singular value decomposition Here, n is the number of spatial degrees of freedom, M + 1 is the number of FOM snapshots, and d ≤ min(n, M + 1).The SVD is then being truncated to the size N d using a retained energy criterion, cf.[36,38,47] e.g., For coupled problems with several PDEs, like poroelasticity, we could again perform the SVD of the snapshot matrix Y consisting of the entire solution vectors but we choose to create separate SVDs for displacement u(t i ) and pressure p(t i ) such that POD bases with various sizes can be used for the different solution components.Instead of building repeatedly a reduced basis from scratch, we suggest updating an already existing truncated SVD (tSVD) or solely its left-singular (POD) vectors according to modifications of the snapshot matrix without recomputing the whole tSVD or requiring access to the snapshot matrix [32,46].The POD becomes incremental by appending additional snapshots to the initial snapshot matrix.This difference between the classical SVD and incremental SVD (iSVD) is illustrated in Fig. 1.In this context, we rely on the general approach of an additive rank-b modification of the SVD, mainly developed by [23,24] and applied to the model-order reduction of fluid flows in [46].
Let Y ∈ R n× m be a given snapshot matrix that includes m > 0 snapshots.Usually, m is equal or connected to the number of already computed time steps.Further, we have the rank-N tSVD T of the matrix Y .Additionally, let b ∈ N newly computed snapshots {U 1 , . . ., U b } be stored in the bunch matrix We aim to compute the tSVD that is updated by the information contained in the bunch matrix B according to

The MORe DWR algorithm
The overall procedure is illustrated in Fig. 2.

Fig. 2 Schematic representation of MORe DWR algorithm
The aim is to solve the reduced-order model and adaptively enrich the reduced basis by means of the iPOD with full-order solutions until a given estimated error tolerance is reached for the chosen goal functional.The approach is designed to work without any prior knowledge or exploration of the solution manifold, while also attempting to minimize the number of full-order solves.The error associated with the reduced solution is estimated on each temporal element.The new full-order solution used for the basis enrichment is computed on the temporal element with the largest error.So, two full-order solves are conducted for each enrichment.We observe that this procedure is perfectly compatible with the adaptive basis selection based on DWR estimates presented by Meyer and Matthies in [52] to reduce the dimension of the reduced space.Thus, if incorporated it would be possible to either enrich or downsize the reduced basis according to the problem statement.The resulting approach is detailed in Algorithm 2. The reduced primal and dual solutions are used to provide a fast error estimation.Thus, the incremental ROM has both a primal and a dual reduced basis as input in the form of the reduced basis matrix.

Initialization
If no prior information is available, the bases are initialized by computing the primal and dual solution snapshots on the first temporal element.If prior reduced bases are available e.g., from a previous simulation on a different parameter configuration or some (very cheap) coarse-grid solution in the context of a multigrid idea, prior bases are used as an initial guess and are then altered by the goal-oriented adaptation.Thus, the MORe DWR approach is totally appropriate for reduced-order modeling of parameterized problems.The method can be an efficient substitute for the full-order model, as it builds a reduced basis tailored to a quantity of interest with a minimum of FOM solves.
Nevertheless, in the scenario without prior knowledge, we observed an issue when exclusively the first dual FOM solution (t = T ) is utilized for the basis generation.We identified a significant discrepancy between the error estimates and the real errors, which is demonstrated and further described in Fig. 8.To mitigate this issue, further assessments utilizing the final dual snapshots for the basis enrichment at the first temporal elements (t = 0) were conducted, which consistently improved the accuracy of the error estimates.
Nonetheless, a practical constraint is posed, as these snapshots are not readily accessible without completing the whole dual FOM solution process.Thus, a computationally efficient surrogate to bypass the need for the full-order solution process is required.For this, the dual FOM is substituted by the dual ROM, which is already computed during the enrichment process and only the last dual solution steps are conducted by the dual FOM.Given that the inconsistency in the estimate predominantly emerges during the initial enrichment iterations, we have optimized the process by confining the additional dual ROM enrichment solely to these early iterations for efficiency.

Adaptive algorithm
At each MORe DWR iteration, the primal and dual reduced bases are enriched for the temporal element m max corresponding to the highest relative error comparing the local relative error associated with any temporal element m.Because the full-order solutions are not available, the error is normalized according to The adaptive MORe DWR algorithm is stopped according to the global relative error estimator with TOL rel a user-chosen threshold.This ensures that the relative error over the entire temporal domain is below a given tolerance, whereas in [32] we ensured that the local relative error is below a given tolerance.However, enforcing the accuracy locally is a more greedy approach, which leads to more full-order-model solves.

Remark 3.2
The additional enrichment of the dual space to provide an accurate estimation of the error is omitted in Algorithm 2 and in Fig. 2 for the sake of clarity.This dual enrichment takes place after line 12 in Algorithm 2 and after enriching the primal and dual reduced bases in Fig. 2. The supplementary enrichment is based on the FOM solution of the dual problem for a few time steps, corresponding to the few last time steps of the dual problem i.e., the few first time steps of the primal problem.As the dual problem is backward in time, it starts from the ROM solution at t = t l with l ∈ N and the solution is computed for t l−1 , t l−2 , …, t 0 .The dual reduced basis is then enriched with these additional snapshots.This procedure is applied only in the few first MORe DWR iterations prior to the update of the reduced system components and the error estimator.

Numerical tests
The MORe DWR framework is numerically substantiated on two different problem configurations.The Mandel problem [1,26,29,51,73] (see also [33,35,40,50,78] and more recently also for nonlinear poroelasticity [73]) is first considered as a benchmark for two-dimensional poroelasticity and a linear goal functional.In Mandel's problem one Compute error estimate: η rel 7: if η rel > tol then 8: Find temporal element with maximum error: Solve primal full-order system on I m max for Update primal reduced basis: Solve dual full-order system on Update dual reduced basis: 13: Update reduced system components and error estimator observes the so-called Mandel-Cryer effect [29,72] of a non-monotonic pressure evolution: first increasing pressure, followed by decreasing pressure.The second numerical test is a three-dimensional footing problem inspired by [34].Our computations have been performed on an AMD EPYC 7H12 64-Core Processor.The FEM codes have been written in FEniCS [4] and the reduced-order modeling has been performed with NumPy [41] and SciPy [74].

Mandel problem in 2D
Let := (0 m, 100 m) × (0 m, 20 m) and I := (0 s, 5 000 000 s) ≈ (0 ,58 d) with boundaries as shown in Fig. 3.The initial and boundary conditions are given by   This results in the spatial function spaces (then modified in Section 2.1) being given by V u := {u ∈ H 1 ( ) 2 | u y = 0 on bottom , u x = 0 on left }, The parameters for Mandel's problem are summarized in Table 1.In space, we use Taylor-Hood elements, i.e., quadratic finite elements for the displacements u and linear finite elements for the pressure p.For this numerical test, a fixed spatial mesh with 80 spatial cells in the x-direction and 16 spatial cells in the y-direction is designed, which leads to an isotropic mesh with 10, 626 DoFs for displacement and 1, 377 DoFs for pressure.For the temporal discretization, the end time is T = 5 000 000 s ≈ 58 days with 5, 000 temporal elements i.e., the chosen time step size k = 1000 s.For the quantity of interest, we choose the time-integrated pressure acting on the bottom boundary i.e., J (U ) := Further, for the MORe DWR enrichment, we choose 1 − 10 −7 and 1 − 10 −11 for the primal displacement and pressure energy tolerances (11).Accordingly, the dual energy tolerances are both set to 1 − 10 −9 .As argued in Sect.3.4, it is crucial to have a dual space large enough for estimating accurately the error.Thus, the reduced dual basis is enriched during the first MORe DWR loops from the snapshots of the full-order dual solution on the first time steps [t 0 , t −1 ] i.e., for the last time points of the dual problem which is  First, the adaptive enrichments of the method are shown in Fig. 4. The true relative error between the reference full-order solution and the reduced-order solution is given by Here, we observe a clear decrease with the iterations of the MORe DWR algorithm; see Fig. 4a.In contrast to the original MORe DWR paper [32], the relative error is now much closer to or even slightly exceeds the error tolerance.This can be explained by the fact that we do not enforce the error tolerance on each temporal element, but only on the whole integrated time domain.The estimated relative error is very close to the true error after iteration 2. Thereafter, only slight deviations can be noticed locally.On the contrary, a large discrepancy can be observed between the true and estimated errors for the first iterations, where the estimate yields significant underestimations due to the inaccuracy of the dual reduced space.The target tolerance of 0.1% is reached after 36 iterations.For larger error tolerances, the algorithm would terminate sooner e.g., for 1% relative error, we need 29 iterations, for 5% error we need 22 iterations and for 20% error only 13 iterations are required (see Table 2).The enrichment of the POD basis is illustrated in Fig. 4b.The primal reduced spaces increased progressively and regularly with the iterations.We observe that the primal pressure solution requires significantly more modes than the primal displacement.This implies that the choice of separating the bases for displacement and pressure is beneficial.Interestingly, the dual displacement solution (37 modes) requires more POD modes than the primal displacement solution (5 modes), Fig. 5 Goal functionals for Mandel's problem and different relative error tolerances whereas the dual pressure solution (26 modes) requires less POD modes than the primal pressure solution (33 modes).We also see that the sizes of the dual bases increase faster than in the primal case during the first 5 iterations as we enforce an additional dual basis enrichment.But, its slope decreases after these initial enrichments.
The time trajectories of the goal functional from the full-order solution U h and the reduced-order solution U N are compared for each temporal element in Fig. 5.For these comparisons, we choose relative error tolerances between 0.1% and 20%.Both quantity of interest trajectories are almost indistinguishable for tolerances between 0.1% and 2%.Some differences appear between the full-order and the reduced-order results for larger relative error tolerances.They might be seen as a too large error.However, for certain real-world applications an error of e.g., 10% might still be acceptable for an appealing computational cost.We can see that choosing the relative error tolerance, the MORe DWR approach allows to control the quality of the approximation of the true quantity of interest, even within a reduced-order solution framework.
Table 2 gives an overview of the results obtained for the Mandel problem using the MORe DWR algorithm with various error tolerances TOL rel between 0.1% and 20%.For each tolerance value, the relative true error e rel is given as well as the computational speedup, the total number of FOM solves, the POD basis sizes for the primal displacement/pressure and the dual displacement/pressure, the effectivity index I eff from (9) and the indicator index I ind defined in (10).Here, the number of FOM solves sums up all primal and dual solves and the basis sizes are shown in the pattern: primal displacement / primal pressure + dual displacement / dual pressure.The relative error as well as the speedup increase with a rise in tolerance.The relative error is rather close to the error tolerance, but sometimes slightly exceeds it.Indeed, here error estimates are employed, not error bounds.This can be explained by the fact that we do not enforce the error tolerance on each temporal element, but only on the whole integrated time domain.In contrast, enforcing the error tolerance on each temporal element leads to much lower error than the error tolerance [32].The speedup is explained by the decreasing amount of FOM solves and smaller POD bases for both the primal and dual problems, as well as displacement and pressure w.r.t. the given tolerance.The effectivity index and the indicator index are close to 1 as expected for a linear PDE and goal functional.This validates the usage of dual-weighted residual error estimates to control the reduced-order-modeling error because the estimated error accurately approximates the true error.
The initial and boundary conditions are given by The material parameters as listed in Table 1 are used again.As before, Taylor-Hood elements, namely quadratic finite elements in space for the displacement u and linear finite elements for the pressure p, are employed.The spatial mesh is fixed, it comprises 16 spatial cells in each direction i.e., an isotropic mesh with 107, 811 DoFs for displacement and 4, 913 DoFs for pressure.The temporal domain from 0 to the end time T = 5 000 000 s ≈ 58 d is discretized with 5, 000 temporal elements, which leads to a time step size k = 1000 s.For the goal functional, we choose the time-integrated pressure acting at the compression boundary i.e., J (U ) := I compression p dx dt.
A linear system of equations with 112, 724 unknowns is obtained, which has to be solved for each temporal element.The previous 2D Mandel problem was solved with a direct solver.Due to the much larger size of the equation system, we now have to resort to an iterative solution scheme for maintaining a low computational cost.The equation system is not symmetric and the memory for the calculations is not limited; thus we choose the generalized minimal residual method (GMRES).Specifically, SciPy [74] implementation of GMRES is utilized with a tolerance for convergence of 5 • 10 −8 .For the preconditioning of the system, we have tested a diagonal Jacobi preconditioner as well as the smoothed aggregation algebraic multigrid (SA-AMG) method utilizing the implementation provided by [13].Both approaches have been benchmarked on the footing problem for the first 500 temporal elements.The SA-AMG ansatz needs a mean iteration count of 6.33 for convergence per temporal element solve, while the Jacobi preconditioner requires 283.81 iterations to reach the same accuracy.Although the number of iterations is significantly lower for the SA-AMG method, the mean wall time per temporal element solve is 93.37 s for the SA-AMG approach, but only 3.89 s for the Jacobi preconditioner.The Jacobi preconditioner is significantly faster than the SA-AMG method in our case, most likely due to easier parallelization capabilities on multicore machines.Thus, the Jacobi preconditioner is employed for this application case.
Similar to Mandel's problem, we choose 1 − 10 −7 and 1 − 10 −11 for the primal displacement and pressure energy tolerance (11).Again, the dual energy tolerances are both set to 1 − 10 −9 .As argued in Sect.3.4, the dual spaces are enriched with additional dual full-order solutions.In this case, we perform this additional enrichment for the first 8 MORe DWR iterations and include the dual solutions for the last five time steps of the dual problem i.e., the time steps [t 0 , t 4 ] of the primal problem discretization.
The adaptive enrichment of the method is illustrated in Fig. 7.The relative error between the full-order and reduced-order solutions, as well as the error estimator is shown in Fig. 7a.We observe a similar behavior for the first iterations as in Mandel's problem: a discrepancy between the estimate and the error, with a severe underestimation of the error.Thereafter, both quantities align for the rest of the iterations.In addition, we see a sharp decline after 18 iterations that is followed by a non-monotonic behavior, highlighting the greedy enrichment process of the MORe DWR method.
The evolution of the size of the primal/dual POD bases for displacement/pressure is shown in Fig. 7b.We see again a steep increase in the dual POD sizes in the first iterations Next, in Fig. 8 we compare the latter results with the MORe DWR method neglecting the additional dual reduced space enrichment proposed in Sect.3.4.A minimum number of 20 iterations is enforced to ensure the consistency of the results because it is expected from the previous study that the error is severely underestimated during the first iterations without extra enrichment (Fig. 7a). Figure 8a compares the relative error between the full-order and reduced-order solutions with its estimate over the course of MORe DWR iterations.In contrast to the results in Fig. 7a, we observe that the first phase of discrepancy between the estimate and the error, with the typical severe underestimation of the error, lasts much longer with more variation.The true error and its estimate align after 35 iterations.Thereafter, the error estimate is as accurate as when using the extra dual enrichment.Thus, the extra enrichment of the dual basis appears necessary for an accurate error estimate from the early stage of the algorithm.The ROM solution that complies to a tolerance of 0.5% is obtained with a total of 104 FOM evaluations (52 primal plus 52 dual evaluations) without the dual enrichment, whereas only 80 FOM evaluations (20 primal plus 20 dual plus 40 extra dual) were required with the enriched dual basis.Thus, although we added dual FOM evaluations in the first iterations of the extra-enriched approach, the trade-off seems beneficial, as this approach reduces the total number of iterations and thus the computational costs.In addition, the development of the POD bases in case of no additional dual enrichment is shown in Fig. 8b.In contrast to Fig. 7b, except for the primal displacement, all the POD bases start growing at nearly the same rate.The dual pressure basis size stagnates after 42 iterations, and the primal displacement has, after 7 iterations, a constant size of 4. Comparing the POD basis sizes for both approaches for a relative tolerance of 0.5%, the final dual POD basis sizes are nearly the same.In contrast, in the case of the extra enrichment, the primal pressure POD size of 19 is considerably smaller than the 47 modes needed when neglecting this enrichment.
The time trajectory of the goal functional from the reduced-order solution is compared with the trajectory given by the full-order solution for error tolerances between 0.1% and 20% in Fig. 9.The trajectories of quantities of interest are nearly identical for a tolerance of 0.1%.Differences emerge between the full-order and reduced-order goal functionals for larger tolerances.In essence, the reduced-order solution reaches the accuracy chosen by the user.
The time trajectories of the goal functional from the reduced-order solution is compared with the trajectory given by the full-order solution for error tolerances between 0.1% and 20% in Fig. 9.The trajectories of both quantities of interest are nearly identical for a tolerance of 0.1%.Differences emerge between the full-order and reduced-order goal functionals for larger tolerances.In essence, the reduced-order solution reaches the accuracy chosen by the user.
In Table 3, we give an overview of the simulation results of the three-dimensional footing problem for different error tolerances between 0.1% and 20%.Tolerances are chosen based on their potential applicability in real-world scenarios.Lower tolerances could necessitate too many POD modes for the reduced-order model, rendering it less efficient.The tolerances are met besides in the case of 0.5%, where the error is slightly underestimated.We observe that with a rise in the tolerance, the relative error as well as the speedup increase.We note that the algorithm has the same behavior for the tolerances 1%, 2% and 5%.Indeed, these tolerances are reached within the same number of MORe DWR iterations, yielding all other quantities to be the same.This behavior can be explained by the sharp decline in the error after 19 MORe DWR iterations, c.f. Figure 7a.The evaluation of the dual-weighted residual error estimates by means of the effectivity and indicator indices shows that the estimate integrated over the temporal domain is almost perfect, as the effectivity index is near 1.0 in all cases.However, the indicator indices vary from 1.00 up to 3.44 which translates to inaccurate temporal error localizations in the case of high indicator indices.

Conclusion
We developed an adaptive incremental reduced-order (MORe DWR) model for singlephase flow problems in porous media, namely the Biot system.To this end, we further developed and extended our prior work [32].Therein, the following ingredients are combined.First, a space-time formulation and Galerkin space-time finite element discretization for the coupled problem.Second, a ROM method based on an incremental POD.Third, a goal-oriented a posteriori error estimator using the dual-weighted residual method for adaptively refining the ROM bases.The efficiency of the methodology has been demonstrated, by using the incremental ROM for the two-dimensional Mandel benchmark and a three-dimensional footing problem.The algorithm is stopped according to the total relative error considering the whole temporal domain.The speed-ups offered by this approach compared with the reference full-order approach vary from 8.5x in two dimensions for the smaller relative tolerances, to 26.2x for the three-dimensional case with large relative tolerances.The accuracy of the error estimator has been improved by additional enrichment of the dual reduced bases.The true error and the error estimator almost coincide with effectivity indices close to one, which makes the error estimates reliable.It is important to note that the reduced-order model starts with no other prior information than a single primal snapshot and a dual FOM snapshot.Hence, the full-order model needs to be queried in the online stage to gather snapshots for basis enrichment.Consequently, our approach could be further sped up by efficient full-order model solvers.Therefore, for the three-dimensional test, we used GMRES as our iterative solver and compared a simple diagonal Jacobi preconditioner and an algebraic multigrid preconditioner.An aspect for future work is the combination of temporal adaptivity and reduced basis adaptivity.The MORe DWR algorithm already requires the solution of the primal and dual problems.This additional computational effort can be reused to refine the temporal mesh adaptively using the dual-weighted residual method and could reduce the number of time steps, thus further speeding up the computations.

Remark 2 . 3
To be consistent with the backward Euler time discretization in Remark 2.2, we again get rid of the timestep size k in displacement related terms and obtain

Fig. 4
Fig. 4 Adaptivity in the MORe DWR method of Mandel's problem

Fig. 7 Fig. 8
Fig. 7 Adaptivity in the MORe DWR method of the 3D footing problem

Fig. 9
Fig.9 Goal functionals for the 3D footing problem, depending on the tolerance in the goal functional 3, the MORe DWR (Model Order Reduction with Dual-Weighted Residual error estimates) concepts are explained.Goal-oriented error control for refining the ROM basis is also discussed.In Sect. 4 two numerical examples are investigated.The first test is the two-dimensional Mandel benchmark.The second test is a three-dimensional footing problem.

Table 1
Parameters in Mandel's problem

Table 2
Performance of MORe DWR method for the Mandel problem depending on the tolerance in the goal functional formulated backward in time.Here, it is enriched for the first five MORe DWR iterations based on the snapshots corresponding with [t 0 , t 4 ] in the primal temporal discretization.

Table 3
Performance of MORe DWR method for the 3D footing problem, depending on the tolerance in the goal functional