A hyperreduced manifold learning approach to nonlinear model order reduction for the homogenisation of hyperelastic RVEs

TL;DR

This paper introduces a hyperreduced manifold learning approach for nonlinear model order reduction in hyperelastic RVE homogenisation, significantly accelerating computations while maintaining accuracy.

Contribution

It integrates hyperreduction methods into a graph-based nonlinear MOR framework, achieving complexity independent of system size and improving online linearisation for robustness.

Findings

Accelerates RVE computations by over two orders of magnitude.

Maintains high accuracy with minimal training data.

Outperforms alternative methods in accuracy-runtime trade-off.

Abstract

In a recent work, we proposed a graph-based manifold learning scheme for the nonlinear Galerkin-reduction of quasi-static solid mechanical problems [1]. The resulting nonlinear approximation spaces can closely and flexibly represent nonlinear solution manifolds. The present work discusses how this nonlinear model order reduction (MOR) approach can be employed to reduce online computational costs by multiple orders of magnitude while retaining high levels of accuracy. We integrate two popular hyperreduction methods into the nonlinear MOR framework and discuss how we achieve an algorithmic complexity which is independent from the original system size. Furthermore, improvements are made to the local online linearisation scheme for the sake of performance and robustness. On an example RVE problem, the MOR scheme accelerates computations by more than two orders of magnitude with little training data and negligible loss of accuracy. Additionally, the algorithm Pareto-dominates alternative approaches in the trade-off between accuracy and runtime on the considered example.

Tables

Table 1. Table 1 : Parameters for the numerical experiments.

parameter	variable	value
Young’s modulus matrix	$E$	1000 Nmm-2
Young’s modulus inclusions	$E$	3000 Nmm-2
Poisson’s ratio	$\nu$	0.2
RVE edge length		6 mm
inclusion radius		1.5 mm
centre inclusion 1 ($x,y,z$)		(2,2,2) mm
centre inclusion 2 ($x,y,z$)		(4,4,4) mm
snapshot number	$s$	200
validation set size		500
step size	$\Delta H_{\text{LP}}$	0.03
perturbation size	$\Delta H_{\text{LS}}$	0.015
independent DOFs	$D$	19,182
element number	$n_{\text{elem}}$	4,821
node number	$n_{\text{node}}$	7,650

Table 2. Table 2 : POD DEIM mean error. Failed simulations highlighted in dark gray, large mean errors highlighted in light gray.

		$m$
		20	25	30	35	40	50	75	100	150	200	250
$d$	9	$\times$	$\times$	0.6360	0.5264	0.4964	$\times$	$\times$	4.9010
	12		$\times$	$\times$	0.4414	2.0346	1.5257	1.8905	1.1872
	15			$\times$	$\times$	$\times$	0.3525	0.2949	0.2922	0.2883	0.2846
	20				$\times$	$\times$	$\times$	$\times$	4.3245	$\times$	3.5516
	30					$\times$	$\times$	$\times$	$\times$	0.1158	0.0923	$\times$

Table 3. Table 3 : POD LEHM mean error. Failed simulations highlighted in dark gray, large mean errors highlighted in light gray.

		$m$
		20	25	30	35	40	50	75	100	150	200	250
$d$	9	$\times$	0.7228	0.5178	0.5018	0.4785	$\times$	3.1818	1.9488
	12		$\times$	0.4452	$\times$	1.1310	1.0947	1.0243	0.8432
	15			$\times$	$\times$	$\times$	$\times$	0.2887	0.2847	0.2820
	20				$\times$	$\times$	$\times$	0.7315	1.9190	2.9067	2.4180
	30					$\times$	$\times$	$\times$	0.1075	0.0852	0.0824	$\times$

Table 4. Table 4 : PM LEHM mean error. Failed simulations highlighted in dark gray, large mean errors highlighted in light gray.

		$m$
		20	25	30	35	40	50	75	100	150	200	250
$d$	9	$\times$	0.3858	$\times$	$\times$	$\times$	$\times$	$\times$	$\times$
	12		$\times$	0.3257	0.2852	0.2602	0.2487	$\times$	$\times$
	15			0.2287	0.2228	0.2204	$\times$	0.1834	0.1784	2.2028	1.5465
	20				$\times$	$\times$	$\times$	0.8375	0.7812	0.7827	1.4868
	30					$\times$	$\times$	$\times$	0.0628	0.0584	2.8437	$\times$

Table 5. Table 5 : LLE LEHM mean error. Failed simulations highlighted in dark gray, large mean errors highlighted in light gray.

		$m$
		20	25	30	35	40	50	75	100	150	200	250
$d$	9	0.2519	0.2254	0.2229	0.2193	0.2135	0.2090	1.1803	1.4237
	12		$\times$	$\times$	$\times$	0.1789	0.1635	2.6949	2.6363
	15			$\times$	0.2285	0.1804	0.1454	0.1309	0.1309	0.5261	2.3642
	20				$\times$	$\times$	0.1575	0.1111	0.1018	0.0989	2.4503
	30					$\times$	$\times$	0.0946	0.0599	4.3670	3.4263	3.4195

Table 6. Table 6 : POD LSPG mean error. Failed simulations highlighted in dark gray, large mean errors highlighted in light gray.

		$m$
		20	25	30	35	40	50	75	100	150	200	250
$d$	9	0.7342	0.6632	0.6190	0.5944	0.5860	0.9419	0.7918	0.6996
	12		0.5502	0.4899	0.4731	0.6901	0.6401	0.5724	0.5320
	15			0.4398	0.4425	0.3963	0.3965	0.3923	0.3881	0.3850	0.3786
	20				1.3323	0.8790	0.6298	0.5070	0.7100	1.0957	0.8303
	30					0.2302	0.1913	0.1517	0.1429	0.1356	0.1323	1.0020

Table 7. Table 7 : LPOD LSPG mean error. Failed simulations highlighted in dark gray, large mean errors highlighted in light gray.

		$m$
		20	25	30	35	40	50	75	100	150	200	250
$d$	9	0.3462	0.3374	0.3376	0.3307	0.3271	0.3258	0.3191	0.3141
	12		0.2854	0.2769	0.2748	0.2699	0.2672	0.5859	0.4999
	15			2.4800	2.0267	1.6828	1.3331	1.2201	1.0216	0.7958	0.6923
	20				0.2315	0.2169	0.2040	0.1933	0.1928	0.6251	0.4960
	30					0.2230	0.1565	0.1263	1.4651	1.0175	1.0903	0.8995

Table 8. Table 8 : PM LSPG mean error. Failed simulations highlighted in dark gray, large mean errors highlighted in light gray.

		$m$
		20	25	30	35	40	50	75	100	150	200	250
$d$	9	0.5979	$\times$	$\times$	$\times$	$\times$	0.7884	$\times$	$\times$
	12		0.3419	0.3410	0.3262	0.3114	0.3123	$\times$	0.5107
	15			0.2693	0.2603	0.2564	0.2559	0.2395	0.2346	0.4214	0.3744
	20				0.2187	0.2216	0.2087	0.4496	0.3558	0.2858	0.5551
	30					0.1726	0.1279	0.1055	0.0980	0.0934	0.4669	0.6089

Table 9. Table 9 : LLE LSPG mean error. Failed simulations highlighted in dark gray, large mean errors highlighted in light gray.

		$m$
		20	25	30	35	40	50	75	100	150	200	250
$d$	9	0.2543	0.2503	0.2504	0.2508	0.2504	0.2496	0.4951	0.5790
	12		0.2141	0.2086	0.2073	0.2072	0.4544	0.4452	0.3916
	15			0.1884	0.1786	0.1726	0.1692	0.1645	0.1630	0.2552	0.5629
	20				0.1670	0.1577	0.1501	0.1401	0.1365	0.1320	0.6165
	30					0.1469	0.1266	0.1007	0.0873	0.6937	0.5472	0.4516

Equations

$${\int }_{\Omega }\mathbf{P}:\delta \mathbf{F}dV-{\int }_{\Omega }\mathbf{b}\cdot \delta \text{u}dV-{\int }_{\partial \Omega }\mathbf{t}\cdot \delta \text{u}dA=0,$$

$$W=\frac{\mu }{2}(I_c-3)+\frac{\kappa }{4}(J^2-1-2\ln J),$$

$$\sum _{e\in E}\delta {\mathbf{u}}^{ekT}{\mathbf{g}}^{ek}=0,\text{with}{\mathbf{g}}^{ek}={\int }_{{\Omega }^e}\mathbf{P}\frac{\partial N^{ek}}{\partial \mathbf{X}}dV-{\int }_{{\Omega }^e}\mathbf{b}{N}^{ek}dV-{\int }_{\partial {\Omega }^e}\mathbf{t}{N}^{ek}dA,$$

$$g_{3K+\alpha }=\sum _{e\in E}{g}_{\alpha }^{ek},$$

$$\delta {\vec{u}}^T\vec{g}(\vec{u};\vec{p})=0,$$

$$\vec{g}(\vec{u};\vec{p})=0.$$

$$\mathbf{K}({\mathbf{u}}_{\text{cur}};\vec{p})\Delta \vec{u}=-\vec{g}({\mathbf{u}}_{\text{cur}};\vec{p}).$$

$$\mathbf{K}({\mathbf{u}}_{\text{cur}};\vec{p})=\frac{\partial \vec{g}({\mathbf{u}}_{\text{cur}};\vec{p})}{\partial \vec{u}},$$

$$K_{3K+\alpha ,3L+\beta }=\sum _{e\in E}{K}_{\alpha \beta }^{ekl}$$

$$K_{\alpha \gamma }^{ekl}={\int }_{{\Omega }^e}\frac{\partial N^{el}}{\partial X_{\beta }}{\mathcal{A}}_{\alpha \beta \gamma \delta }\frac{\partial N^{ek}}{\partial X_{\delta }}dV.$$

$$\bar{\mathbf{P}}:\Delta \bar{\mathbf{F}}=\frac{1}{V}{\int }_{\partial \Omega }\mathbf{t}\cdot \Delta \mathbf{x}dA.$$

$$\bar{\mathbf{F}}=\frac{1}{V}{\int }_{\partial \Omega }\mathbf{x}\otimes \mathbf{N}dA,\text{and}\bar{\mathbf{P}}=\frac{1}{V}{\int }_{\Omega }\mathbf{P}dV,$$

$$\tilde{\text{u}}=\mathbf{x}-\bar{\mathbf{F}}\mathbf{X},$$

$$\bar{\mathbf{P}}=\frac{1}{V}\sum _e{\int }_{{\Omega }^e}\mathbf{P}dV.$$

$${\bar{\mathbf{A}}}^v=\frac{1}{V}\sum _e{\int }_{{\Omega }^e}{\mathbf{A}}^{v}dv+\frac{1}{V}{\mathbf{L}}^T\frac{\partial \Delta \tilde{\mathbf{u}}}{\partial \Delta {\mathbf{F}}^v},$$

$${\mathcal{L}}_{ϵ\alpha \beta }^{ek}={\int }_{{\Omega }^e}{\mathcal{A}}_{\alpha \beta ϵ\zeta }\frac{\partial N^{ek}}{\partial X_{\zeta }}dV.$$

$$\delta {\tilde{\mathbf{u}}}^T(\mathbf{K}\Delta \tilde{\mathbf{u}}+\mathbf{L}\Delta {\bar{\mathbf{F}}}^v)=0,\text{as}\frac{\partial \Delta \tilde{\mathbf{u}}}{\partial \Delta {\mathbf{F}}^v}=\mathbf{S}=-{\mathbf{K}}^{-1}\mathbf{L},\text{and thus}\mathbf{K}\mathbf{S}=\mathbf{L}.$$

$$\vec{y}={\mathbf{\psi }}^T\vec{u},\text{and}\bar{\mathbf{u}}=\mathbf{\psi }\vec{y},$$

$$\bar{\mathbf{u}}=R(\mathbf{y})=\bar{\mathbf{V}}\mathbf{y}+\tilde{\mathbf{V}}\mathbf{\Xi }\mathbf{q}(\mathbf{y}),$$

$$\underset{\bar{\mathbf{V}},\tilde{\mathbf{V}},\mathbf{\Xi },{\mathbf{y}}_i}{\min }\sum _i\Vert {\mathbf{u}}_i-R({\mathbf{y}}_i){\Vert }^2,$$

$$\underset{\mathbf{W}}{\min }\sum _i\Vert {\vec{u}}_i-\sum _{j\in N_i}{W}_{ij}{\vec{u}}_j{\Vert }^2.$$

$$\underset{\mathbf{Y}}{\min }\sum _i\Vert {\vec{y}}_i-\sum _{j\in N_i}{W}_{ij}{\vec{y}}_j{\Vert }^2.$$

$$\overrightarrow{\stackrel{ˉ}{u}}=\mathbf{\phi }\vec{y}+{\vec{u}}^0,$$

$$\underset{\mathbf{\phi },{\vec{u}}^0}{\min }f(\mathbf{\phi },{\vec{u}}^0)=\underset{\mathbf{\phi },{\vec{u}}^0}{\min }\sum _n\sum _i(\sum _j{\phi }_{ij}{Y}_{Njn}+u_i^0-U_{Nin}{)}^2$$

$$\mathbf{\phi }={\mathbf{U}}_N{\mathbf{W}}_N{\mathbf{Y}}_N^T({\mathbf{Y}}_N{\mathbf{W}}_N{\mathbf{Y}}_N^T{)}^{-1}.\text{with}{\mathbf{W}}_N={\mathbf{I}}_N-\frac{1}{N}{1}_N.$$

$$\overrightarrow{\stackrel{ˉ}{y}}={\bar{\mathbf{\phi }}}^T\vec{u},\text{and}\bar{\mathbf{u}}=\bar{\mathbf{\phi }}\overrightarrow{\stackrel{ˉ}{y}},$$

$$\tilde{\mathbf{\phi }}={\bar{\mathbf{Y}}}_N{\mathbf{W}}_N{\mathbf{Y}}_N^T({\mathbf{Y}}_N{\mathbf{W}}_N{\mathbf{Y}}_N^T{)}^{-1},$$

$$\delta {\bar{\mathbf{u}}}^T\mathbf{g}(\bar{\mathbf{u}};\mathbf{p})=0,\text{subject to}\bar{\mathbf{u}}\in {\mathcal{M}}_{\bar{\mathbf{u}}},\delta \bar{\mathbf{u}}\in {\mathcal{T}}_{\bar{\mathbf{u}}}{\mathcal{M}}_{\bar{\mathbf{u}}},$$

$$\delta {\mathbf{y}}^T{\mathbf{\phi }}^T\mathbf{g}(\bar{\mathbf{u}};\mathbf{p})=0,\text{subject to}\bar{\mathbf{u}}\in {\mathcal{M}}_{\bar{\mathbf{u}}},$$

$${\mathbf{g}}_{\text{red}}(\bar{\mathbf{u}};\mathbf{p})={\mathbf{\phi }}^T\mathbf{g}(\bar{\mathbf{u}};\mathbf{p})=0,\text{subject to}\bar{\mathbf{u}}\in {\mathcal{M}}_{\bar{\mathbf{u}}}$$

Full text

1 Introduction

Engineering tasks ranging from uncertainty quantification [2] over design optimisation [3, 4, 5], the tackling of inverse problems [6, 7, 8, 9, 10] and simple parameter studies to multiscale modeling [11, 12, 13, 14] necessitate repeated evaluations of parameterised simulation models. When an engineering task requires hundreds or even thousands of queries to such a simulation model, the use of traditional high-fidelity simulation techniques, e.g. the Finite Element Method (FEM), may result in prohibitive computational costs. The parametric setting, however, permits the use of parametric model order reduction (MOR) techniques: a reduced order model (ROM) can be constructed based on empirical snapshot data and a priori knowledge about quantities of interest in an offline training phase. The ROM can then be deployed at reduced computational cost in the online application phase.

In this article, we primarily deal with the acceleration of computations on representative volume elements (RVEs) in the context of quasi-static, hyperelastic multiscale modeling, though the methods deployed to this end are intended to be sufficiently general to be applied in different parametric settings as well. In a previous work [1], we outlined a range of approaches to accelerating solid-mechanical multiscale simulations: these include simplified physical models [15, 16], methods which simplify system equations by exploiting the specific physics [17, 18, 19, 20, 21, 22], projection-based approaches which simplify equations in a general, data-based manner [23, 24, 25, 26, 27, 28, 29, 13, 30], and data-based surrogate models [14, 31, 32, 27]. Which approach is most appropriate strongly depends on the problem at hand. The choice of method is generally informed by trade-offs between that method’s offline training cost and data hungriness, online accuracy, online computational cost, robustness and generality, as well as the ease of use and implementation.

As discussed in [1], projection-based MOR techniques strike an attractive balance between these performance characteristics: they require very little data and enable considerable reductions in computational cost while maintaining high degrees of accuracy, robustness and generality. This is achieved at the cost of intrusiveness and comparatively high implementational complexity111Note that, in this work, we take “projection-based MOR” to only denote intrusive methods which to some extent retain the integration of the local, underlying balance laws and thus remain closely tied to the underlying physics. For non-intrusive techniques which strike a different balance between the aforementioned performance criteria, see e.g. [33, 34, 35, 36].. Important structural features of the underlying full-order simulation scheme are retained, while the algorithmic complexity of expensive operations, such as large linear system solves and integration of balance laws over finely resolved calculation domains, is significantly reduced.

Conceptually, projection-based MOR capitalises on the observation that solutions $\bm{u}(\bm{p})\in\mathbb{R}^{D}$ to a well-posed, quasi-static, parametric problem $\bm{g}(\bm{u};\bm{p})=\bm{0}$ which is discretised in terms of $D$ unknowns and formulated in terms of parameters $\bm{p}\in\mathbb{R}^{\delta}$ lie on a $\delta\ll D$-dimensional solution manifold $\mathcal{M}_{\bm{u}}$. Projection-based ROMs attempt to look for solutions not in the high-dimensional solution space $\mathbb{R}^{D}$, but a $d\ll D$-dimensional approximation space (or approximation manifold) $\mathcal{M}_{\bm{\bar{u}}}$, which in the ideal case is a tight superset of the low-dimensional, but possibly strongly nonlinear solution manifold. The Galerkin projection of the unknown variables and system equations onto this approximation space can significantly reduce the computational cost of solving the linear equation systems arising from, e.g. a Newton-Raphson solution scheme.

Intrusive projection-based MOR techniques additionally rely on a hyperreduction methodology to reduce the computational cost of integrating the underlying local balance laws to assemble the aforementioned linear equation systems [37]. These methods generally restrict the evaluation of the underlying local physical equations to a restricted set of mesh entities such as residual degrees of freedom [38, 39, 40, 41, 42, 43, 44, 37, 45, 46, 47, 48], elements [49], Gauss points [25] or, more abstractly, points in strain space [50, 30]. The contribution of mesh entities which are discarded in favour of computational speed can then be approximated explicitly or implicitly. This low-order approximation is possible because the reduced residual is a function of the location of the sought solution on the solution manifold $\bm{u}\in\mathcal{M}_{\bm{u}}$, as well as of the location of the current iterate in the approximation space $\bm{\bar{u}}\in\mathcal{M}_{\bm{\bar{u}}}$ – and thus also lies on a low-dimensional manifold.

Hyperreduced, projection-based MOR algorithms can achieve truly significant reductions in runtime and memory costs when care is taken that their computational complexity is no longer dependent on the size of the underlying full-order problem. In the following, we outline established Galerkin- and hyperreduction techniques which can be utilised to this end. Additionally, further discussion of algorithmic complexity is integrated into the remainder of this work.

In the context of intrusive, projection-based MOR for quasi-static solid mechanics, the Proper Orthogonal Decomposition (POD) is a popular approach to performing the Galerkin projection. It defines the approximation space as a linear subspace of the original solution space [51, 52, 53, 54]. While the POD is attractive on account of its simplicity and robustness, relatively high model sizes $d>\delta$ might be required to accurately represent a strongly nonlinear solution manifold. The local basis method (LPOD) thus aims for a tighter representation of the solution manifold by localising the POD to regions of solution- or parameter space [55, 56, 57, 58, 59, 4, 5]. Alternatively, polynomial manifold (PM) approaches attempt to achieve a smoother, tighter approximation via a polynomial Ansatz [60, 61, 33, 62]. Artificial Neural network (ANN) approximation spaces (using Feedforward Neural Networks or Autoencoders) have also been used with a view to improved flexibility and generality, but are relatively data-hungry [63, 64, 65, 66].

In a recent proof-of-concept [1], we instead proposed a graph-based manifold learning (GML) approach for the nonlinear Galerkin-reduction of quasi-static solid-mechanical problems. GML techniques such as Locally Linear Embedding (LLE) [67] strike an attractive balance between the flexibility of an ANN Ansatz and the data-economy of the POD. As a result, they have successfully been applied to the surrogate modeling of RVEs [14, 31, 32] and fluid problems [68] as well as to non-intrusive model reduction in elastodynamics [69]. As opposed to the approaches outlined above, the LLE and related methods do not utilise an explicit Ansatz for the approximation space, instead only providing a low-dimensional embedding of the snapshot data. Approximation spaces have to be constructed a posteriori. [69] used max-Ent interpolants to this end, [70] a LPOD-like linearisation, and [1] a forward-Euler-like scheme.

While research on nonlinear approximation spaces is quite well-developed for the dynamic case, particularly in fluid mechanics [71, 55, 56, 59, 72, 73, 74, 61, 62, 63, 75, 76, 77, 78, 79, 66, 70], their application to (quasi-static) solid mechanics in general and computational homogenisation in particular are in a stage of relative infancy. Dynamic solid problems have been tackled with LPOD approaches in [80, 59, 81, 72], with the PM in [73, 61], and with a (non-intrusive) manifold learning scheme in [69]. The LPOD has also been applied to quasi-static solid mechanics in [82, 5] and a time-dependent computational homogenisation problem in [83]. The PM has been applied to a quasi-static problem by [84], but its performance in high-dimensional parameter spaces and in combination with hyperreduction techniques remain to be explored. Similarly, we applied a proof-of-concept GML approach to an RVE problem in [1], but the hyperreduction, as well as the homogenisation of this approach remain to be addressed.

It is important to note that the context of employment results in some specific challenges for nonlinear MOR techniques: in quasi-static applications, the search for a static solution, rather than the evolution of a dynamic one, is constrained to the approximation space [1]. Data can be scarce, and the accurate approximation of the reduced system Jacobian, as well as the convergence of the reduced scheme, are critical. In the context of multiscale modeling, the parameter space can relatively large (at a minimum, $\delta=6$), and an efficient and accurate post-processing of homogenised stresses is paramount.

By comparison, there has been extensive research on hyperreduction techniques in general and for the quasi-static solid case in particular. Approaches which assemble a reduced set of entries of the residual vector and approximate the reduced residual in an approximated Galerkin scheme include the Gappy POD [38, 39], Best Points Interpolation [40], Missing Point Estimation [41], Discrete Empirical Interpolation Method (DEIM) [42], unassembled DEIM [43], and S-OPT [44]. While these approaches are attractive on account of their simplicity, a lack of preserved structure in the resulting reduced equation systems has been noted to result in suboptimal robustness [26, 85]. Methods for structure-preservation have been proposed for specific problem classes (usually in dynamics) [86, 87, 88], but this may come at the cost of reduced generality and increased implementational effort. In contrast, collocation-like Petrov-Galerkin methods solve the underlying equation systems at the selected entries of the residual (either in a weighted or least-squares manner) more directly. Such methods include the eponymous hyperreduction (HR) technique from [37], Gauss-Newton with Approximate Tensors (GNAT) [45, 46], Least-Squares Petrov Galerkin (LSPG) [47], and Reduced Over-Collocation [48]. Some of these methods strike a promising balance between simplicity, robustness [26, 85] and generality. In the context of the FEM, Energy Preserving Weighting and Sampling (ECSW), which considers the contribution of a restricted set of elements to the reduced residual, provides a natural, robust, structure-preserving alternative [49]. The Empirical Cubature Method generalises this approach to integration points and augments ECSW with additional physical constraints [25]. Finally, promising strain-based hyperreduction techniques such as the Empirically Corrected Cluster Cubature were recently proposed in [50, 30, 30]. These approaches select integration points at which to evaluate the material law as statistical representatives of the set of integration points in strain space and allow for natural and highly efficient hyperreduction especially in the context of computational homogenisation.

As mentioned above, such hyperreduction techniques are usually applied in tandem with the POD in the context of computational homogenisation. For example, approaches assembling a reduced set of residual degrees of freedom are applied in [89, 23, 90, 91, 26, 92, 85, 93] while a reduced set of elements are assembled by [24, 25, 92, 94, 27, 28, 29, 13, 95, 96], and [97, 50, 30, 30] apply strain-based approaches. Additionally, the LPOD has been coupled with an element-based approach by [83]. In quasi-static solid mechanics more generally, the POD has been applied with residual-based hyperreduction approaches in [98, 82, 99, 85, 100, 101, 102, 103, 104] and element approaches in [105, 106, 107, 108, 109]. Furthermore, the LPOD has been deployed on quasi-static problems with DEIM-like approaches by [82, 80, 5]. On dynamic problems, the LPOD has also been variously applied with the DEIM [59], gappy POD [59], GNAT [55, 56], and ECSW [81, 72]. The PM, meanwhile, has been coupled with DEIM-like approaches [73, 74], and the ECSW [61, 62]. Meanwhile, ANN approaches have been deployed with S-OPT [76], HR [75], GNAT [63, 76], LSPG [77, 78, 79], and ECSW [66, 76]. There have also been some efforts to exploit nonlinear correlations in the hyperreduction step itself, see e.g. [110, 111, 112].

In this work, we build on the proof-of-concept GML approach to nonlinear MOR proposed in [1]. In order to reduce the computational cost of assembling reduced equation systems, we augment the scheme with two popular hyperreduction techniques – namely, the DEIM and the LSPG. We emphasise aspects of the implementation which are necessary to achieve significant reductions in runtime and memory costs, and discuss algorithmic complexity where appropriate. Additionally, we extend the nonlinear MOR scheme to allow for the efficient computation of homogenised stresses and stiffnesses for multiscale applications. Finally, we compare the performance of the resulting hyperreduced nonlinear MOR scheme on a simple example RVE problem with several alternative techniques: the POD, LPOD, and PM. On the example RVE problem, the hyperreduced GML approach is shown to yield speedups of two orders of magnitude with respect to full FEM simulations while retaining high levels of accuracy. Furthermore, an LLE approach is shown to be able to Pareto-dominate alternative methods in terms of the tradeoff between computational efficiency, accuracy, and robustness.

The remainder of this work is structured as follows: Sections 2 and 3 outline the fundamentals of the underlying quasi-static solid mechanical- and computational homogenisation problems, respectively. Section 4 then discusses the nonlinear projection-based MOR techniques to be investigated, including their application to the Galerkin-reduction of a Newton-Raphson solver scheme and to the homogenisation of stress and stiffness. Section 5 covers the hyperreduction techniques. Finally, Section 6 features numerical comparisons of the MOR and hyperreduction techniques on an example RVE problem.

2 Quasi-static solid-mechanical problems

On a solid domain $\Omega$ consisting of points $\bm{X}$ in the reference configuration, the weak form of the quasi-static balance of linear momentum can be written as

[TABLE]

with $\bm{P}$ denoting the first Piola-Kirchoff stress, $\bm{F}$ the work-conjugate deformation gradient, $\bm{b}$ a volumetric load, u the displacement, $\bm{t}$ a nominal traction on boundary $\partial\Omega$, and $dA$ and $dV$ infinitessimal surface and volume elements, respectively. $\delta\bm{\text{\bf{u}}}$ and $\delta\bm{F}$ refer to kinematically admissible variations in u and $\bm{F}=\frac{\partial(\bm{X+\text{\bf{u}}})}{\partial\bm{X}}$, respectively, under which the variation of the virtual work must vanish [113, p.82].

A material-dependent constitutive relation $\bm{P}=\bm{P}(\bm{F})$ links strain and stress; in the hyperelastic case with which the current work is concerned, the first Piola-Kirchhoff stress is assumed to be derivable from a stored energy function $W$ as the first derivative $\bm{P}=\frac{\partial W}{\partial\bm{F}}$ [114, p.207]. For the numerical examples in Section 6, we consider a simple compressible neo-Hooke potential

[TABLE]

with $I_{c}=F_{\beta\alpha}F_{\beta\alpha}$, $J=\text{det}(\bm{F})$, and $\mu$ and $\kappa$ being the shear and bulk moduli. The fourth-order nominal stiffness tensor $\bm{\mathcal{A}}=\frac{\partial\bm{P}}{\partial\bm{F}}$ can then also be derived from this potential as the second derivative $\bm{\mathcal{A}}=\frac{\partial^{2}W}{\partial\bm{F}^{2}}$.

The discretisation of Eq. (1) via a standard isoparametric Galerkin Ansatz with ${\text{{u}}}^{e}(\bm{X})=\sum_{k}N^{ek}(\bm{X})\bm{u}^{ek}$ over a set of elements $e\in E$ yields

[TABLE]

where $N^{ek}$ denotes the shape function corresponding to node $k$ in element $e$, with $\delta\bm{u}^{ek}$ being the corresponding nodal displacement and $\bm{g}^{ek}$ the associated nodal residual [113, p.101-125]. When all nodal variations $\delta{{u}}_{\alpha}^{ek}$ are collected in a vector $\delta\bm{u}$ without duplicating nodal degrees of freedom shared between elements, a global residual $\bm{g}$ can be defined. $\bm{g}$ can be assembled with a complexity of roughly $\mathcal{O}(|E|\mathcal{C}_{e}^{g})$ where $|E|$ is the total number of elements and $\mathcal{C}_{e}^{g}$ the cost of evaluating Eq. (3), as

[TABLE]

where $K$ denotes the global node identifier associated with node $k$ in element $e$. Thus, the weak form becomes

[TABLE]

and since this must vanish for any kinematically admissible $\delta\bm{u}$,

[TABLE]

Here, $\bm{p}\in\mathbb{R}^{\delta}$ denotes a vector of parameters which specifies an instance of the parametric problem class and implicitly determines a concrete solution $\bm{u}(\bm{p})$. $\bm{p}$ might for example specify the material constants $\kappa(\bm{X})$ and $\mu(\bm{X})$ or boundary conditions, e.g. via the traction $\bm{t}$ on boundary $\partial\Omega$.

For sufficiently well-behaved problems, roots of Eq. (5) can be sought using a classical Newton-Raphson scheme, in which case Eq. (5) is linearised around the current iterate $\bm{u}_{\text{cur}}$ and an increment $\Delta\bm{u}$ computed via [113, p.148]

[TABLE]

The Jacobian or stiffness matrix $\bm{K}$ appearing in the above is defined as

[TABLE]

which can be assembled analogously to $\bm{g}$ at $\mathcal{O}(|E|\mathcal{C}_{e}^{K})$, with the cost of computing the element stiffness $\mathcal{C}_{e}^{K}$, as

[TABLE]

in terms of an element stiffness matrix $\bm{K}^{ekl}$ (in index notation, using Einstein summation convention)

[TABLE]

3 Computational homogenisation with RVEs

Computational homogenisation techniques such as the FE2 method model microstructural mechanical processes and their effect on macroscopic behaviour via a characteristic RVE. The evaluation of a material law on the macroscale $\bm{\bar{P}}(\bm{\bar{F}})$ is replaced by the solution of an RVE problem subject to suitable boundary conditions and the subsequent computation of the average stress and the associated stiffness. In the following, we review aspects of computational homogenisation which are important for our discussion of Galerkin- and hyperreduction later; please consult the reference literature for more detail [115, 116]. Below, macroscale quantities will be denoted by an overbar $\bar{\cdot}$, while microscale quantities remain undecorated.

The Hill-Mandel condition constitutes a physically sensible scale-coupling relation: the work done on a point $\bm{\bar{X}}$ on the macroscale via a strain perturbation $\Delta\bm{\bar{F}}$ must equal the work done on an associated microscale RVE via an associated perturbation $\Delta\bm{x}$ of the current position $\bm{x}(\bm{X})$ on the microscale, i.e. [115]

[TABLE]

Therein, the macroscopic deformation gradient $\bm{\bar{F}}$ and first Piola-Kirchhoff stress $\bm{\bar{P}}$ can be defined as

[TABLE]

where $\bm{x}$ is the deformation on the microscale, assuming no tractions act on pores within the RVE [117, 115].

In this work, we employ classical periodic boundary conditions, which satisfy Eq. (9) in a physically sensible way [115]. Thus, the displacement fluctuation

[TABLE]

is assumed to be periodic with the RVE acting as a unit cell. We do not discuss implementational details here, since boundary handling is not the focus of this work. The vector $\bm{\tilde{u}}\in\mathbb{R}^{D}$, which results from the discretisation of $\tilde{\text{{u}}}$ and the application of boundary conditions, then becomes the primary unknown of the RVE problem.

To the accuracy of the FE discretisation, the macroscopic first Piola-Kirchhoff stress $\bm{\bar{P}}$ can be postprocessed from the stress field on the RVE by volume averaging with around $\mathcal{O}(|E|\mathcal{C}_{e}^{P})$, with $\mathcal{C}_{e}^{P}$ being the cost of computing the element-level stress contribution (see Eq. (10))

[TABLE]

The macroscopic nominal stiffness $\bm{\mathcal{\bar{A}}}$ is defined as the rate of the change of $\bm{\bar{P}}$ with $\bm{\bar{F}}$, i.e. $\bm{\mathcal{\bar{A}}}=\frac{\partial\bm{\bar{P}}}{\partial\bm{\bar{F}}}$ [115]. When written in Voigt notation and when the FE discretisation is substituted, this stiffness becomes

[TABLE]

with the sensitivity coefficient $\bm{L}\in\mathbb{R}^{D\times 9}$, which can be assembled analogously to $\bm{g}(\bm{\tilde{u}};\bm{F})$ and $\bm{K}(\bm{\tilde{u}};\bm{F})$ from element contributions, in this case

[TABLE]

For further details, please consult the reference literature, e.g. [118, 115].

The first term in Eq. (13) can be assembled with around $\mathcal{O}(|E|\mathcal{C}_{e}^{A})$, and and the degrees of freedom of $\bm{L}$ scale with about $\mathcal{O}(|E|\mathcal{C}_{e}^{L})$, with $\mathcal{C}_{e}^{A}$ and $\mathcal{C}_{e}^{L}$ denoting the cost of the respective element-level operations. The sensitivity $\frac{\partial\Delta\bm{\tilde{u}}}{\partial\Delta\bm{F}^{v}}$ appearing in the structural softening term can be computed as solutions to the variational problem

[TABLE]

$\bm{S}\in\mathbb{R}^{D\times 9}$ is thus computed as the solution to 9 equation systems, where the matrix $\bm{K}$ is the same for all 9 cases, meaning that an appropriate sparse factorisation must be performed only once [115]. The complexity of this step is of around $\mathcal{O}(D^{2})$. The operations and computational costs involved in solving an RVE problem subject to periodic boundary conditions and homogenising the stress and stiffness are summarised in Alg. 1. For the sake of notational simplicity, no further distinction will be made between $\bm{u}$ and $\bm{\tilde{u}}$ in the following; the vector of the primary unknowns will be denoted as $\bm{u}$ for the sake of generality.

4 Nonlinear Model Order Reduction

The computational cost of the iterative solver procedure outlined in Alg. 1 is dominated by the assembly procedure with $\mathcal{O}(|E|\mathcal{C}_{e})$, where $\mathcal{C}_{e}=\mathcal{C}_{e}^{g}+\mathcal{C}_{e}^{K}$ and the solution of the linear system with around $\mathcal{O}(D^{2})$. The increments also scale with $\mathcal{O}(D)$, albeit with a small constant coefficient. As a practical contributor to computational cost, the number of iterations until convergence is also consequential, since runtime scales roughly linearly with it. The cost of post-convergence homogenisation, meanwhile, is dominated by the assembly of $\bm{L}$ with $\mathcal{O}(|E|\mathcal{C}_{e}^{L})$, the solution of 9 linear problems with around $\mathcal{O}(D^{2})$, stress and stiffness integration with $\mathcal{O}(|E|\mathcal{C}_{e}^{P})$ and $\mathcal{O}(|E|\mathcal{C}_{e}^{K})$, and (less significantly) the homogenised stiffness computation at $\mathcal{O}(D)$. The cost of element-level operations $\mathcal{C}_{e}^{g},\mathcal{C}_{e}^{K},\mathcal{C}_{e}^{P}$ and $\mathcal{C}_{e}^{A}$ are similar in magnitude, and, while vanishing compared to system-level operations, should not be underestimated: the evaluation of the constitutive relation, Ansatz functions, tensor operations, and numerical integration accumulate costs and runtime which are incurred for every element $e\in E$. Classical implementational considerations such as assembling $\bm{g}$ and $\bm{K}$ as well as integrating $\bm{P}$ and $\bm{A}^{v}$ in tandem mitigate this slightly, but not qualitatively.

In a parametric multi-query context – e.g., computational homogenisation – mitigating computational costs is crucial not only to accelerate investigations, but to make certain research feasible in the first place. As argued in the Introduction and in [1], projection-based MOR methods applied in combination with hyperreduction techniques strike an attractive balance between computational cost reductions, accuracy, data-economy, low offline cost, and generality in this circumstance. Broadly speaking, these approaches truncate the problem size by reducing the size of the solution space and the domain of integration. If constructed well, this results in algorithms with a computational cost which no longer scales with the size of the original problem. In the remainder of this Section, we discuss the first of the two major steps necessary to this end, while hyperreduction is covered in Section 5.

4.1 Dimension reduction and representation learning problem

The set of all solutions $\bm{u}$ to a parametric, quasi-static solid-mechanical problem defines a solution manifold $\mathcal{M}_{\bm{u}}=\{\vec{{u}}\mid\exists\,\bm{p}:\bm{g}(\vec{{u}};\vec{p})=\vec{0}\}$ for parameters $\vec{p}\in\mathbb{R}^{\delta}$ [57, 58]. If the problem is well-posed, this solution manifold is $\delta$-dimensional, meaning that solutions lie in a significantly lower-dimensional nonlinear subspace of the $D$-dimensional solution space. In the case of a hyperelastic RVE problem, for example, $\bm{p}$ contains the degrees of freedom of the macroscopic deformation gradient $\bm{\bar{F}}\in\mathbb{R}^{3\times 3}$ which are not responsible for rigid body rotations, such that $\delta=6$ [1]. Projection-based MOR techniques attempt to exploit this observation by searching for solutions in a $d$-dimensional approximation space $\mathcal{M}_{\bm{\bar{u}}}$, where $d\ll D$. In the ideal case, $\mathcal{M}_{\bm{\bar{u}}}\supset\mathcal{M}_{\bm{u}}$, such that all possible solutions lie in $\mathcal{M}_{\bm{\bar{u}}}$.

The solution manifold $\mathcal{M}_{\bm{u}}$ is of course not known a priori, such that projection-based MOR techniques have to have recourse to $s$ discrete snapshot solutions $\bm{U}=[\vec{u}_{1},..,\vec{u}_{s}]\in\mathbb{R}^{D\times s},\vec{u}_{i}\in\mathcal{M}_{\bm{u}}$ gathered from the underlying high-fidelity model in an offline training phase. The dimension reduction problem at the heart of nonlinear projection-based MOR could then be phrased as follows: given the snapshot data, find an embedding map $M:\mathbb{R}^{D}\rightarrow\mathbb{R}^{d}$ such that the embedding $\bm{y}=M(\bm{{u}})$ retains as much as possible of the structure of the solution manifold. Nonlinear projection-based MOR techniques can then seek solutions in the low-dimensional reduced space $\bm{y}\in\mathbb{R}^{d}$. To this end, a reconstruction map $R:\mathbb{R}^{d}\rightarrow\mathbb{R}^{D}$ defining approximate solutions $\bm{\bar{u}}=R(\bm{y})$ in the original solution space is also required, which turns the dimension reduction problem into a representation learning problem: find embedding and reconstruction maps $M:\mathbb{R}^{D}\rightarrow\mathbb{R}^{d}$ and $R:\mathbb{R}^{d}\rightarrow\mathbb{R}^{D}$, such that, for $d\ll D$, the reconstruction error $\sum_{i}\|\bm{\bar{u}}_{i}-\bm{u}_{i}\|=\sum_{i}\|R(M(\bm{u}_{i}))-\bm{u}_{i}\|$ is minimised. The set of reconstructions $\mathcal{M}_{\bm{\bar{u}}}=\{\bm{\bar{u}}\mid\bm{\bar{u}}=R(\bm{y}),\,\bm{y}\in\mathbb{R}^{d}\}$ then defines the aforementioned approximation space which is parameterised via the reduced variables $\bm{y}\in\mathbb{R}^{d}$.

4.2 Approaches to defining approximation spaces

In the following, we briefly cover several established approaches to defining approximation spaces, before recalling our manifold learning approach. For details, the interested reader is referred to the respective source material and to our previous work [1].

Proper Orthogonal Decomposition (POD)

As motivated in the Introduction, the POD [119, 54, 120] is a popular approach to defining low-dimensional approximation spaces in the context of quasi-static solid mechanics. It defines linear embedding and reconstruction mappings

[TABLE]

where the mode matrix $\bm{\psi}$ can e.g. be obtained via a singular value decomposition $\bm{U}=\bm{L}\bm{\Sigma}\bm{R}$, as the leading $d$ left singular vectors $\bm{\psi}=[\bm{L}_{1},..,\bm{L}_{d}]\in\mathbb{R}^{D\times d}$ of the snapshot matrix. Fig. 2 illustrates this approach on a low-dimensional data set.

Local Basis Method (LPOD)

As noted in the Introduction, the POD might require comparatively high-dimensional reduced- $\bm{y}\in\mathbb{R}^{d}$ and approximation spaces $\mathcal{M}_{\bm{\bar{u}}}$ to accurately represent a highly nonlinear, $\delta$-dimensional solution manifold $\mathcal{M}_{\bm{u}}$. The LPOD [55, 56] aims to reduce the gap between $d$ and $\delta$ by localising the POD: snapshots $\bm{u}_{i}$ are clustered and local POD bases defined for each cluster. This conceptual approach is visualised in Fig. 2. For the variant of the LPOD algorithm used here, as well as for parameter studies which were used to inform our choices for the (copious) model parameters, see [121]. Despite its successes (for examples in quasi-static solid mechanics, see e.g. [82, 80, 5]), the LPOD is subject to some limitations, especially in the data-poor regime. As we noted in previous work [1, 121], in addition to the challenges of handling local basis transitions [122, 56, 57, 58, 121], choosing parameters properly and assuring cluster quality is nontrivial and locally linear approximations might not represent a solution manifold as closely as may be possible.

Polynomial manifold approach (PM)

The PM [60, 61, 33, 62] aims to avoid some of the shortcomings of the LPOD with a continuously nonlinear approximation space: the reconstruction map is defined via a polynomial Ansatz

[TABLE]

where $\bm{\bar{V}}\in\mathbb{R}^{D\times d}$ and $\bm{\tilde{V}}\in\mathbb{R}^{D\times\tilde{d}}$ are orthonormal basis matrices. $\bm{q}(\bm{y})\in\mathbb{R}^{p}$ contains polynomial terms of $\bm{y}$, e.g. monomials or full Kronecker products. Here, we use (vectorised) quadratic Kronecker products, i.e. $\bm{p}(\bm{y})=\text{vec}(\bm{y}\otimes\bm{y)}$. $\bm{\Xi}\in\mathbb{R}^{\tilde{d}\times p}$ is a low-dimensional coefficient matrix for this nonlinear term, with $p$ being the number of polynomial terms. The coefficients $\bm{\bar{V}},\bm{\tilde{V}}$, and $\bm{\Xi}$, as well as the embedding $\bm{y}_{i}$ of the snapshots $\bm{u}_{i}$ can be obtained by fitting the Ansatz to the snapshot data

[TABLE]

e.g. via the alternating minimisation approach from [33].

4.3 Manifold learning approach

While the continuously nonlinear approximation space obtained by the PM is very attractive in the context of nonlinear MOR, the use of a specific polynomial Ansatz curtails the flexibility of this approach and to some extent limits the PM to solution manifolds which can be represented well globally by low-order polynomials. Of course, an ANN Ansatz could be substituted instead, but these are comparatively data-hungry and hence not suitable for the data-poor context in which we are interested here [60, 61, 33, 62].

In [1], we instead proposed a proof-of-concept for a graph-based manifold learning approach to the nonlinear MOR of quasi-static solid mechanical systems. Such GML schemes work by constructing a graph adjacency matrix $\bm{G}\in\mathbb{B}^{s\times s}$ for the snapshot data $\bm{U}\in\mathbb{R}^{D\times s}$, e.g. via a k-nearest neighbours approach. Then, an embedding $\bm{Y}\in\mathbb{R}^{d\times s}$ which conserves essential structural characteristics of $\bm{U}$ and $\bm{G}$ is found. This embedding implicitly defines the reduced space $\mathbb{R}^{d}$ in which a projection-based MOR scheme can subsequently search for solutions. For a sufficiently well-behaved underlying (solution) manifold, it is often possible to successfully obtain $\bm{Y}$ based on only a few snapshots $s$, without recourse to a specific Ansatz [67, 123, 124, 125]. GML schemes thus combine the flexibility of an ANN approach with the data-economy of the POD.

Suitable GML schemes for our setting include Locally Linear Embedding (LLE) [67], ISOMAP [123], Laplacian Eigenmaps (LEM) [124] and Local Tangent Space Alignment (LTSA) [125]. In this work, we exclusively employ LLE. After defining a local neighbourhood $N_{i}=\{j\mid G_{ij}=1\}$ of snapshot $i$, LLE constructs a reconstruction weight matrix $\bm{W}$ which optimally reconstructs each $\bm{u}_{i}$ from its neighbours $\bm{u}_{j}$, i.e.

[TABLE]

LLE then computes the embedding $\bm{Y}$ which optimally respects these reconstruction weights in the reduced space, i.e.

[TABLE]

Note that solutions to these optimisation problems can be computed at low cost [67]. For further detail, the interested reader is referred to [67, 1]. An illustration of dimensionality reduction via the LLE can be found in Fig. 4.

LLE achieves flexibility in the construction of the embedding at the cost of not obtaining explicit embedding $M:\mathbb{M}^{D}\rightarrow\mathbb{R}^{d}$ and reconstruction $R:\mathbb{R}^{d}\rightarrow\mathbb{R}^{D}$ maps. In [1], we obtained a reconstruction via a simple data-based local linearisation scheme deployed in the online phase. By linearising around the current value of the reduced variable $\bm{y}_{\text{cur}}$, we obtain a local model

[TABLE]

with tangent $\bm{\varphi}\in\mathbb{R}^{D\times d}$ based on the reduced $\vec{y}_{i}\in\mathbb{R}^{d},i\in N_{y}$ and original coordinates $\vec{u}_{i}\in\mathbb{R}^{D},i\in N_{y}$ of the $N$ snapshots $N_{y}$ nearest to ${\vec{y}}_{\text{cur}}$. The neighbour search among $\bm{y}_{i}\in\mathbb{R}^{d},i\in 1,..s$ can be computed at a low cost of $\mathcal{O}(sN)$. Minimising the least squares error

[TABLE]

yields

[TABLE]

An illustration of the local linearisation scheme can be found in Fig. 5.

In [1], we performed this linearisation in each Newton iteration in an Euler forward-like scheme. Here, we instead note that at the beginning of each load step within an RVE computation, the macroscopic deformation gradient $\bm{\bar{F}}$ at the targeted solution is known (see Alg. 1). More generally, in a quasi-static solid mechanical problem, the parameters $\bm{p}$ determining the sought solution $\bm{u}$ are available at the beginning of each load step. Thus, we can perform the neighbour search among the snapshot parameter values $\bm{\bar{F}}_{i},i\in 1,..,s$ rather than among the reduced snapshots $\bm{y}_{i}\in\mathbb{R}^{d},i\in 1,..s$, and obtain an approximate Euler backward linearisation scheme. Unsurprisingly, this improves convergence and reduces runtime, since linearisations within the Newton scheme are avoided. We have also investigated more sophisticated strategies, such as an approximate Crank-Nicolson-like linearisation, but did not find these to be necessary in this work.

Finally, we deploy the embedding and reconstruction operations in a two-stage manner. Firstly, a linear compression to an intermediate, $\bar{d}$-dimensional space can be performed via the POD

[TABLE]

followed by dimensionality reduction via the LLE into the reduced space $\mathbb{R}^{d}$. If $\bar{d}=s$, the first stage becomes a lossless compression into the span of the snapshots. Then, a local online linearisation can be performed to the intermediate space in the online phase at $\mathcal{O}(\bar{d}Nd)$ via

[TABLE]

where $\bm{\bar{Y}}$ denotes snapshots in the intermediate space, after compression via the POD. A two-stage local linearisation with $\bm{\varphi}=\bm{\bar{\varphi}}\bm{\tilde{\varphi}}$ can thus be computed at computational cost which no longer scales with $\mathcal{O}(D)$, with no (or, if $d<\bar{d}<s$, little) additional error.

4.4 Galerkin-reduced solution procedure

With the solution to the representation learning problem at hand and the nonlinear approximation space $\mathcal{M}_{\bm{\bar{u}}}$ defined, we can turn our attention to the search for solutions within this approximation space. We aim to search for solutions to Eq. (4) in $\mathcal{M}_{\bm{\bar{u}}}$, i.e.

[TABLE]

with the variation $\delta\bm{\bar{u}}$ in the tangent space of the approximation manifold. With the approximation space parameterised in terms of $\bm{y}$ via the reconstruction $\bm{\bar{u}}=R(\bm{y})$, we have \delta\bm{\bar{u}}=\frac{\partial R(\vec{y})}{\partial\vec{y}}\Bigr{|}_{\vec{y}}\delta\bm{y}. Approximating the tangent to the approximation space via a local linearisation, i.e. \bm{\varphi}\approx\frac{\partial R(\vec{y})}{\partial\vec{y}}\Bigr{|}_{\vec{y}}, we obtain

[TABLE]

which defines a reduced equation system

[TABLE]

with the reduced residual $\bm{g}_{\text{red}}$. As in the full-order system, we can use a Newton-Raphson scheme to seek roots of Eq. (22), i.e.

[TABLE]

The reduced Jacobian $\bm{K}_{\text{red}}$ can be obtained by linearisation as

[TABLE]

As discussed in [1], the cost of solving the reduced equation system in Eq. (23) now scales with $\mathcal{O}(d^{3})$ instead of around $\mathcal{O}(D^{2})$, which slashes computational cost significantly. The assembly of $\bm{g}$ and $\bm{K}$, however, still scales with $\mathcal{O}(|E|\mathcal{C}_{e})$. Less significantly, the computation of $\bm{g}_{\text{red}}$ via Eq. (22) and $\bm{K}_{\text{red}}$ via Eq. (24) scale with $\mathcal{O}(D)$ if performed naïvely. The reduction of these costs using hyperreduction techniques is addressed in Sect. 5. Before this can be done, however, the Galerkin projection of the homogenisation procedure must be addressed.

4.5 Galerkin-reduced homogenisation

The macroscopic algorithmically consistent stiffness can be computed via Eq. (13). With the (nonlinear) Galerkin approximation, this becomes

[TABLE]

with $\bm{L}_{\text{red}}=\bm{\varphi}^{T}\bm{L}$. The sensitivity $\bm{S}_{\text{red}}=\frac{\partial\Delta\bm{y}}{\partial\Delta\bm{F}^{v}}$, meanwhile, can be computed via Eq. (15)

[TABLE]

such that

[TABLE]

Again, the solution of the linear equation systems now scales with $\mathcal{O}(d^{3})$ rather than around $\mathcal{O}(D^{2})$. However, the assembly of $\bm{L}$ still scales with $\mathcal{O}(|E|\mathcal{C}_{e}^{L})$ and the projection with $\mathcal{O}(D)$ if done naïvely. The computation of $\bm{\bar{P}}$ and $\bm{\bar{A}}^{v}$ via Eq. (12) and Eq. (13) also still scale with $\mathcal{O}(|E|\mathcal{C}_{e}^{P})$ and $\mathcal{O}(|E|\mathcal{C}_{e}^{A})$, respectively. Thus, an additional hyperreduction step also is required for these homogenisation operations.

5 Residual vector-based hyperreduction

In simple terms, the Galerkin projection reduces the computational cost of solving the equation systems to which the Newton-Raphson solver and the homogenisation procedure give rise, but not the cost of evaluating constitutive relations and integrating local balance laws to assemble these equation systems. As motivated in the Introduction, we further need to reduce the domain of integration via a second, hyperreduction step to achieve independence from the problem sizes $D$ and $|E|$ in the computational complexity and to achieve truly significant runtime reductions [37].

In the context of the FEM, classical hyperreduction techniques such as the Gappy POD [38, 39], the Best Points Interpolation Method [40], Missing Point Estimation [41], Reduced Integration Domain [37], Discrete Empirical Interpolation Method (DEIM) [42, 43], S-OPT [44] Least Squares Petrov Galerkin [45, 47] or Reduced Over-Collocation [48], Gauss-Newton with Approximate Tensors [46], Energy Conserving Weighting and Sampling (ECSW) [49], Empirical Cubature [25], Statistically Compatible Hyperreduction [50], and Empirically Corrected Cluster Cubature [30] do this by reducing the effective number of elements $|E_{m}|<|E|$ or Gauss points which are considered for numerical integration. A low-order approximation can be successful because the reduced residual is a function of the location of the sought solution on the solution manifold $\bm{u}\in\mathcal{M}_{\bm{u}}$, as well as of the location of the current iterate in the approximation space $\bm{\bar{u}}\in\mathcal{M}_{\bm{\bar{u}}}$ – and thus also lies on a low-dimensional manifold. Here, we consider only hyperreduction methods which accomplish this by assembling only few $m\ll D$ rows of the original nonlinear equation systems [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48]. In keeping with existing literature in the context of the popular DEIM hyperreduction technique, we denote this set of degrees of freedom $P_{m}=\{m_{1},...m_{m}\}$ of $\bm{K}$, $\bm{g}$, and $\bm{L}$ which are assembled exactly as magic points.

5.1 Magic points selection

The DEIM [42] (for application to quasi-static solid mechanics, see also [98, 99]) uses a POD of residual snapshots to characterise the residual manifold via modes. Then, some characteristic entries of these modes (i.e., magic points) are selected via a greedy procedure. In this Subsection, we outline the magic point selection, which can also be employed in different hyperreduction schemes; the definition of the resulting hyperreduced equation systems are discussed in the following Subsections.

To this end, $s_{g}$ residual snapshots $\bm{G}=[\bm{g}_{1},..,\bm{g}_{s_{g}}]\in\mathbb{R}^{D\times s_{g}}$ are gathered from a Galerkin ROM. In contrast to the time-dependent case, in which physical dynamics are to be approximated, we later aim to characterise the “dynamics” of a Newton-Raphson solver on the approximation space via the residual approximation. To faithfully represent the residual manifold, it is therefore paramount to collect residual snapshots from the intermediate, non-converged steps of the Newton-Raphson scheme. A POD of these residual snapshots $\bm{G}$ then yields $m$ modes $\bm{\Omega}\in\mathbb{R}^{D\times m}$ and mode coefficients $\bm{y}^{\bm{g}}\in\mathbb{R}^{m\times s_{g}}$ such that

[TABLE]

with $\bm{\Omega}=[\bm{\omega}_{1},\bm{\omega}_{2},...,\bm{\omega}_{m}]$, where $m$ is the desired number of magic points.

Degrees of freedom which are predictive of these modes are selected iteratively. The first magic point is the maximal degree of freedom of the first POD mode

[TABLE]

We also define a sampling matrix $\bm{Z}$, the first column of which is unit vector in $m_{1}$-th coordinate direction $\bm{Z}=[\bm{e}_{m_{1}}]$. For subsequent entries, we find the mode activity coefficients for the previous modes $\bm{\tilde{\Omega}}=[\bm{w}_{1},..,\bm{\omega}_{j-1}]$ which best predict the magic points of the current mode

[TABLE]

and subtract this contribution from the current mode

[TABLE]

to define the component $\tilde{\bm{\omega}}_{j}$ of the current mode $\bm{\omega}_{j}$ which is not predicted well by previous magic points. Subsequent magic points are the maximal degrees of freedom of this, as of yet unpredicted component

[TABLE]

and the selection matrix is augmented with $\bm{Z}=[\bm{Z},\bm{e}_{m_{j}}]$.

The set of elements $E_{m}$ featuring magic points and the associated nodes and degrees of freedom $I_{m}$ define a reduced integration domain, via which the residual at the magic points $\bm{g}_{m}=\bm{Z}^{T}\bm{g}$ can be evaluated. Then, the hyperreduced assembly of $\bm{g}_{m}$ scales roughly with $\mathcal{O}(|E_{m}|\mathcal{C}_{e}^{g})$, rather than with $\mathcal{O}(|E|\mathcal{C}_{e}^{g})$.

Meanwhile, the product of the magic rows of the stiffness matrix with the modes $\bm{\psi}$ can be computed on the element level. Denoting the product of $\bm{K}$ and $\bm{\psi}$ as $\widehat{\bm{K}\bm{\psi}}$, the assembly can be written as

[TABLE]

where $m_{j}$ and $m_{j}^{e}$ denote the global and element indices of the magic points, respectively. $\bm{\psi}^{e}$ denotes the element-level projection matrix. This way, element stiffness matrices need only be computed for elements containing magic points, meaning that the hyperreduced assembly of $\widehat{\bm{K}_{m}\bm{\psi}}$ also scales roughly with $\mathcal{O}(|E_{m}|\mathcal{C}_{e}^{K})$. Furthermore, only the degrees of freedom of the $\bm{\psi}$ belonging to elements featuring magic points are required. Finally, only degrees of freedom of $\bm{u}$ belonging to elements in the reduced integration domain are required for the assembly, meaning that the required reconstruction also scales with $\mathcal{O}(|I_{m}|)$.

5.2 Discrete Empirical Interpolation Method (DEIM)

Once the magic points of the residual $\bm{g}_{m}$ have been assembled, the DEIM [42] seeks to estimate the residual $\bm{g}$ via its modes $\bm{\Omega}$ and mode coefficients $\bm{y}^{\bm{g}}$. $\bm{y}^{\bm{g}}$ can be estimated via an interpolation constraint: when reconstructing $\bm{g}$ from $\bm{y}^{\bm{g}}$, the reconstructed values at the magic points should equal the known $\bm{g}_{m}$, i.e.

[TABLE]

The estimated mode coefficient therefore become

[TABLE]

meaning that we can estimate residual and reduced residual as

[TABLE]

Consequently, the consistent hyperreduced stiffness matrix is given by

[TABLE]

In the linear case, the left factor of $\bm{K}_{m}$ and $\bm{g}_{m}$, $\bm{\psi}^{T}\bm{\Omega}(\bm{Z}^{T}\bm{\Omega})^{-1}=\bm{\psi}^{T}\bm{M}\in\mathbb{R}^{d\times m}$ can be preprocessed offline, such that the left projection only scales with $\mathcal{O}(m)$. In the nonlinear case, the product of $\bm{M}$ with the transpose of the first-stage linear compression matrix, $\bm{\bar{\varphi}}^{T}\bm{M}\in\mathbb{R}^{\bar{d}\times m}$, can similarly be preprocessed, meaning that the DEIM allows us to achieve the desired reductions in computational costs.

5.3 Linear Extrapolation Hyperreduction Method (LEHM)

We introduce a slight modification of the DEIM, which we call “Linear Extrapolation Hyperreduction Method”, i.e. LEHM, here, as this led to improvements in robustness in combination with nonlinear MOR techniques in the numerical experiments outlined below. The DEIM reconstructs $\bm{g}$ via POD modes $\bm{\Omega}$ using an interpolation constraint. Instead, the reconstruction matrix $\bm{M}$ can also be computed more directly, via a least squares problem on the snapshot data $\bm{g}_{j},j\in 1,..,s_{g}$, i.e.

[TABLE]

The solution to this is given by

[TABLE]

where $\bm{G}^{m}$ denotes the rows of the snapshots matrix corresponding to magic points only. Now, the residual and reduced residual can be approximated as

[TABLE]

Consequently, the reduced stiffness matrix is given by

[TABLE]

The linear system for the reduced increment, once approximated by extrapolation from the reduced integration domain, becomes

[TABLE]

which, just as in the case of the DEIM, can be assembled with significantly reduced computational complexity.

5.4 Least-Squares Petrov Galerkin (LSPG)

The DEIM and DEIM-like hyperreduction techniques assemble the magic points of the residual and estimate the remainder of the entries based on these via the reconstruction matrix $\bm{M}$. The hyperreduced residual is then computed by left projection onto the mode matrix $\bm{\varphi}$, meaning that these methods effectively approximate a Galerkin scheme [38, 39, 40, 41, 42, 43, 44]. They do not, however, preserve the advantageous structure of such a Galerkin scheme (e.g. symmetries), which can lead to a loss of robustness [26, 85].

Instead, the LSPG and LSPG-like schemes bypass the approximate reconstruction of the full and the reduced residual entirely [45, 46, 47, 48]. Instead, they solve a reduced equation system defined on the magic points in a weighted or least-squares manner. The LSPG simply performs an overdetermined collocation at the magic points, iteratively solving the least squares problem [45, 46, 47, 48]

[TABLE]

Standard solvers can be used to this end, with a computational cost of $\mathcal{O}(md^{2})$. The resulting LSPG algorithm effectively approximates a left basis $\bm{K}\bm{\varphi}$ via a left projection with $\bm{\varphi}^{T}\bm{K}_{m}^{T}$ being applied to the magic points, hence justifying the classification as a Petrov-Galerkin scheme.

5.5 Hyperreduced Homogenisation

Recall that, as discussed in Subsections 5.2 and 5.3, the non-magic rows of the residual and stiffness matrix can be reconstructed using a reconstruction matrix $\bm{M}$ in a DEIM-like scheme. Observing Eqs. (8) and (14), we further note that the rows of $\bm{L}$ correlate with each other just like the rows of $\bm{K}$, meaning that $\bm{M}$ can also be used to reconstruct $\bm{L}$, i.e.

[TABLE]

Consequently

[TABLE]

or, for practical computation

[TABLE]

The inversions required for computational homogenisation can thus also be computed with $\mathcal{O}(d^{3})$.

However, the computation of the macroscopic stress via Eq. (12) and the Voigt term necessary in the algorithmically consistent tangent stiffness in Eq. (13) by integration scale with $\mathcal{O}(|E|\mathcal{C}_{e}^{P})$ and $\mathcal{O}(|E|\mathcal{C}_{e}^{A})$. To reduce these computational cost, we may integrate only over elements in reduced integration domain $E_{m}$, and extrapolate from this reduced set of elements to the overall integral in an ECSW-like manner, i.e.

[TABLE]

with element coefficients ${\xi_{e}}$, at a cost of $\mathcal{O}(|E_{m}|\mathcal{C}_{e}^{P})$ [49]. The vector of element coefficients $\bm{\xi}$ can be estimated from snapshots of the homogenised stress $\bm{\bar{P}}^{vh}\in\mathbb{R}^{9s}$ and element stress $\bm{P}^{ve}\in\mathbb{R}^{9s\times|E|}$. Here, the nine entries of stress snapshots in Voigt notation are stacked underneath each other. To this end, we can minimise

[TABLE]

via standard nonlinear nonnegative least squares solver. The resulting positive weights can be interpreted as larger effective relative volumes with which each element is multiplied; each element from the reduced integration domain thus additionally represents some share of the volume neglected in integration. The Voigt stiffness term can similarly be approximated via

[TABLE]

such that its evaluations scales with $\mathcal{O}(|E_{m}|\mathcal{C}_{e}^{A})$.

With this accelerated homogenisation procedure in place, we have now hyperreduced all computationally expensive operations in Alg. 1. The resulting solution- and homogenisation-algorithm, which is outlined as pseudocode in Algs. 2 and 3, does not feature operations with scale with the original problem size, i.e. with $\mathcal{O}(D)$ or $\mathcal{O}(|E|)$, any more.

6 Numerical experiments

In this Section, we apply the MOR methods discussed above to a simple hyperelastic homogenisation problem. To this end, we consider an artificial RVE with two inclusions, as illustrated in Fig. 7. Both the matrix and the inclusions are modeled via the neo-Hookean stored energy function in Eq. (2), with Young’s moduli of $E=1000$ for the matrix and $E=3000$ for the inclusions and Poisson’s ratios of $\nu=0.2$ for both. As summarised in Tab. 1, the RVE has an edge length of $6$, while origin of the RVE coordinate system lies in a corner, such that $x_{\text{min}}=y_{\text{min}}=z_{\text{min}}=0$ and $x_{\text{max}}=y_{\text{max}}=z_{\text{max}}=6$. The inclusions each have a radius of $1.5$ and are centered at $(2,2,2)$ and $(4,4,4)$ in the RVE coordinate system, respectively. The microstructure and material behaviour result in a moderately nonlinear behaviour of the RVE. The microstructure is discretised using $4,821$ quadratic tetrahedral elements, resulting in $7,650$ nodes and $19,182$ independent degrees of freedom once boundary conditions are applied. Tab. 1 summarises further problem parameters.

Periodic boundary conditions are applied as outlined in Section 3, with the macroscopic deformation gradient $\bm{\bar{F}}$ being prescribed. We define load paths consisting of $10$ load steps each: each load step consists of a step with length $\Delta{\bar{F}}_{\text{LP}}=0.03$ along a randomly sampled load direction $\bm{N}_{\text{LP}}$ which remains constant along the whole load path, as well as a perturbation of $\Delta{\bar{F}}_{\text{LS}}=0.015$ along a direction $\bm{N}_{\text{LS}}$ which is sampled again for each step, i.e.

[TABLE]

In Fig. 6, we illustrate three components of the macroscopic deformation gradients $\bm{\bar{F}}$ generated this way. Additionally, Fig. 7 features deformation states of the RVE at the end of three example load paths, shown to scale.

For the numerical investigations, we first generate $s=200$ snapshot solutions $\bm{U}\in\mathbb{R}^{D\times s}$ along $20$ load paths via full FEM simulations. These training load paths are highlighted as red points in Fig. 6. The gradual Eigenvalue decay of the snapshot data correlation matrix, which is shown in Fig. 8, indicates that the snapshot data, and thus the solution manifold do not lie in a low-dimensional linear subspace, since otherwise, a rapid drop to zero would be observed. Meanwhile, the scale-dependent correlation dimension which is shown in Fig. 9 indicates empirically that the snapshots solutions indeed lie on a low-dimensional manifold of around $\delta=6$, since this measure approaches a value of around $8$ in the low length-scale limit.

Solution snapshots are then used to train nonlinear Galerkin ROMs using the techniques outlined in Section 4. Then, the resulting Galerkin ROMs are used to compute solutions along the same $20$ load paths to obtain residual snapshots $\bm{G}\in\mathbb{R}^{D\times s_{g}}$ at the intermediate, i.e. non-converged, steps of the Newton-Raphson solver scheme, such that $s_{g}>s$ (alternatively, residual snapshots gathered from the full FEM model could be projected onto the relevant approximation spaces post-hoc). This residual training data is used to train hyperreduced models using the techniques outlined in Section 5.

Finally, the hyperreduced nonlinear ROMs are used to compute approximate solutions $\bm{\bar{u}}_{i}^{\text{ROM}}$ along the $30$ validation load paths (alongside the $20$ training load paths). Validation solutions $\bm{u}_{i}^{\text{VAL}}$ are obtained using a full FEM model for reference. Then, the mean error in the displacement field with respect to these validation solutions

[TABLE]

is computed as a measure of ROM solver accuracy. Additionally, an equivalent error measure in the homogenised stress $\bm{\bar{P}}$ is computed as an indicator of ROM homogenisation accuracy.

The performance of all MOR techniques investigated in this work is of course subject to the choice of several algorithmic parameters. The most important of these are the reduced model size $d$ – i.e. the size of the reduced space – and the hyperreduced model size $m$ – i.e. the number of rows of $\bm{g}$ which are assembled exactly. We investigate the performance, measured in terms of accuracy, online speedup, and robustness, obtained using the MOR methods outlined in Sections 4 and 5 for a range of values of $d$ and $m$ in the following Subsection. We do not investigate the influence of further parameters here. For all methods, these were chosen in view of preliminary parameter studies such as those communicated in [121, 1]. Parameters to which all methods are subject, such as those of the Newton-Raphson solver scheme, were of course chosen to be identical for all methods. All numerical investigations in this work were performed using an in-house python FEM and MOR library.

6.1 Results

Tabs. 2, 3, 4, and 5 present the mean error over all $500$ solutions obtained using the POD, PM, and LLE in combination with DEIM-like hyperreduction techniques, for $d\in[9,3],m\in[20,250]$. High values of $E_{\text{mean}}>1\%$ are highlighted in light gray. Parameter combinations for which at least one of the $500$ simulations failed to converge are marked by a $\times$-symbol and highlighted in dark gray; no mean error was calculated in this case.

Tab. 2 makes the robustness issues of the POD-DEIM apparent; these were also highlighted by [26, 85]. For $19$ parameter combinations, at least one of the $500$ simulations fails to converge, and $7$ combinations yield unacceptable errors of $E_{\text{mean}}>1\%$. When the MOR scheme does converge to reasonable solutions, error levels are low, but this only happens for $11$ parameter combinations. Generally, higher values of $d$ and $m$ promote accuracy, but these trends are nonuniform. Additionally, there is no large region in the investigated algorithmic parameter space in which the POD-DEIM yields predictably robust results, such that it is not trivial to find parameters for which this method works as desired on the RVE problem considered here.

In combination with the POD, the LEHM modification to the DEIM outlined above does not yield significant improvements. $13$ parameter combinations yield converging simulations with reasonable results, the error is unacceptably high in $8$ cases, and in $15$ cases at least one simulation fails to converge. Error levels and trends are very similar, and the desirable regions of the algorithmic parameter space are no more contiguous.

Meanwhile, the LPOD did not yield any parameter combinations for which all $500$ simulations converged in combination with DEIM-like hyperreduction methods, suggesting significant robustness issues222When not paired with any hyperreduction scheme, no robustness issues appeared. We attempted to couple the LPOD to the DEIM and LEHM, either with multiple localised, and with one global, hyperreduction model, to no avail. That the LPOD works as desired with the LSPG (see below) suggests that the specific combination of LPOD and DEIM is responsible for this lack of convergence..

The PM, meanwhile, yields lower mean errors when all simulations converge, and there are only $4$ parameter combinations leading to excessive error levels. However, $18$ parameter combinations lead to at least one of $500$ simulations not converging; in particular, there is only one value of $m$ for which a PM model with $d=9$ runs robustly.

The results obtained by the LLE-LEHM are more encouraging. The mean error levels are lower than those obtained by all other methods – around half those obtained via the POD, and ca. $20-30\%$ lower than those obtained by the PM at many parameter combinations. While $9$ parameter combinations yield a mean error of $E_{\text{mean}}>1\%$, only $8$ parameter combinations yield one out of $500$ diverging simulations or more. Finally, there is a contiguous range of parameters in which the LLE-LEHM performs as desired in terms of accuracy and robustness on this problem, which means that it might be more straightforward to reliably select parameters for similar problems.

Next, Fig. 10 highlights the runtime consumed by different operations within the full-, reduced-, and hyperreduced FEM RVE homogenisation algorithms. Here, we only discuss results for the LLE-LEHM with $d=15$ and $m=100$, but similar observations hold for other methods and parameters. As is to be expected, the runtime of the full FEM simulation is dominated by the linear solves in the Newton-Raphson scheme, followed closely, on account of the moderate problem size $D$, by the assembly. All other operations within the Newton-Raphson scheme only account for marginal computational costs. The runtime consumed by the homogenised stress- and stiffness computations break down similarly, though they account for less of the overall runtime since they only need to be performed once after the Newton-Raphson scheme converges. The Galerkin-reduced MOR scheme eliminates the cost of the linear solver nearly completely. Assembly costs increase slightly, since the stiffness matrix and residual still need to be assembled fully and since, for some macroscopic deformation gradient values $\bm{\bar{F}}$, one more iteration is required until convergence. The hyperreduced ROM slices these assembly costs significantly, though this reduction is not as severe as that in the cost of the linear solver operations. Thus, the assembly costs dominate the runtime of the hyperreduced model.

Next, Fig. 11 highlights the tradeoff between relative runtime reductions and relative errors (in the displacement $\bm{u}$ and the homogenised stress $\bm{\bar{P}}$) which can be achieved via different $d$ and $m$ for various MOR techniques coupled with DEIM-like hyperreduction methods. The LLE-LEHM Pareto-dominates the other techniques in the trade-off between runtime and accuracy particularly in $\bm{u}$. At a fixed computational budget, the LLE can roughly halve the error obtained by the POD and achieve an improvement of ca. $20\%$ over the PM, in addition to affecting robustness improvements. For example, a $50$-fold speedup can be achieved with a $0.3\%$ error using the POD, a $0.25\%$ error using the PM, and a $0.1\%$ error using the LLE. The advantage is palpable especially at low runtimes, and speedups of $200$ with mean errors of $E_{\text{mean}}<1\%$ (in particular, $E_{\text{mean}}\approx 0.25\%$) can only be achieved using the LLE.

In the homogenised stress $\bm{\bar{P}}$, the advantage is less pronounced, as is to be expected since the averaging via Eq. (12) might mitigate local errors in the displacement. At low runtimes, there is a more significant reduction in the relative error; e.g., a $100$-fold speedup can be achieved with an error of around $0.3\%$ using the LLE and $0.5\%$ using the POD. The PM actually yields higher errors at a fixed computational budget than the POD, but this is mainly due to the PM not running robustly for smaller model sizes here.

Next, Tabs. 6, 7, 8, and 9 present the mean error $E_{\text{mean}}$ in the displacement $\bm{u}$ obtained by the POD, the LPOD, the PM and the LLE in tandem with the LSPG hyperreduction method. Again, parameter combinations for which at least one out of $500$ simulations failed to converge are highlighted in dark gray and parameter combinations yielding an excessively high mean error in light gray.

The POD performs considerably more reliably in tandem with the LSPG than with DEIM-like hyperreduction methods on the considered example. While error levels are slightly higher for equivalent parameter combinations, no simulations fail to converge and only three parameter combinations yield mean error levels of $E_{\text{mean}}>1\%$. Performance trends are still not uniform in $d$ and $m$, but are far more predictable by comparison, making for more reliable parameter setting.

In tandem with the LSPG, the LPOD yields encouraging results. There are no diverging simulations for any parameter combination, and the mean error levels are significantly below those obtained via the POD, often by about $50\%$. $9$ parameter combinations yield mean errors of $E_{\text{mean}}>1\%$ and performance trends fluctuate more strongly in $d$ and $m$, meaning that parameter setting is not entirely straightforward.

The PM produces similar mean error levels as the LPOD; sometimes being outperformed and sometimes outperforming by some margin. The PM appears to have a higher accuracy ceiling in the investigated parameter range on the considered problem and performance trends seem slightly more predictably, but at least one out of $500$ simulations fails to converge for $7$ parameter combinations, mostly at $d=9$.

The performance of the LLE-LSPG is very encouraging: on the investigated problems, no simulations fail or yield a mean error level of $E_{\text{mean}}>1\%$. Additionally, the error in the displacement field $\bm{u}$ is lower than that obtained by other techniques, and performance trends are generally predictable.

Finally, Fig. 12 highlights the tradeoff between mean errors in $\bm{u}$ and $\bm{\bar{P}}$ which can be achieved using different $d$ and $m$ using the POD, LPOD, PM, and LLE in combination with LSPG. Again, the LLE Pareto-dominates competing methods in the tradeoff between speed and accuracy in $\bm{u}$, outperforming the best of the alternatives by around $20-30\%$ in terms of relative error at a fixed computational budget. The LPOD and PM, meanwhile, perform very similarly on these measures. Again, the benefit conferred by LLE is particularly pronounced at low runtimes. A 200-fold speedup can be achieved with an error of around $0.25\%$, a 100-fold speedup with around $0.15\%$, and a 50-fold speedup with around $0.1\%$.

In terms of the error in $\bm{\bar{P}}$, the LLE slightly outperforms the POD and PM, but the LPOD Pareto-dominates here. The more continuously nonlinear approximation spaces obtained by the LLE and the PM thus do not seem to confer any advantage in the computation of the homogenised stress on the considered example. The observation made in the context of DEIM-like hyperreduction is accentuated here: continuously nonlinear approximation spaces yield an advantage when the underlying solution field is of interest, but this advantage is diminished when the homogenised quantities are more relevant.

7 Summary and Outlook

In a recent work [1], we proposed a projection-based nonlinear MOR scheme which uses graph-based manifold learning techniques to obtain flexible nonlinear approximation spaces. In this work, we show how this approach can be employed to reduce computational costs by multiple orders of magnitude while retaining high levels of accuracy. To this end, we extended the nonlinear MOR method with two hyperreduction methods: a DEIM-like approach and LSPG. Additionally, we improved the robustness and performance of the local online linearisation with an approximate Euler-backward-like scheme and a two-stage approach for the DEIM-like hyperreduction. Finally, the NLMOR scheme was extended to the homogenisation of stresses and stiffnesses based on RVE solutions. The resulting, hyperreduced algorithm no longer scales with the sizes $D$ and $|E|$ of the original problem.

On the example problem considered above, the hyperreduced manifold learning approach proposed in this work facilitates speedups of two orders of magnitude with negligible mean validation errors of $E_{\text{mean}}\approx 0.1\%$ while requiring only $s=200$ snapshots. When used in tandem with a DEIM-like hyperreduction scheme, the graph-based manifold learning approach Pareto-dominates alternative approaches in terms of the tradeoff between runtime and accuracy in the predicted displacement $\bm{u}$ and homogenised stress $\bm{\bar{P}}$. The advantage over competing methods is more pronounced in the displacement predictions, yielding a benefit of around $50\%$ over the POD and $20-30\%$ over the PM in the investigated RVE problem. The LPOD faces severe robustness issues, with at least one out of $500$ simulations diverging for every investigated algorithmic parameter combination. The POD and PM work as desired for a comparatively narrow and unpredictable parameter range, while the LLE yields good results for a somewhat broader, contiguous region in the algorithmic parameter space.

When used with LSPG, all methods perform much more robustly, though the mean error levels obtained with this method are slightly higher. The POD, LPOD, and PM yield unacceptably high mean errors or fail to converge only for $3$, $9$, and $7$ parameter combinations only, while the LLE never does. Again, the LLE Pareto-dominates all other methods in the tradeoff between speed and accuracy in the displacement $\bm{u}$, outperforming the best competing method by around $20-30\%$. In terms of the homogenised stress $\bm{\bar{P}}$, however, the LLE actually yields worse results than the LPOD.

For the example RVE considered here, nonlinear MOR techniques based on graph-based manifold learning methods such as LLE thus yield an advantage in the tradeoff between accuracy and speed especially when the underlying solution field is of interest, while the benefit to using such methods for the computation of the homogenised stress is less pronounced. It would be intriguing to investigate whether this observation generalises to other, more complex microstructures with larger phase contrasts and more nonlinear material behaviour. A systematic investigation over a range of RVEs could help identify where nonlinear approximation spaces might be employed profitably, and which methods to generate nonlinear approximation spaces might be most suitable for a particular class of multiscale problems.

Additionally, continuously and flexibly nonlinear approximation spaces might be of interest for history-dependent or multiphysical multiscale problems involving e.g. plasticity or damage indicators for which the evolution of variables on the microscale is critical. Of course, further developments are necessary to facilitate the application of these nonlinear MOR techniques to coupled and history-dependent problems, particularly if these involve localisation phenomena. A promising direction for research might be to generalise the manifold learning approaches explored in this work to multiple manifolds of coupled reduced solution variables.

If these applications to more complicated problem classes prove promising, it would be sensible to attempt to push the nonlinear MOR techniques closer to their performance and development ceilings with further theoretical and implementational work, especially in terms of computational cost. With this in mind, it is important to note that the assembly, rather than linear system solutions, constitute the bottleneck in terms of runtime in the hyperreduced model. Consequently, efforts toward further improving the speed of hyperreduced models ought to be targeted mainly at reducing $m$ rather than $d$. This to some extent motivates the considerable body of research effort expended on developing more advanced hyperreduction techniques. Note, however, that nonlinear Galerkin-MOR techniques, which work with smaller approximation spaces, might be able to yield improvements here, too. Smaller approximation spaces result in smaller reduced residual manifolds; and ideally fewer mesh entities $m$ being required to parameterise them.

Acknowledgements

We would like to thank Rudy Geelen for his helpful comments.