Space–time adaptive hierarchical model reduction for parabolic equations

Background Surrogate solutions and surrogate models for complex problems in many fields of science and engineering represent an important recent research line towards the construction of the best trade-off between modeling reliability and computational efficiency. Among surrogate models, hierarchical model (HiMod) reduction provides an effective approach for phenomena characterized by a dominant direction in their dynamics. HiMod approach obtains 1D models naturally enhanced by the inclusion of the effect of the transverse dynamics. Methods HiMod reduction couples a finite element approximation along the mainstream with a locally tunable modal representation of the transverse dynamics. In particular, we focus on the pointwise HiMod reduction strategy, where the modal tuning is performed on each finite element node. We formalize the pointwise HiMod approach in an unsteady setting, by resorting to a model discontinuous in time, continuous and hierarchically reduced in space (c[M(\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf{M}$$\end{document}M)G(s)]-dG(q) approximation). The selection of the modal distribution and of the space–time discretization is automatically performed via an adaptive procedure based on an a posteriori analysis of the global error. The final outcome of this procedure is a table, named HiMod lookup diagram, that sets the time partition and, for each time interval, the corresponding 1D finite element mesh together with the associated modal distribution. Results The results of the numerical verification confirm the robustness of the proposed adaptive procedure in terms of accuracy, sensitivity with respect to the goal quantity and the boundary conditions, and the computational saving. Finally, the validation results in the groundwater experimental setting are promising. Conclusion The extension of the HiMod reduction to an unsteady framework represents a crucial step with a view to practical engineering applications. Moreover, the results of the validation phase confirm that HiMod approximation is a viable approach.

models. Surrogate solutions are generally formalized with a reduction of the size of the finite dimensional solution, as in the reduced basis approach [2], or in the proper orthogonal decomposition (POD) [3] and proper generalized decomposition (PGD) methods [4,5].
Surrogate models directly replace the reference model via a simplified formulation as with a geometric multiscale modeling [6,7] or with compressed sensing [8]. This is usually accomplished by taking advantage of specific features of the problem at hand, such as a prevalent direction in the involved dynamics rather than in the geometry of the computational domain. This is exactly the criterion exploited to settle the hierarchical model (HiMod) reduction proposed in [9,10]. The HiMod technique derives enriched 1D surrogate models to describe phenomena characterized by a leading dynamics albeit in the presence of locally significant transverse features. In particular, the description properties of purely 1D models are enhanced by keeping track of the transverse dynamics in the reduced model. This is achieved by enriching a finite element discretization of the mainstream with a modal representation of the secondary dynamics. This strategy leads to a 1D finite element model with ad-hoc coefficients that implicitly include the generally non-constant description of the transverse dynamics. The possibility of locally tuning the modal expansion to match spatial heterogeneities represents one of the main strengths of the HiMod approach [11].
In this paper, we focus on the pointwise HiMod reduction strategy proposed in [12], where the modal tuning is performed on the finite element nodes. For this reason, the pointwise approach turns out to be the most flexible one among the available HiMod procedures [13], being suited to model both localized and widespread dynamics. In particular, with a view to practical applications, we extend the pointwise HiMod formulation to an unsteady setting by resorting to a discretization discontinuous in time. We generalize the cG(s)-dG(q) formulation in [14][15][16] to the HiMod setting, by defining a reduced solution that we denote by c[M(M)G(s)]-dG(q) approximation. We replace the full model with a solution that is continuous in space and discontinuous in time. It is obtained via a Galerkin spatial approximation that combines finite elements of degree s with the modal expansion identified by the index M, and via discontinuous piecewise polynomials of degree q in time.
The selection of the modal distribution as well as of the space-time discretization represents a crucial step of the HiMod reduction. For this reason, we introduce a preprocessing phase to automatically identify the HiMod solution, for fixed values of s and q. The final outcome of this phase is a table that identifies the time partition and then, for each time interval, selects the corresponding 1D finite element mesh together with the associated modal distribution. We call this table HiMod lookup diagram. To this purpose, we resort to an adaptive procedure based on an a posteriori analysis of the global (modeling plus space-time discretization) error. We rely upon a goal-oriented setting [17][18][19], so that the prediction of the c[M(M)G(s)]-dG(q) model is driven by a physical quantity of interest.
The estimator for the global error consists of a modeling and of a discretization contribution, which are preserved distinct [11,[20][21][22]. This represents a crucial property with a view to a global adaptation algorithm. In particular, the modeling estimator is a generalization of the goal-oriented hierarchical a posteriori error estimator derived in [11] The full setting We introduce the reference parabolic model we aim at reducing via an adaptive spacetime model reduction procedure. A standard notation is adopted for the Sobolev spaces associated with the spatial independent variable only, as well as for the space of the functions bounded almost everywhere [28]. Concerning a space-time dependence, we introduce the spaces L 2 (0, T ; W ) = v : (0, T ) → W : T ] → W continuous : ∀t ∈ [0, T ], �v(t)� W < +∞ , where W denotes a generic Hilbert space, with � · � W the associated norm [29].

The problem
We select as model to be reduced the unsteady problem where ⊂ R d (d = 2, 3) is the computational domain, Ŵ D and Ŵ N constitute a measurable non-overlapping partition of ∂� such that ∂� = Ŵ D ∪ Ŵ N and • ŴD ∩ • ŴN = ∅, I = (0, T ) is the time window of interest, and L is a generic second-order elliptic operator with diffusive contribution given by −∇ · (D∇u) so that D∇u · n ≡ ∂ ν u is the conormal derivative of u, n being the unit outward normal vector to ∂�. Concerning the data, we choose the source f ∈ L 2 (0, T ; L 2 (�)), the diffusivity tensor D = [d ij ] ∈ [L ∞ (�)] d×d such that the uniform ellipticity condition holds, the initial datum u 0 ∈ L 2 (�), and the Neumann datum g ∈ L 2 (0, T ; L 2 (Ŵ N )). In the next section, further requirements are added on the computational domain as well as on the boundary conditions in view of the HiMod procedure.
We consider the weak formulation associated with (1), given by: find u ∈ V = L 2 (0, T ;  with u(x, 0) = u 0 (x), and where a(·, ·) : H 1 Ŵ D (�) × H 1 Ŵ D (�) → R is the bilinear form associated with operator L, here assumed continuous and coercive. Problem (2) represents the full problem, with u the full solution.

The computational domain
Problems suited to a HiMod reduction are defined on domains characterized by a prevalent dimension and the leading dynamics is aligned with such a dimension.
Thus, we assume to coincide with the d-dimensional fiber bundle � = x∈� 1D {x} × γ x , where 1D is the supporting 1D fiber described by the independent variable x and aligned with the dominant dynamics, while γ x ⊂ R d−1 denotes the transverse fiber that is, in general, a function of x and parallel to the transverse dynamics. For the sake of simplicity, we assume 1D ≡]x 0 , x 1 [ to be rectilinear and we refer to [30] for the more general case of a curved supporting fiber. We partition the boundary ∂� of into three disjoint sets, Ŵ 0 = {x 0 } × γ x 0 , Ŵ 1 = {x 1 } × γ x 1 and Ŵ * = x∈� 1D ∂γ x , such that ∂� = Ŵ 0 ∪ Ŵ 1 ∪ Ŵ * (see Remark 2 for further details). Now, we map the domain into a reference bundle , where the computations are easier, free from undetermined constants, and are carried out once and for all. To this aim, for any x ∈ 1D , we introduce the map ψ x : γ x → γ d−1 between the generic fiber γ x and the reference fiber γ d−1 ⊂ R d−1 . Maps ψ x are instrumental to define the global map : → , where � = x∈� 1D {x} × γ d−1 denotes the reference computational domain (see Fig. 1 for an example of map ). Regularity assumptions are introduced on the maps ψ x and . In particular, we assume ψ x to be a C 1 −diffeomorphism, for all x ∈ 1D , and to be differentiable with respect to z (essentially to exclude any kinks along Ŵ * ).
We also demand that the supporting fiber 1D is preserved by map , so that the generic point z = (x, y) ∈ � is mapped into z = �(z) = ( x, y), with x ≡ x and y = ψ x (y) . Finally, without reducing the generality, we assume 1D to be the subset of with y = 0, i.e., 1D exactly coincides with the centerline of .
Remark 1 In a 2D setting, we may always select ψ x as a linear transformation, so that y = ψ x (y) = y/L(x), with L(x) = meas(γ x ). In 3D a similar choice is possible only for   specific configurations, for instance when is a cylindrical domain. In this case L(x) coincides with the diameter of the pipe along the centerline.

HiMod reduction
The HiMod technique has been proposed in [9,10] with the idea of exploiting the fiber structure demanded on , or, likewise, the preferential dynamics of the phenomenon at hand. Currently, three versions of HiMod reduction have been investigated, from both a theoretical and a numerical viewpoint (see [13] for a survey on the different approaches). Independently of the selected technique, the idea is to manage in a different way the dependence of the solution on the leading and on the transverse dynamics. In particular, since HiMod aims at providing enriched 1D models to be associated with the dominant direction, only the dominant dynamic is discretized via a standard finite element scheme, while getting information on the transverse dynamics via a modal expansion. In this section, we consider two of the available HiMod formulations.

Uniform HiMod reduction
The distinguishing feature of a uniform HiMod formulation is the adoption of a unique level of detail (i.e., the same number of modal functions) in modeling the transverse dynamics. For the sake of simplicity, we start from a steady setting. The function space associated with a uniform HiMod approach is where m ∈ N + is a given integer, V 1D ⊆ H 1 (� 1D ), and B = {ϕ j } j∈N + is a modal basis of functions in H 1 ( γ d−1 ), orthonormal with respect to the L 2 ( γ d−1 )-scalar product. The boundary conditions assigned on Ŵ 0 and Ŵ 1 are taken into account by the space V 1D , while the boundary data on Ŵ * are included in B. Space V m represents the hierarchy of models. We complete definition (3) by adding a conformity (V m ⊂ V ) and a spectral approximability (lim m→+∞ inf v m ∈V m �v − v m � V = 0, for any v ∈ V ) hypothesis on V m [9,10].

Remark 2
The analysis below is completely general with respect to the boundary data. So far the robustness of the HiMod reduction has been verified when either homogeneous Dirichlet or homogeneous Neumann boundary conditions are assigned on Ŵ 0 , Ŵ 1 , Ŵ * , or when non-homogeneous Dirichlet data are enforced on Ŵ 0 and Ŵ 1 . In general, the critical point is the identification of a basis B matching Robin boundary conditions or non homogeneous data on Ŵ * . A new strategy with respect to this issue has been recently proposed in [31].
With a view to unsteady problems, we introduce a time partition of the time window I into N subintervals I n = (t n−1 , t n ] of width k n = t n − t n−1 , for n = 1, . . . , N, with k = max n k n , t 0 ≡ 0 and t N ≡ T. This partition induces a subdivision of the cylinder Q into N space-time slabs S n = × I n , with n = 1, . . . , N. Notice that partition is not necessarily uniform, to match the possible time heterogeneities of the problem. Now, we look for an approximate solution to (2) coinciding, on each space-time slab S n , with a polynomial of degree at most q in time, with q ∈ N + , and with an element of V m in space, i.e., a function of the reduced space The boundary conditions in (2) according to the definition of Ŵ D , while functions ϕ j,r belong to the modal basis B. Moreover, since 0 � ∈ I 1 , the value v m (x, y, 0) has to be specified separately.
A priori, functions in V N m may exhibit a discontinuity at each time level, with continuity from the left. As a consequence, a different number of modal functions can be selected on each time interval I n (see Fig. 2). This choice leads to replace in (4) the modal index m with the index m n ∈ N + with n = 1, . . . , N. In such a case we adopt the term space-time slabwise uniform HiMod reduction and we change the notation in (4) is the vector that collects the number of modes used on each interval I n , with v m the generic function in V N m . The possible time discontinuity in V N m leads us to distinguish between the values v n,+ m = lim t→0 + v m (x, y, t n + t) and v n,− m = lim t→0 + v m (x, y, t n − t), and to define the temporal jump [v m ] n = v n,+ m − v n,− m at the generic time t n , for n = 0, . . . , N − 1. Notice that this jump is identically equal to zero for functions in V. This remark allows us to provide a weak formulation for problem (1) equivalent to (2): find u ∈ V such that (4)  where, for any w, ζ ∈ V , with u 0,+ = u 0,− = u 0 (x, y) and ∂Q n N = Ŵ N × I n for n = 1, . . . , N. The space-time slabwise uniform HiMod formulation can thus be stated: find u m ∈ V N m such that, for any v m ∈ V N m , The jump terms in (6) provide now an actual contribution, and we distinguish between the HiMod approximation u 0,− m ∈ V N m | I 1 of the initial datum u 0 and u 0,+ m that is unknown. The conformity and the spectral approximability hypotheses are now added slabwise to guarantee the well-posedness of formulation (8). Indeed, due to the discontinuity in time, we can only expect that V N m S n ⊂ V S n , while V N m � ⊂ V . Concerning the spatial discretization, following [9,10], we consider a finite element discretization of the function dependence on x, after introducing a subdivision, not necessarily uniform, of 1D into subintervals. The time discontinuity allows to employ a different 1D mesh on each space-time slab (see Fig. 2). In particular, we denote by T h n = {K n l } M n l=1 the spatial partition associated with S n for n = 1, . . . , N, with K n l = (x n l−1 , x n l ) the generic subinterval of width h n l = x n l − x n l−1 for l = 1, . . . , M n , with h n = max l h n l and x n 0 ≡ x 0 , x n M n ≡ x 1 . Then, we furnish each S n with the space X 1D,s h n of the conforming finite elements of degree s associated with T h n , and with dim(X 1D,s h n ) = N h n < +∞. A standard density hypothesis in V 1D is advanced on each finite element space. Thus, the discrete counterpart of formulation (8) (6) and (7) to V N m,h taking advantage of the slabwise splitting. By generalizing the notation used in [14][15][16] to denote finite elements that are continuous in space and discontinuous in time, we refer to V N m,h as to the HiMod c[M(m) G(s)]-dG(q) space (and, analogously, to (9) as to the c[M(m)G(s)]-dG(q) HiMod formulation). We mean that, on each S n , the full solution is replaced by a reduced solution continuous in space and discontinuous in time, obtained via a Galerkin approximation based on finite elements of degree s combined with the modal expansion associated with the multi-index m to discretize the space, and piecewise polynomials of degree q for the time discretization.
The finite element discretization along 1D allows us to further expand the Fourier coefficient v n,h j,r in (10)

HiMod versus PGD
Following the classification proposed in [5], both HiMod reduction and PGD can be categorized as a priori approaches since they do not rely on any solution to the problem at hand as, for instance, a POD strategy. Both the methods involve the weak form of the full problem and are based on a classical separation of variables. Nevertheless, while HiMod reduction applies this separation only to the space-time coordinates, PGD involves also problem parameters, such as boundary conditions or material properties, thus increasing the dimension of the space of the unknowns. HiMod applies a different discretization to the variables based on the physics of the problem. The accuracy for each variable may be tuned locally via a posteriori arguments. A PGD approach replaces in (3) the known modal function ϕ j (ψ x (y)), with y = (y, z) ′ , with a term of the form F j (y)G j (z) , with F j and G j unknowns. The two procedures both lead to 1D algebraic systems. In PGD they are intrinsically nonlinear and, in general, of large dimension. PGD requires therefore specific methods for the nonlinearity. In addition, the construction of the PGD approximation via a successive enrichment of an initial solution closely resembles the heuristic approach initially used in HiMod for selecting the number of transverse modes [10]. In this respect, the automatic selection of the HiMod approximation in [11] may represent an important evolution in a PGD setting as well. 1 To simplify notations, with the super-index h we understand both the space and time discretizations.

Pointwise HiMod reduction
A fixed number of modal functions on the whole may be too restrictive in the presence of spatial heterogeneities. This justifies the formalization of HiMod strategies alternative to the uniform approach, where a different number of modes is adopted in different subdomains of (via a piecewise HiMod reduction [10,11]), rather than in correspondence with each finite element node (thanks to a pointwise HiMod formulation [12]). We focus on the last approach. The numerical verification in [12] identifies the pointwise method as the best-performing one in the presence of either widespread or localized transverse dynamics.
The idea exploited in a pointwise HiMod expansion consists in rewriting (11) by emphasizing the sum on the finite element nodes, as and then in making the modal index m n dependent on the nodal index l. Space V N m,h is thus replaced by the new space where M n = [m n,1 , . . . , m n,N hn ] ′ ∈ N + N hn is the modal nodewise vector collecting the number of modes used at each finite element node of the slab S n for n = 1, . . . , N, whereas M is just the subindex used to denote a pointwise HiMod approximation. The nodewise tuning of the number of modes leads to an algebraic system sharing the same sparsity pattern as for the uniform case, but with a smaller dimension [12]. The formulation related to space V N M,h coincides with a space-time pointwise HiMod reduction and will be denoted by c[M(M)G(s)]-dG(q) form. It reads exactly as (9), simply by replacing space V N m,h with V N M,h . Notice that, since definition (12) strictly depends on the finite element discretization, there does not exist a weak counterpart of the pointwise HiMod formulation.

Uniform versus pointwise HiMod reduction: an example
We compare the uniform and the pointwise HiMod approaches on the steady test case 4 in [10], where the transport of oxygen in a wavy channel, representing a Bellhouse oxygenator for extra-corporeal circulation [32], is modeled. This problem is characterized by a widespread dynamics, that is suited to be reduced via both the HiMod techniques.   As far the HiMod reduction, we discretize the dependence of u on x via affine finite elements after introducing a partition of uniform step h = 0.1 on 1D . The transverse dynamics are described with a basis B of sinusoidal functions. To evaluate the integrals of the modal functions, we resort to Gaussian quadrature formulas based on four quadrature nodes per wavelength, at least. No stabilization scheme is used. We first apply the uniform HiMod approach by resorting to 11 modal functions (see Fig. 3, right). Indeed, as shown in [10], at least 11 modes are required to obtain a sufficiently reliable HiMod approximation.
As second assessment, we build the pointwise HiMod approximation u h M associated with the modal distribution M in Fig. 4, center. By comparing Fig. 3, right with Fig. 4, left we recognize that the two reduced solutions are very similar. In particular, the innermost contour lines associated with u h M are more accurate, despite the lower number of dof involved by the pointwse approximation (48,400 dof characterize u h 11 to be compared with 28,282 dof for u h M , see the corresponding sparsity pattern in Fig. 4, right). In accordance with [12], results in Figs. 3 and 4 show the improved modeling capabilities of the pointwise HiMod method vs the uniform approach, for a fixed computational effort. The main issue related to a pointwise formulation is the selection of the nodewise modal distribution. This corroborates the need for an automatic modal selection.

Adaptive HiMod reduction
Due to its significant impact on practical applications, we consider a goal-oriented framework (see, e.g., [17][18][19]), so that the predicted reduced model fits a goal functional representing a physical quantity of interest (e.g., mean or pointwise values, fluxes across sections or regions, the energy of the system, the vorticity of a turbulent flow). We denote by J the selected functional and we assume it is linear. We aim at approximating, within a prescribed tolerance TOL, the value J(u), with u solution to the full problem (2), via J (u h M ), where u h M is the reduced solution identified by a preprocessing phase. At this stage, we use a uniform and sufficiently fine discretization x n l , t n M n l=1 N n=1 on 1D × I so that we can neglect the error due to the space-time discretization.

The a posteriori modeling error analysis
We generalize the error analysis in [11] to an unsteady setting, to automatically produce the HiMod lookup diagram that provides the number of modes to be switched on at (see Fig. 7, left for an example). The a posteriori analysis is carried out on the slabwise uniform HiMod formulation, while the pointwise approximation u h M constitutes the output of the adaptive procedure in the next section.
According to a goal-oriented approach, we introduce the dual problem associated with where j is the density function associated with the goal functional J. Notice that, since V N m � ⊂ V , J has to be defined on V ∪ V N m and analogously for J cGdG . A null final condition, z N ,+ m = 0, allows to get rid of the first integral in (14), whereas boundary contributions may modify the definition of J cGdG when functional J involves a control on the boundary. The assignment of boundary conditions to the dual problem is a crucial issue that is usually tackled via the Lagrange identity. (6) can be alternatively rewritten integrating by parts the time derivative and recombining the jump terms as for any w, ζ ∈ V ∪ V N m . This form better fits the dual setting due to the reverse time scale.

Remark 3 The bilinear form
To derive the a posteriori modeling error estimator, we need to introduce the enriched primal and dual slabwise uniform HiMod problems, The analysis derived in [11] can be applied to the slabwise uniform HiMod formulations, to state the following m and e m + = u − u m + ∈ V ∪ V N m + be the modeling error associated with the reduced formulation (8) and (15), respectively for Thanks to the requirement on the dual final data, we have that J cGdG ≡ J. Result (18) identifies the estimator η mm + for the modeling error J (e m ) with the value |J (δu mm + )| , while guaranteeing the efficiency and the reliability of η mm + via the lower and upper bound, respectively. Following [11], to evaluate estimator η mm +, we can adopt three equivalent formulas, i.e., denote the weak primal and dual residual associated with formulation (8) and (13), respectively. Moreover, to make computable η mm +, we replace the reduced primal and dual solutions with corresponding discrete approximations. Estimator η mm + exhibits the structure typical of a hierarchical error estimator, yet in a goal-oriented framework. We refer to [11] for further computational remarks and for some considerations on hypothesis (17) that represents a generalization of the standard saturation assumption [34][35][36] to a goal-oriented setting.

Construction of the HiMod lookup diagram
Estimator η mm + is now used to automatically select the pointwise HiMod approximation u h M for problem (2) that guarantees the desired accuracy TOL on the functional error J (u − u h M ). To start the adaptive algorithm, we assign two initial (possibly small) values to the uniform modal indices m and m + . Then, we resort to the following five-stage procedure: (S1) we compute the discrete uniform reduced primal and dual solutions, u h m , u h m + , z m , z h m + , on the whole space-time cylinder Q; (S2) we evaluate the modeling estimator η n mm + = η mm + S n localized to each spacetime slab S n ; (S3) we apply the adaptive procedure outlined in Fig. 5 on each slab S n to predict the corresponding nodewise modal distribution M n , i.e., to build the HiMod lookup diagram (see below for all the details); (S4) we compute the discrete pointwise reduced primal and dual solutions,  η MM + ≤ TOL, the procedure stops, providing the HiMod lookup diagram in (S3) as final outcome. Vice versa, if η MM + > TOL, we come back to (S2).
Before detailing the adaptive procedure at stage (S3), some remarks are in order.
The computational effort associated with stage (S1) takes advantage of the time discontinuity of the c[M(M)G(s)]-dG(q) scheme. More sophisticated approaches such as checkpointing [37] may be adopted to further reduce the computational costs. The modeling estimator can obviously be evaluated in correspondence with any HiMod approximation [uniform as in (S2), slabwise uniform as in (19), pointwise as in (S5)]. Indeed, via the first definition in (19), it suffices to properly evaluate the bilinear form (6). Concerning the localization of the estimator to S n at stage (S2), by exploiting again the first definition in (19), we have Finally, the HiMod pointwise approximations u h M , z h M and u h M + , z h M + at stage (S4) are the solutions to problems (5), (13) and (15), (16) settled in the space V N M,h and V N M + ,h , respectively. In particular, we assume that V N M,h and V N M + ,h share the same spatial partitions T h n for n = 1, . . . , N, so that M + identifies reduced solutions with a pointwise larger number of modes with respect to u h M and z h M . Let us focus now on the adaptive procedure devised to commute the slabwise evaluations of η mm + into the lookup diagram predicted at stage (S3). We focus on the generic space-time slab S n and on the case of linear finite elements: (S3_1) we assign a number of modes equal to m to each node and to each subinterval of partition T h n ; (20) (S3_2) we evaluate the estimator η n,l mm + = η n mm + K n l localized to each interval K n l of T h n , for l = 1, . . . , M n ; (S3_3) we invoke an equidistribution criterion on the slabs as well as on the subintervals K n l . If η n,l mm + > TOL δ 1M /(N M n ), we increase by one the modal index associated with K n l (model refinement); if η n,l mm + <TOL δ 2M /(N M n ), we decrease by one such an index (model coarsening); otherwise, we preserve the current modal index; (S3_4) we update the number of modes associated with each finite element node by assigning to the generic node x n l , for l = 1, . . . , M n − 1, a number of modes equal to m n,l = min(m K n l , m K n l+1 ), with m K n l the number of modes assigned on the interval K n l . In particular, to avoid an abrupt variation of modes on consecutive nodes, the actual value m * n,l associated with x n l coincides with max(0.5 m n,l−1 + 0.5 m n,l+1 − 3, m n,l ) . The endpoints of 1D are updated separately as m n,0 = m K n 1 and m n,M n = m K n Mn if Dirichlet boundary conditions are not imposed on Ŵ 0 and on Ŵ 1 , respectively. The assignment of the modal indices m n,l predicts the modal multi-index M n = [m n,1 , . . . , m n,N hn ] ′ for the slab S n .
The procedure in (S3) is exemplified in Fig. 5 for a partition T h n of 1D consisting only of three subintervals K n l (l = 1, 2, 3). Of course, steps (S3_1)-(S3_4) are replayed on the enriched modal index m + , with a view to the evaluation of the modeling error estimator at stage (S5).
The adaptive modal algorithm includes both model refinement and coarsening. A minimum number of modes constrains the modal coarsening, while a maximum number of adaptive iterations is fixed to avoid too restrictive demands on TOL. The tuning parameters δ 1M and δ 2M at stage (S3_3) make the adaptive algorithm more robust, while increasing the corresponding computational efficiency. We set δ 1M = 0.5, δ 2M = 1.5. Finally, the modal update at step (S3_4) plays a crucial role since it explains how to build a pointwise HiMod approximation u h M starting from a HiMod lookup diagram, where the modes are associated with the subintervals.

Numerical verification
The numerical verification is carried out in a 2D setting. Moreover, to select the discrete HiMod space, we choose q = 0 and s = 1, i.e., we use linear finite elements to discretize the leading dynamics and functions piecewise constant in time. It can be checked that the adopted time discretization is equivalent to a modified backward Euler scheme [15].

Reliability of the adaptive HiMod reduction procedure
We approximate problem (1) on the rectangular domain � = (0, 3) × (0, 1) for t ∈ I = (0, 1), and by choosing Lu = − u + c · ∇u, with c = [10, 0] ′ . Besides the directionality induced by the advective field, we introduce a local heterogeneity via the source term f ≡ 10χ D , with χ D the characteristic function associated with the elliptic region D = {(x, y) : (x − 1.5) 2 + 4(y − 0.25) 2 ≤ 0.01}. Concerning the boundary conditions, homogeneous Dirichlet data are assigned on ∂�\Ŵ N , with Ŵ N = {(3, y) : 0 ≤ y ≤ 1}, where a homogeneous Neumann datum is enforced. Finally, a null initial datum u 0 is chosen. Figure 6, left shows at five different times, the contour plots of the full solution u approximated with FreeFem++ via a standard 2D cG(1)-dG(0) scheme on a uniform unstructured mesh of 10252 triangles. As expected, the convective field acts on the purely diffusive phenomenon by horizontally bending the contour lines. From a modeling viewpoint, we are simulating, for instance, the process of convection and diffusion of a pollutant emitted by a chimney localized at D, in the presence of a moderate horizontal wind. In this context, the full solution u(t) represents the pollutant concentration in the domain at a certain time t ∈ I We aim at controlling the mean value of the full solution on the whole at the final time T = 1, i.e., we select the goal functional J as J mean, The choice of a localized functional is challenging with a view to the modeling adaptive procedure. The dual problem is characterized by the differential operator L * z = − z − c · ∇z, with source term given by the density function j(x, y, t) = [meas(�)] −1 δ T associated with J mean,T , where δ T denotes the Dirac distribution associated with the final time. On Ŵ N a homogeneous Robin boundary condition is imposed, while a homogeneous Dirichlet data is assigned on ∂�\Ŵ N . A null final value z N ,+ m is selected. Both the primal and dual problems are computed by discretizing the supporting fiber Finally, the modeling tolerance TOL is set to 10 −2 , while the uniform modal indices m and m + are set to 1 and 3, respectively. The adaptive algorithm converges after 21 iterations and provides as output the HiMod lookup diagram in Fig. 7, left. The diagram coincides with the space-time rectangle 1D × I, where 1D and I exhibit the corresponding partition of uniform size h and k, respectively. A certain number of modal functions is associated with each cell K n l × k for l = 1, . . . , M n and n = 1, . . . , N. Thus, by resorting to the procedure in Fig. 5, (S3_4) it is possible to build the HiMod pointwise approximation u h M n for n = 1, . . . , N, i.e., the reduced solution u h M that guarantees the estimate |J mean,T (u) − J mean,T (u h M )| <TOL. The HiMod diagram in Fig. 7, left shows that few modes are demanded on the whole space-time domain, except for the two last time intervals, where a larger number of modes is switched on in correspondence with the localized source and the downstream region. More quantitative information are provided by the plot in Fig. 7, center of the number of modes associated with node x = 1.5 as a function of time. Only three modes are used on the whole time interval except for the subintervals I N −1 and I N , when five and 13 sine functions are required, respectively. The modal distribution predicted by the lookup diagram is completely coherent with a goal-oriented approach. Since we are interested in the mean value of the solution only at the final time, it is reasonable to expect a reliable approximation of the full solution only in correspondence with the last time intervals. This trend is confirmed by the corresponding pointwise HiMod approximation which reproduces more closely the full one during the last times of the simulation (compare Fig. 6, left and right).
In Fig. 7, right we show the value of η MM + on the same space-time structure of the HiMod diagram. The boxes associated with the largest values of the estimator identify a pattern similar to the one in Fig. 7, left.

Sensitivity of the adaptive HiMod reduction procedure to the goal-functional
We re-run the adaptive procedure by preserving all the input parameters, but for J = J  Fig. 8, right. In agreement with a goal-oriented approach, the reduced solution is far from the full one at T in Fig. 6, left-bottom. The mean value is controlled in an area where the full solution is extremely smooth so that a single mode is enough.

Robustness of the HiMod lookup diagram
The computational effort demanded by the adaptive procedure is justified by the possibility to employ the lookup diagram associated with a specific setting to hierarchically reduce a variant of this. Figure 9 performs this check on three variants of the test-case in Fig. 6. In more detail, we adopt the HiMod lookup diagram in Fig. 10, top-left to build the HiMod approximation for three new advection-diffusion problems, characterized by a different choice of the source term, namely, (Fig. 9, middle) and f 3 ≡ 10χ D 3 with (Fig. 9, bottom), respectively. Figure 9, left shows the HiMod approximations thus obtained. To check the reliability of the obtained solutions, we apply the HiMod adaptive algorithm directly to the three new problems. The corresponding lookup diagrams are gathered in Fig. 10, whereas the associated HiMod approximations are collected in Fig. 9, right. The matching between the contour plots in the two columns of Fig. 9 is substantial, in particular for the two-source test case which represents the most sigificant variant with respect to reference configuration.
As last investigation, we check the robustness of the HiMod lookup diagram with respect to the shape of the computational domain. This is a challenging issue due to the crucial role played by the maps ψ x and in the HiMod reduction. For this reason, we focus on a steady setting, in particular on the one in Fig. 3. We consider two variants of the wavy-shaped domain. With the first one, we simply reduce the height of the sinusoidal sections (see Fig. 11, top), whereas in the second variant we add a rectangular channel at the beginning and at the end of the original geometry (see Fig. 11, bottom). Figure 11 compares the solution computed on the modal distribution in Fig. 4, center with the HiMod approximation provided by the adaptive HiMod procedure. As expected, the matching between the two approximations is very good in Fig. 11, top. Despite the abrupt change in the domain shape, we get acceptable results also for the second geometric variant, thus confirming the good robustness of the HiMod lookup diagram with respect to possible changes of the original problem.

Combined HiMod reduction and space-time adaptation
Goal of this section is to enrich the information provided by the HiMod lookup diagram by predicting also the space-time partition of 1D × I. Consequently, we remove any assumption on the finite element discretization as well as on the time partition {I n }. In practice, we expect to replace the diagram in Fig. 7, left with a new diagram characterized by a non uniform horizontal (spatial) and vertical (temporal) spacing.

The a posteriori estimator for the global error
With a view to a global adaptation, following, e.g., [11,[20][21][22], we derive an a posteriori estimator for the global error E h m = e m + e h m , where the contributions of the modeling (e m = u − u m ) and of the discretization (e h m = u m − u h m ) errors remain distinct. In particular, since we are interested also in an adaptive selection of the space and time step size, we expect that the estimator for the discretization error consists of a spatial contribution separate from the temporal one [23][24][25][26].
As for the adaptive HiMod reduction, we carry out the a posteriori analysis in a slabwise uniform HiMod setting. The following statement plays a crucial role in the definition of the global error estimator.    contribution takes into account the error associated with both the spatial and the temporal discretizations. The main effort of this section will be to explicitly estimate this term, with the additional requirement of distinguishing the space from the time contribution. As in [11], we modify the standard goal-oriented analysis to tackle the intrinsic dimensionally hybrid nature of a HiMod reduced formulation.
Concerning hypothesis (22), it essentially coincides with a sufficient grid resolution requirement since establishing a ratio between the modeling and the discretization errors. With a view to estimate |J (e h m )|, we preliminarily prove the following Galerkin orthogonality property for the discretization error e h m . {x} × γ x , with K n l the generic subinterval of T h n , while we denote the interface between R n τ and R n τ +1 by ζ n τ , for τ = 1, . . . , M n − 1 and n = 1, . . . , N, and with ζ n 0 ≡ Ŵ 0 and ζ n M n ≡ Ŵ 1 . Finally, S R n l = R n l × I n denotes the space-time prism associated with R n l , while L R n l = ∂R n l × I n identifies the corresponding lateral surface. We introduce now the spatial and temporal local residuals. For a fixed time interval I n and for any R n l , we consider the internal residual and the boundary residual . By definition, the projection error v − T n v is orthogonal to any function c constant in time, so that whereas the estimate can be proved [16]. Notice that no constant is involved in this result. Concerning the Clément quasi-interpolant, the estimates hold, for any v ∈ H 1 (� 1D ), where K denotes a generic interval of 1D , K is the associated patch of elements, and with C 1 and C 2 constants depending on the relative size of the elements constituting K [38].

Lemma 1 For any v h
We are now ready to prove the following result:  (34) and (35), on q and on max n m n , where the residuals are defined by with r R n l = T n r R n l , j R n l = T n j R n l , h n l and k n the length of the generic subinterval K n l and I n , respectively for l = 1, . . . , M n and n = 1, . . . , N, and with δ 1,n the Kronecker symbol associated with the first slab S 1 , while the weights are given by with the patch associated with the subinterval K n l , L(x) = meas(γ x ), z n j,r and z n,h j,r the modal coefficients associated with the dual solution z m and with the corresponding discretization z h m , respectively.
Proof We start from the dual problem (13) by choosing v m = e h m and we apply the orthogonality relation (25). It follows that, for any v h m ∈ V N m,h , Thanks to definitions (29) and (30), we have We consider separately the four terms (I)-(IV). In particular, we choose v h m coinciding with z h m + T n (I 1 (z m − z h m )), with z h m the discrete HiMod approximation of the dual . solution. In particular, the Clément operator involves only the x-dependent modal coefficients since it is one-dimensional. Notice that, since we estimate slabwise the terms (I)-(IV), all the functions in V N m and V N m,h have to be meant restricted to I n , for each n = 1, . . . , N. Function v h m is extended to zero outside I n when considered as a function of V N m,h . To exploit the projection and the interpolation estimates in (33)-(35), we consider the following splitting Let us focus on term (I). Using the splitting above, the definition of the averaged residual r R n l and of the projection operation T n , and by combining results (32) and (33) with the Cauchy-Schwarz inequality, we obtain We now consider separately the norm associated with the interpolation error. Let w m be a generic element in V N m . By exploiting the modal expansion for w m and the orthonormality of the modal basis, together with interpolation estimate (34), we obtain (40) |� w n j,r | 2 where D(x, ψ −1 x ( y)) = L(x) −1 denotes the Jacobian associated with the map ψ x , and with γ 1 the reference fiber for the two-dimensional setting. Via this estimate, we obtain the following bound for the term (I) in (38): with C a constant depending on C 1 in (34), q and m n . From now on, C denotes a constant whose value may change from line to line. Term (II) can be bounded analogously to contribution (I), by restricting the computations on the lateral surface L R n l of S R n l . This yields Inequality (40) is replaced by a corresponding trace estimate, obtained essentially by invoking result (35) instead of (34), to have for any w m ∈ V N m . Combining this result with (42), we attain the following control for the second term in (38): where constant C depends on C 2 in (35), q and m n . We focus now on term (III) and, first of all, we apply again splitting (39): Now, thanks to the mean value theorem, we remark that, for any function w m ∈ V n m , (43) �w m − I 1 (w m )� 2 � w n j,r � 2 with t * n ∈ (t n−1 , t n ), as well as equality �J n−1 � L 2 (S R n l ) = k 1 2 n �J n−1 � L 2 (R n l ) trivially holds. Moving from these results and by exploiting the definition of the projection operator T n , the Cauchy-Schwarz inequality and estimate (40), we derive the final bound for (III): with C as in (41). The last term in (38) can be controlled by repeating the same computations adopted for (III), by replacing the temporal residual J n−1 with the initial error e h,0,− m and by focusing on the first time interval. We achieve the following estimate with C as in (43). Now, result (36) follows by properly combining the individual estimates obtained for terms (I)-(IV).
Moving from (36), we propose as error estimator for the discretization contribution in (24) the value so that the estimator for the global functional error, |J (ε h m )|, coincides with η h mm + = η mm + + η h , with η mm + as in (19). In particular, since it is straightforward to distinguish in η h the space from the time contribution given by respectively, it is immediate to decompose η h mm + into a modeling, a space and a time contribution, as This splitting will be crucial with a view to the global adaptive procedure. Both the estimators η h S and η h T share the structure characterizing a goal-oriented analysis, i.e., they coincide with the product of a residual depending on the primal solution and a weight related to the dual solution. In addition, we remark that, due to the HiMod procedure, the contribution along the x-and y-direction in the weights is split.
Some computational remarks on estimator η h are now in order.
To make computable the weights, we replace the dual solution z m with a computable discrete counterpart z * ,h m . A possibility is to resort to the discrete enriched dual solution z h m + . Nevertheless, since the temporal weights involve the time derivative of z m , we resort to a temporal recovery procedure yielding an approximation z * ,h m that is at least linear in time. In particular, we follow the approach in [25,26]. The dependence of the weights on the dual discretization error rather than on the dual solution is optimal in terms of convergence. Moreover, the time averaged residuals r R n l and j R n l make the estimator more reliable since �w� L 2 (I n ) ≤ �w� L 2 (I n ) as well as �w − w� L 2 (I n ) ≤ �w� L 2 (I n ) , for any function w ∈ L 2 (I n ). An extra care has to be devoted to the computation of the temporal residual J n−1 that combines solutions associated with two different meshes. We use an interpolation operator from the degrees of freedom of T h n onto the ones associated with T h n+1 . Finally, the analysis in Proposition 3 may be generalized to a 3D framework provided that map ψ x is properly chosen. In particular, the orthonormality of basis B may be exploited to derive estimates (40) and (43) only if D −1 (x, ψ −1 x ( y)) does not depend on y. This has to be explicitly demanded in a 3D setting while it always holds in a 2D framework.

Building the space-time HiMod lookup diagram
Goal of this section is to keep the global functional error below a fixed tolerance TOL via an automatic selection of the modal distribution and now also of the space-time mesh K n l , I n . Different strategies are followed in the literature to combine model with mesh adaptation [11,20,21,39]. The approach we propose iteratively alternates model with spacetime mesh adaptation, by advantageously exploiting the additive structure of the global error estimator (46). For this reason, we distinguish a model (TOL_MODEL) and a mesh (TOL_MESH) tolerances, such that TOL_MODEL + TOL_MESH = TOL. Then, we follow the procedure outlined in Fig. 12. We distinguish two main modules, ADMOD devoted to model adaptation and ADMESH dealing with the space-time mesh adaptation. The module ADMOD exactly implements the five-stage adaptive procedure (S1)-(S5). Concerning the space-time mesh adaptation, the algorithm set by ADMESH is very straightforward, due to the one-dimensional nature of both the spatial and temporal meshes. In particular, while the space adaptation includes both mesh refinement (via bisection) and coarsening (gluing two consecutive intervals where η h S is below tolerance), the time adaptive algorithm deals only with mesh refinement. This suggests to start the adaptive procedure on a sufficiently coarse time partition. Error equidistribution drives both the space and time adaptation. A maximum value constrains the number of iterations as well as tuning parameters δ 1H (=0.5), δ 2H (=1.5) limit the spatial mesh refinement and coarsening to the worst and to the best subintervals, respectively.
After a preliminary check on the accuracy of the global error estimator associated with the initial uniform modal distribution and the initial uniform space-time grid, model adaptation takes place till the accuracy TOL_MODEL is met by estimator η MM +. Then, we check if model adaptation suffices to provide the global tolerance TOL without any space-time mesh adaptation. If not the module ADMESH is activated. In particular, we apply the spatial rather than temporal adaptation depending on which of the estimators η h S , η h T is the greatest one. When η h < TOL_MESH, we come back to the initial check on the global accuracy.
A maximum number of iterations ensures the end of the whole adaptive procedure. We remark that each time the space-time partition is updated, a projection of the primal and dual solutions involved in the evaluation of the error estimator is demanded. As for the choice of the tolerances, we resort to a convex combination of the two tolerances, by selecting TOL_MODEL = θ TOL and TOL_MESH = (1 − θ) TOL, with 0 ≤ θ ≤ 1 [11]. The parameter θ settles a relation between model and discretization error, in accordance with requirement (22).
Finally, we refer to the outcome of the whole adaptive algorithm as to the space-time HiMod lookup diagram. Some instances of this table are provided in the next section.

Numerical verification
In this section we assess the reliability of the global adaptive procedure.

Reliability of the space-time adaptive HiMod reduction procedure
The test case used to validate the modeling adaptive procedure for J = J mean,T is now tackled by activating the mesh adaptation as well. We preserve the same values of the previous run for TOL, for the initial uniform modal indices m and m + , and for the initial space-time mesh. Then, we set θ = 0.5.
The adaptive procedure converges after 50 iterations, with 23 model iterations followed by nine and eight adaptations of the spatial and of the temporal mesh, respectively and by ten additional model adaptations. The final outcome of the adaptive procedure is the HiMod lookup diagram in Fig. 13, top-left. A comparison between this table and the one in Fig. 7, left shows a similar trend for the modes, i.e., a gradual increment of the number of modes as we approach the final time and in correspondence with the source location and the downstream areas. Nevertheless, the combination of model with mesh adaptation reduces from 3 to 1 the number of modes used in the first phase of the test case (compare Fig. 7, center with Fig. 14, left). Concerning the spatial adaptation, a coarse mesh consisting of less than 20 subintervals and refined around x = 1.5 is predicted for the first time intervals. Then, this number increases with an abrupt variation in the last time interval when it reaches its maximum (see Fig. 14, center). The monotone trend characterizing the model and the spatial mesh adaptation is qualitatively the same, exhibiting a refinement of the modes and of the finite element partition confined to the last time intervals, in accordance with the goal quantity.
On the contrary, the time adaptation yields a non monotone prediction for the time step distribution, as depicted in Fig. 14, right. Essentially we recognize two phases when the initial time step is considerably reduced, the first one around the initial time and the second one just before time T. A strong refinement of the initial grid is recurrent in mesh adaptation and here it likely balances the initial rough modal and spatial discretizations. The second refinement occurs when the control of the mean value becomes more relevant. At time t = 0.8, both the modal discretization and the space-time mesh are considerably refined to ensure the imposed tolerance. Probably, a complex interplay among the three discretizations takes place during the last time intervals, so that the severe demand on the time step can be then relaxed before reaching the final time. Figure 13 gathers the distribution of the three error estimators on the space-time lookup diagram. The choice made for the tolerances leads to values of the same order of magnitude for η MM + and η h S , while the error estimator associated with the time discretization assumes larger values.
As shown in Fig. 15, the c[M(M)G(1)]-dG(0) HiMod solution associated with the diagram in Fig. 13, top-left is qualitatively different from the one in Fig. 6, right. The adoption of a single mode till t = 0.7 identifies a reduced solution which is initially very far from the full one. Nevertheless, the time steps predicted by the adaptive algorithm are enough to refine, during the last time intervals, the number of modes as well as the partition along 1D so that solution u h M becomes fully comparable with the full one at the final time.

Assignment of Neumann boundary conditions
We challenge the whole adaptive procedure by modifying the boundary conditions in the previous test case. We assign a homogeneous Neumann condition on the whole boundary, except for the edge Ŵ D = {(0, y) : 0 ≤ y ≤ 1} where we preserve the homogeneous Dirichlet data. The new condition along the horizontal sides leads to select a new modal basis. After identifying the reference fiber γ 1 with the interval [0, 1], we choose B = {ϕ j ( y) = √ 2 cos(π j y)} j∈N . Figure 16, left shows the cG(1)-dG(0) full solution at four different times, computed with FreeFem++ on a uniform unstructured mesh of 10,252 elements. In particular, the new flux-free configuration erases the horizontal dynamics in Fig. 6, pushing the pollutant to contaminate also the northeast and the southeast areas. If we set the global adaptive procedure to control J mean,T , we do not expect much benefit from the modal basis since all the cosine functions have a null mean except ϕ 0 . A completely different prediction is performed by selecting the goal functional The global tolerance TOL = 10 −2 is now guaranteed after 30 model iterations, followed by seven spatial and nine temporal mesh adaptations, plus a final model adaptation. The space-time adaptive HiMod lookup diagram yielded by the adaptive procedure is shown in Fig. 18, top-right. The number of cosine functions is gradually increased to eight in correspondence with D. Additional modes are now demanded also upstream the source location in contrast to Fig. 13, top-left. The modal as well as the spatial mesh cardinality trend is very similar to the one in Fig. 14, whereas three refinements of the time step now occur (see Fig. 17, bottom). The additional refinement about in the middle of the time window corresponds to the phase when the pointwise HiMod solution starts to become similar to the full one. Indeed, as shown in Fig. 16, right solution u h M is initially far from the full one (and similar to the approximation in Fig. 18). Then, from t = 0.5, u h M becomes more and more similar to the full solution till, at the final time, the two solutions are almost identical.

Robustness of the space-time HiMod lookup diagram
We replicate the test performed for model adaptation, by checking the robustness of the space-time HiMod diagram with respect to possible variants of the reference problem. To this aim, we employ the space-time lookup diagram in Fig. 19, left tailored on the problem in Fig. 6 to build the HiMod approximation for the problem identified by the source term f 3 in Fig. 9, bottom. Figure 20 compares the approximation thus obtained (left) with the HiMod approximation yielded by the global adaptive procedure (right), whose space-time HiMod lookup diagram is provided in Fig. 19, right. The agreement between the two solutions is satisfying. The adaptive procedure optimizes the computational costs by predicting a lower number of modes and a coarser mesh. Nevertheless, the possibility of exploiting a previously computed HiMod diagram thus avoiding the cost of the adaptive procedure is a sufficient motivation to exploit the precomputed diagram.

Computational saving
Goal of this section is to verify the benefits due to the HiMod adaptation procedure in terms of CPU times 2 with respect to a full and a uniform HiMod approximation. For the sake of simplicity, we consider a steady problem. We solve on the rectangular domain � = (0, 2π) × (0, π) the advection-diffusion problem − u + c · ∇u = f , with c = (10, 0) ′ , by assigning homogeneous Dirichlet data on ∂�\Ŵ N , with Ŵ N = {(2π , y) : 0 ≤ y ≤ π}, and a homogeneous Neumann data on Ŵ N . Then, we choose the source term such that the exact solution coincides with u(x, y) = sin y sin 0.01y(x 3 − 12π 2 x) (see Fig. 21, top-left). We first investigate the advantages due to a uniform HiMod reduction with respect to a standard 2D finite element approximation. We fix a number of dof around 190 and we compute the L 2 (�)norm of the error associated with the full approximation and with the uniform HiMod solution based on 17 modes and a uniform subdivision of the supporting fiber into 11 subintervals (see Fig. 21, top-right). As Table 1 shows, we gain an order of accuracy via the HiMod reduction. A comparison in terms of CPU time is not reasonable in such a case since the HiMod code is not yet optimized. By resorting to a modal adaptivity and for a comparable number of dof, we obtain a HiMod approximation more accurate with respect to the uniform one (compare the contour plots in Fig. 21, top-right and bottomleft and the values in Table 1) with a similar CPU time (in s). The modal distribution yielded by the adaptive procedure is shown in Fig. 22, left. A number of modes less than 2 All the experiments have been performed using Matlab 2011a 64-bit on a Lenovo ThinkPad T430 equipped with a Intel Core i5 3320M 2x 2.6-3.3 GHz processor and 4 GB of RAM. 17 is demanded on the whole domain except for the last three nodes. Concerning the CPU time, we quantify only the seconds demanded to build the HiMod approximation from the predicted modal distribution, since we have verified the robustness of the HiMod diagrams.
We now add the adaptivity of the spatial mesh. Table 2 compares the accuracy of a full with a HiMod approximation for about the same number of dof. We adopt two different tolerances to drive the global adaptive procedure. The corresponding HiMod diagrams are shown in Fig. 22, center and right. The accuracy characterizing the adapted HiMod solution is higher in both the cases and the computational times remain contained. Figure 21, bottom-right shows the HiMod approximation characterized by 622 dof. The maximum number of modes predicted by the adaptive procedure is still 17 but it is evident that the employment of an adapted mesh improves the reliability of the reduced solution as it is qualitatively evident by comparing the two contourplots in Fig. 21, bottom.

Validation of the HiMod reduction
This is a first attempt of validation for the HiMod reduction procedure. For this purpose, we focus on the experimental and modeling analysis provided in [27] dealing with a reactive transport in homogeneous porous media.
We consider the experimental setting outlined in Fig. 23. It consists of a rectangular laboratory flow cell of dimension 2.5 dm ×1 dm ×0.08 dm along the x-, y-and z-direction, respectively. The cell is filled with a porous media with measured porosity equal to 0.375 and it is initially saturated with an aqueous solution. Segment Ŵ inlet = {(0, y, z) : 0.5 ≤ y ≤ 1, 0 ≤ z ≤ 0.08} coincides with an inlet boundary, where a constant concentration, modeling the injection of a reactive component, is assigned. Simultaneously, a flow rate of 12 ml/h is set at the outlet Ŵ outlet = {(2.5, y, z) : 0 ≤ y ≤ 1, 0 ≤ z ≤ 0.08}, resulting in an average water velocity of about 0.404 dm/h at the equilibrium. We remark that the set-up of the experiment is designed to have a pseudo-1D flow, parallel to the x-axis. Finally, ten sampling ports are located in the cell, to collect measurements of the reactive fluid concentration. Sampling is performed four times during each experiment. The concentration measurements represent the data we aim at matching via a HiMod reduced modeling in the same spirit of the analysis in [27]. The reactive transport experiment is conducted for 60 h, though a stationary state is reached already after 15 h from the beginning of the experiment, so that we restrict the time window of investigation to (0, 30).
For all the further experimental data we refer to [27] since a greater level of detail on the experimental setting is beyond the purposes of the paper.
From a modeling viewpoint, since the setting is invariant along the z-axis, we can simulate the experiment in an effective way as a two-dimensional flow. In particular, we adopt the unsteady equation with � = (0, 2.5) × (0, 1), to model the process of advection and diffusion of the reactive component. Notice that (47) represents a simplified version of the original model in [27]. A preliminary tuning of the model parameters has been carried out to make the solution of the two models as close as possible in the considered experimental context. In more detail, we adopt a constant diffusive coefficient whose value is set, via a trial and error procedure, to replicate the action of the diffusive tensor used in [27]. Moreover, following [27], we select the value for the flux velocity by solving an additional Darcy problem. Figure 24, left shows the full solution computed with FreeFem++ on a uniform unstructured mesh of 13,078 triangles at t = 5, 11, 15, 19 h. The reactive fluid gradually spreads into the flow cell.
We now test the HiMod reduction procedure. We first resort to a uniform HiMod approximation and we use 20 modal functions to describe the transverse dynamics. We adopt a uniform space-time discretization along 1D and (0, 30), with step h = 0.05 and k = 0.5, respectively. In Fig. 24  reduction of the (spatial) dof (1000 vs 13,078). Now, we focus on the actual validation phase. For this purpose, in Fig. 25, we compare the measured (circle symbols) with the simulated concentrations (diamond symbols for the uniform HiMod approximation and star symbols for the 2D finite element discretization) in correspondence with eight of the ten sampling ports in Fig. 23. We refer only to one of the two sets of data available in [27]. Qualitatively, at each port, we recognize a first phase of about 8 h when the chemical breakthrough, characterized by a sigmoid shape curve, occurs; successively, the steady state is reached and each curve exhibits a plateau. The agreement between simulated and measured concentrations is good and comparable with the one of Fig. 3 in [27]. In particular, the full approximation improves the tracking of the data in correspondence with port A4. On the contrary, the prediction at ports A3, C3, D3 is more reliable when resorting to the HiMod approximation, despite the reduced number of dof.
As last test, we assess the reliability of the modeling adaptive procedure in a validation context. We aim at evaluating the reactive fluid concentration at t = 15 h via the c[M(M)G(1)]-dG(0) HiMod solution predicted by the modeling adaptive procedure. We consequently choose functional J as J 15 (ζ ) = [meas(�)] −1 � ζ(x, y, 15) d�. The expectation is to obtain a value for the concentration similar to the one provided by u h 20 and not so far from the experimental data. We set the adaptive algorithm with TOL = 10 −3 , m = 1, m + = 3. Concerning the space-time discretization, we fix a uniform space-time subdivision of 1D × I, with h = 0.05 and k = 0.5. Finally, we reduce the time window to (0, 15) due to the stationary regime of the flow in the interval (15,30). The modeling adaptive algorithm converges after 599 iterations and provides the HiMod lookup diagram in Fig. 26, left characterized by the space-time distribution of η MM + in Fig. 26, right. Both the diagrams corroborate the complexity of this experiment. In contrast to a more localized phenomenon such as the convection-diffusion of a pollutant in the previous sections, the refinement of the number of modes now gradually involves the whole 1D as we approach time t . The non uniform trend of the estimator highlights the demanding work performed by the adaptive procedure to guarantee tolerance TOL. Despite these difficulties, the maximum number of modal functions required by the lookup diagram is 12 to be associated with the area closer to the inlet and with the time intervals immediately preceding the steady state. The pointwise HiMod approximation u h M generated by the online phase is depicted in Fig. 27, for t = 5, 7, 11, 15 h. The trend of the adapted solution becomes more and more similar to the one in Fig. 24, as t approaches t .
Finally, we examine the concentration values predicted by the adapted HiMod solution at t = 15 h in correspondence with the eight ports in Fig. 25 (see the square symbols). It is evident the good matching of the simulated concentrations between u h 20 and u h M , with a slight different prediction at ports C3 and D3.