A hyper-reduction strategy for elastoplastic hydro-mechanical problems in partially saturated soil

1 Introduction

• Repeated parameter evaluations are needed for dam monitoring and decision support, while each coupled simulation requires nonlinear solution across hundreds of time steps.The motivation includes data assimilation, optimization, control, and uncertainty quantification.Source: §1 · pp. 1–2
• The target problem combines elastoplastic soil deformation and partially saturated flow, with strong nonlinearities arising simultaneously in both governing equations.The paper identifies this class of hyper-reduction as previously unaddressed to the authors’ knowledge.Source: §1 · p. 2
• Standard finite-element assembly in FEniCSx favors element-based quadrature; ECSW’s non-negative weights preserve Galerkin-projection structure and support robustness.The method requires no constitutive-formulation-specific modification or additional full-order simulations beyond those used to construct the POD basis.Source: §1 · p. 3
• Yield zones and regions near the phreatic surface localize different nonlinear behaviors, while mechanical and flow residuals can differ in magnitude and nonlinear content.A cubature trained on the complete residual would also spend selected elements reproducing linear terms that can instead be computed exactly.Source: §1 · p. 3
• Field-wise POD bases and separate ECSW rules approximate nonlinear corrections relative to exact affine baselines, reducing the reduced-mesh size at equal accuracy.The mechanical baseline is an elastic stress response; the hydraulic baseline freezes retention coefficients at a mean-pressure state.Source: §1 · p. 3
• 7×–77× speed-ups are reported across the experiments while hyper-reduced solutions remain within the reduced-basis model’s accuracy range.The tests simulate two- and three-dimensional dams through a raise and subsequent consolidation, varying friction angle, Young’s modulus, saturated conductivity, and raise height.Source: §1 · p. 3

2 Coupled elastoplastic hydro-mechanical model

• The model solves for solid displacement and pore-water pressure, with pressure positive in saturated soil and negative under unsaturated suction.It assumes constant atmospheric gas pressure, neglects mass transfer and inertia, and treats solid grains as incompressible.Source: §2 · p. 4

2.1 Momentum balance equation and mechanical constitutive model

• Drucker–Prager elastoplasticity governs stress from elastic strain, while the method treats the stress law as a generic history-dependent function and does not depend on this constitutive choice.The test formulation assumes perfect plasticity and uses friction angle and cohesion to set model parameters.Source: §2.1 · pp. 4–5

2.2 Mass balance equation and retention curve

• The van Genuchten retention model makes saturation and water content pressure-dependent, and hydraulic conductivity varies with effective saturation, producing strongly nonlinear flow.The water balance also couples to mechanics through the rate of volumetric strain acting as a source.Source: §2.2 · p. 5

2.3 Initial and boundary conditions

• Each field has its own Dirichlet, Neumann, and Robin boundary partitions; mechanical Robin conditions model elastic ground support, while hydraulic Robin conditions describe seepage.The coupled initial state is obtained by solving the steady-state problem under the initial boundary data.Source: §2.3 · p. 6

3.1 Weak form and time discretization

• At each time step, the fully implicit weak form couples displacement and pressure through history-dependent elastoplastic stress and pressure-dependent hydraulic coefficients.Backward Euler evaluates the pressure-dependent coefficients at the new pressure.Source: §3.1 · pp. 6–7

3.2 Spatial discretization

• The spatial discretization uses quadratic displacement and linear pressure Taylor–Hood elements, producing a nonlinear algebraic system assembled from element contributions.The residual and Jacobian inherit the sum-over-elements structure used by the hyper-reduction.Source: §3.2 · p. 7

3.3 Newton–Raphson solution of the time step

• A monolithic Newton–Raphson iteration updates displacement and pressure using the residual Jacobian and local elastoplastic return mapping.The iteration stops when relative residual and both fields’ relative increments fall below 10^-6.Source: §3.3 · pp. 7–8
• Repeated full-mesh residual and Jacobian assembly, including constitutive evaluations at every quadrature point, is a central cost of the full-order solve.Solving the sparse linear system is the other operation identified as a major cost.Source: §3.3 · p. 8

4 Reduced-order model

• The HR-RB model combines field-split POD and ECSW, using the same full-order snapshots as the RB model without requiring additional simulations.It uses one ECSW rule per governing equation, trained on the nonlinear correction relative to an affine baseline.Source: §4 · p. 8

4.1 Field-split POD reduced basis

• Separate displacement and pressure POD bases are built from centered snapshots, avoiding cross-field normalization of quantities with different magnitudes.The retained modes are selected using the POD truncation tolerance.Source: §4.1 · pp. 8–9

4.2 Decomposed ECSW with an affine baseline

• Independent ECSW rules for the mechanical and flow equations can select different, non-overlapping elements because their nonlinearities act in different mesh regions.A common rule would evaluate both equations on elements where only one nonlinearity acts.Source: §4.2 · pp. 9–10
• The affine baseline is projected exactly, leaving ECSW to approximate only the nonlinear correction and confining cubature error to that correction.The baseline uses linear elastic stress and pressure-dependent coefficients evaluated at mean training pressure; its projection is precomputed and reduced to a small matrix-vector product online.Source: §4.2 · p. 10

4.3 Element selection algorithm

• A greedy sparse NNLS rule selects elements to fit correction residuals on training snapshots, stopping at tolerance τ to trade accuracy for sparsity.The offline QR approach reduces the selection problem’s memory burden, although its triangular factor grows with the square of the element count.Source: §4.3 · pp. 11–12

• Both examples simulate dam raising as a distributed load on the crest and upstream face, followed by consolidation.Source: §5 · p. 12
• The two-dimensional plane-strain example is used to study hyper-reduction behavior under mesh refinement.Source: §5 · p. 12
• The fully three-dimensional valley-dam example tests whether the conclusions carry over to a full-order model an order of magnitude larger.Source: §5 · p. 12

5.1.1 Problem setup

• The 10 m-high embankment is simulated through a raise and subsequent consolidation, with a 10-day raise followed by 20 days of consolidation.Four ROM inputs vary: dam friction angle, Young’s modulus, saturated conductivity, and raise height.Source: §5.1.1 · pp. 12, 14, §5.1.2 · p. 14

5.1.2 Parameters and training set

• Forty Latin-hypercube parameter points provide 12,520 snapshots per field, with ten additional points reserved for testing; training cases span 3%–77% plastified quadrature points.The largest plastic zones occur with low friction angle and large raise height.Source: §5.1.2 · p. 15

5.1.3 Reduced bases and reduced meshes

• 34 unknowns replace 4,763–45,907; as mesh elements rise tenfold, ECSW selections grow only from 48 to 55 mechanically and 64 to 75 for flow.The selected regions remain tied to yielding and desaturation zones rather than mesh resolution.Source: §5.1.3 · pp. 16–17, Figure 4 · p. 16
• With matched accuracy, omitting the affine baseline grows the flow selection from 64 to 176 elements, while the mechanical selection changes from 48 to 52.Without the baseline, flow elements spread across the dam to represent the mean-state response.Source: §5.1.3 · p. 17, Table 4 · p. 17

5.1.4 Accuracy and speed-up

• About 0.1% global error and 30–43× speed-ups on the finest mesh; maximum errors stay below 0.2%, with RB and HR-RB results essentially identical.Errors remain around 10^-3 through raising and consolidation.Source: §5.1.4 · pp. 18–19, Table 5 · p. 18, Table 6 · p. 20
• Speed-ups reach 7–9×, 14–19×, and 30–43× from coarsest to finest mesh, and increase with the extent of plastic deformation.Hyper-reduction cuts return mapping, assembly, and projection costs; reduced-model Newton iterations remain within 16% of the full-order count.Source: §5.1.4 · pp. 19–20

5.2.1 Problem setup

• The 50 m crest dam sits in a V-shaped valley, with compliant valley support and hydrostatic pressure on the wet upstream face.The downstream face is traction-free and a seepage face; remaining boundaries are impervious.Source: §5.2.1 · p. 21
• Two half-height layers are deposited in five-day sweeps along the dam axis; inclined fronts smooth loading before consolidation.The 30-day transient begins from the steady state before the raise.Source: §5.2.1 · p. 21
• The 117,883-tetrahedron P2–P1 model has 551,116 unknowns, and training cases span 5–39% plastified quadrature points.The valley walls support part of the load, so yielding is less extensive than in plane strain.Source: §5.2.1 · p. 21, Table 7 · p. 21

5.2.2 Reduced-order results

• The 110-unknown reduced problem replaces a full-order model with 551,116 unknowns.The POD bases use 34 displacement modes and 76 pressure modes.Source: §5.2.2 · p. 23
• The mechanical and flow rules select 110 and 611 of 117,883 elements, concentrating near valley contacts and the base, and along the phreatic surface.Their locations are shown in Figures 9 and 10.Source: §5.2.2 · p. 23, Figure 9 · p. 24, Figure 10 · p. 24
• Global errors are about 0.3% for displacement and 0.5% for pressure; per-step errors stay below 1.9%.The HR-RB errors sit on the RB error floor, so cubature adds no visible error.Source: §5.2.2 · p. 23, Table 8 · p. 23
• 65–77× speed-ups cut three-dimensional runs from 11–14 hours to 10–11 minutes, versus 2.0–2.2× for RB.The HR-RB times include precomputation of the reduced baseline operators.Source: §5.2.2 · p. 23
• With refinement, hyper-reduced costs grow more slowly than full-order costs because the reduced dimensions and cubature sizes track physics rather than mesh size.This makes the benefit largest where the full-order model is most expensive.Source: §5.2.2 · pp. 23–24

6 Conclusions

• The method uses separate displacement and pressure POD bases and ECSW rules, each trained on its residual's nonlinear correction relative to an exact affine baseline.The construction uses the POD training snapshots and assembles original variational forms on selected subdomains in FEniCSx.Source: §6 · p. 24
• In 2D, speed-ups rise from 7–9× on the coarsest mesh to 30–43× on the finest, with errors around 0.1%.The HR-RB errors sit on the RB error floor, so accuracy is limited by basis truncation rather than cubature.Source: §6 · p. 26
• The selected reduced meshes track plastified and unsaturated regions whose extent is set by physics, while full-order cost grows with mesh size.The paper describes this spatial selection as an adaptive strategy that permits uniform meshes fine enough to resolve the physics.Source: §6 · p. 26
• In 3D, displacement and pressure errors are about 0.3% and 0.5%, while speed-ups reach 65–77×, reducing 11–14-hour runs to 10–11 minutes.The 3D dam model is an order of magnitude larger and has a fully three-dimensional deposition response.Source: Table 9 · p. 25, §6 · p. 26
• The construction does not depend on the specific mechanical law, retention curve, or relative-permeability model.The employed Drucker–Prager, van Genuchten, and cubic conductivity models enter through assembled residuals.Source: §6 · p. 26

A Blocks of the Jacobian

• At each time step, the coupled problem is solved for displacement and pressure, using the previous step's state as known data.The Jacobian blocks follow by differentiating the residuals with respect to the two unknowns.Source: A Blocks of the Jacobian · p. 26
• The Jacobian blocks are bilinear forms in displacement and pressure increments and their corresponding test functions.All coefficients are evaluated at the current Newton–Raphson iterate.Source: A Blocks of the Jacobian · p. 26
• The seepage term is restricted to saturated parts of the seepage boundary, while Ctan is the consistent tangent of the return mapping.A prime denotes differentiation with respect to pressure.Source: A Blocks of the Jacobian · p. 27