For projection-based ROMs:

is a "large" (POD) subspace always an "enemy"?

Sometimes, it can be your friend...

RAMSES - Dec. 2021 SISSA, Italy

This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. This work was funded by the Advanced Simulation and Computing program and the Laboratory Directed Research and Development program at Sandia National Laboratories, a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525 

Francesco Rizzi

(NexGen Analytics)

(in collaboration with E. Parish, P. Blonigan, J.Tencer - Sandia National Labs)

Follow along/slides at:

fnrizzi.github.io/ramses-12-2021

Why projection-based ROMs?

  • A priori training data and a projection process
    • ​Leverages the governing equations
    • Explainability, physics is already there
    • Error bounds and models
  • To consider:
    • Intrusivity, hyper-reduction for nonlinear models 
    • Linear time-invariant (LTI) systems: theory is quite mature
  • Various:
    • Interplay between DL and pROMs?
    • Can pROMs benefit from DL? Definitely, but must be practical

Follow along/slides at: fnrizzi.github.io/ramses-12-2021/

 

  • Historically, a key focus of pROMs work has been:
    "finding the smallest subspace that can represent/solve a problem"
     

  • Intuitively            : small system, more convenient to compute

  • Mathematically  : intriguing but hard

  • Computationally: is this really the best approach?

    • What if we can formulate the problem such that we don't need to reduce it so much while being efficient?
       

  • This talk aims to provide a counterargument

Why this talk?

Follow along/slides at: fnrizzi.github.io/ramses-12-2021/

  • Focus on pROMs for LTI systems

    • ​Weren't they a solved problem...?

    • Emphasis on computational aspects​

    • ​Little math, error bounds, ML, deep learning (sorry!)

  • ​Disclaimer: this work might seem "obvious" (depending on whom you are talking to)

  • Format: this is going to be a "story"

    • Walk through how this work started and developed

    • Finally, we talk about generalization

Content and format

Follow along/slides at: fnrizzi.github.io/ramses-12-2021/

About a year ago...

  • Started reading about seismic waves

  • ​​Obviously a very important topic, no pROMs work on this

    • ​Relevant: pROMs for acoustic waves: V. Pereyra et al., Electron. Trans. Numer. Anal. (2008)

  • Surface waves: travel at the Earth's surface

  • Body waves: travel through the Earth

    • Affected by the material properties (density, modulus)

    • Primary (P-waves): compressional, longitudinal

    • Secondary or shear (S-waves): transversal (particles oscillate perpendicularly to the direction of wave propagation)

  • Interesting topic, so we started on this

 Formulation

  • Approximate the Earth (or generic planet) as a sphere
  • Simplify via axisymmetric approximation 

Surface

\partial / \partial \phi = 0
\partial / \partial \phi = 0

 Elastic shear waves

  • Velocity-stress axisymmetric formulation*:
\newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \small{ \dens \dvpdt = \dsrpdr + \frac{1}{r} \dstpdtheta + \frac{1}{r} \left(3 \srp + 2 \stp \cot{\theta} \right) + \force }
\newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \small{ \dstpdt = \shMod \left( \frac{1}{r}\dvpdtheta - \frac{\cot{\theta}}{r}\vp \right) }
\newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \small{ \dsrpdt = \shMod \left( \dvpdr - \frac{1}{r}\vp \right) }

  • Shear effects negligible in liquid, the core is not considered

Surface

Core-mantle

\newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} f(r, \theta, t): \text{forcing}
\newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \shMod(r, \theta): \text{shear modulus}
\newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \dens(r, \theta): \text{density}
\newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \vp (r, \theta, t): \text{velocity}
\newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \srp(r,\theta,t), \stp(r,\theta,t): \text{stresses}

Given:

Find:

* H. Igel, M. Weber,  Geophys. Res. Lett. 22 (6) (1995)  

Semi-discrete system

%% FOM \newcommand{\state}{\boldsymbol x} \newcommand{\systemMat}{\mathbf{A}} \newcommand{\forcing}{\boldsymbol f} \newcommand{\fomDim}{N} \newcommand{\tfinal}{t_{final}} \newcommand{\paramsA}{\boldsymbol \eta} \newcommand{\paramsF}{\boldsymbol \mu} \newcommand{\paramDomainA}{\mathcal D_{\eta}} \newcommand{\paramDomainF}{\mathcal D_{\mu}} \newcommand{\nParamsA}{N_{\paramsA}} \newcommand{\nParamsF}{N_{\paramsF}} \newcommand{\paramsMat}{ \mathcal{M} } \newcommand{\stateRef}{\state_\text{ref}(\paramsA, \paramsF)} \newcommand{\stateRefNoPars}{\state_\text{ref}} %% ROM \newcommand{\romBasis}{\boldsymbol{\Phi}} \newcommand{\romState}{\hat{\state}} \newcommand{\romStateTensor}{\hat{\stateTensor}} \newcommand{\stateTensorRef}{{\stateTensor}_{\text{ref}}} \newcommand{\stateTensorRefWithPars}{{\stateTensor}_{\text{ref}}(\paramsA, \paramsMat)} \newcommand{\stateTensor}{\boldsymbol X} \newcommand{\approxStateTensor}{\tilde{\boldsymbol X}} \newcommand{\romSystemMat}{\hat{\systemMat}} \newcommand{\romDim}{K} \newcommand{\nRuns}{M} \newcommand{\forcingTensor}{\boldsymbol F} \newcommand{\approxState}{\tilde{\state}} %\newcommand{\stateRef}{{\state}_{\text{ref}}} \newcommand{\romBasisCol}{\boldsymbol \phi} \newcommand{\trialSubspace}{\mathcal{V}} \newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \frac{d \textcolor{red}{\state_{\vp}(t)}}{dt} = \systemMat_{\vp} \textcolor{gray}{\state_{\stresses}(t)} + \forcing_{\vp}(t)
%% FOM \newcommand{\state}{\boldsymbol x} \newcommand{\systemMat}{\mathbf{A}} \newcommand{\forcing}{\boldsymbol f} \newcommand{\fomDim}{N} \newcommand{\tfinal}{t_{final}} \newcommand{\paramsA}{\boldsymbol \eta} \newcommand{\paramsF}{\boldsymbol \mu} \newcommand{\paramDomainA}{\mathcal D_{\eta}} \newcommand{\paramDomainF}{\mathcal D_{\mu}} \newcommand{\nParamsA}{N_{\paramsA}} \newcommand{\nParamsF}{N_{\paramsF}} \newcommand{\paramsMat}{ \mathcal{M} } \newcommand{\stateRef}{\state_\text{ref}(\paramsA, \paramsF)} \newcommand{\stateRefNoPars}{\state_\text{ref}} %% ROM \newcommand{\romBasis}{\boldsymbol{\Phi}} \newcommand{\romState}{\hat{\state}} \newcommand{\romStateTensor}{\hat{\stateTensor}} \newcommand{\stateTensorRef}{{\stateTensor}_{\text{ref}}} \newcommand{\stateTensorRefWithPars}{{\stateTensor}_{\text{ref}}(\paramsA, \paramsMat)} \newcommand{\stateTensor}{\boldsymbol X} \newcommand{\approxStateTensor}{\tilde{\boldsymbol X}} \newcommand{\romSystemMat}{\hat{\systemMat}} \newcommand{\romDim}{K} \newcommand{\nRuns}{M} \newcommand{\forcingTensor}{\boldsymbol F} \newcommand{\approxState}{\tilde{\state}} %\newcommand{\stateRef}{{\state}_{\text{ref}}} \newcommand{\romBasisCol}{\boldsymbol \phi} \newcommand{\trialSubspace}{\mathcal{V}} \newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \frac{d \textcolor{gray}{\state_{\stresses}(t)}}{dt} = \systemMat_{\stresses} \textcolor{red}{\state_{\vp}(t)}
  • Staggered grid discretization, 2nd-order centered differences
  • Coupled system:

Sparse coefficient matrices

(depend on material prop, not on time)

Velocities

Stresses

  • BC: reflective on surface and core-mantle boundary, zero stress on symmetry axes
  • This is a linear time-invariant system: good! LTI systems* are well-known and studied

* P. Benner, S. Gugercin, K. Willcox, SIAM, 2015

LTI but complex dynamics

Contour plots of velocity field: Ricker wavelet source T = 60 sec, depth = 640 Km

Interference

Reflection

Refraction (from discontinuities)

PREM Earth model

time = 250 sec

time = 1000 sec

time = 2000 sec

%% FOM \newcommand{\state}{\boldsymbol x} \newcommand{\systemMat}{\mathbf{A}} \newcommand{\forcing}{\boldsymbol f} \newcommand{\fomDim}{N} \newcommand{\tfinal}{t_{final}} \newcommand{\paramsA}{\boldsymbol \eta} \newcommand{\paramsF}{\boldsymbol \mu} \newcommand{\paramDomainA}{\mathcal D_{\eta}} \newcommand{\paramDomainF}{\mathcal D_{\mu}} \newcommand{\nParamsA}{N_{\paramsA}} \newcommand{\nParamsF}{N_{\paramsF}} \newcommand{\paramsMat}{ \mathcal{M} } \newcommand{\stateRef}{\state_\text{ref}(\paramsA, \paramsF)} \newcommand{\stateRefNoPars}{\state_\text{ref}} %% ROM \newcommand{\romBasis}{\boldsymbol{\Phi}} \newcommand{\romState}{\hat{\state}} \newcommand{\romStateTensor}{\hat{\stateTensor}} \newcommand{\stateTensorRef}{{\stateTensor}_{\text{ref}}} \newcommand{\stateTensorRefWithPars}{{\stateTensor}_{\text{ref}}(\paramsA, \paramsMat)} \newcommand{\stateTensor}{\boldsymbol X} \newcommand{\approxStateTensor}{\tilde{\boldsymbol X}} \newcommand{\romSystemMat}{\hat{\systemMat}} \newcommand{\romDim}{K} \newcommand{\nRuns}{M} \newcommand{\forcingTensor}{\boldsymbol F} \newcommand{\approxState}{\tilde{\state}} %\newcommand{\stateRef}{{\state}_{\text{ref}}} \newcommand{\romBasisCol}{\boldsymbol \phi} \newcommand{\trialSubspace}{\mathcal{V}} \newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \frac{d \romState_{\vp}(t)}{dt} = \romBasis_{\vp}^T \systemMat_{\vp} \romBasis_{\stresses} \romState_{\stresses}(t) + \romBasis_{\vp}^T \forcing_{\vp}(t)
%% FOM \newcommand{\state}{\boldsymbol x} \newcommand{\systemMat}{\mathbf{A}} \newcommand{\forcing}{\boldsymbol f} \newcommand{\fomDim}{N} \newcommand{\tfinal}{t_{final}} \newcommand{\paramsA}{\boldsymbol \eta} \newcommand{\paramsF}{\boldsymbol \mu} \newcommand{\paramDomainA}{\mathcal D_{\eta}} \newcommand{\paramDomainF}{\mathcal D_{\mu}} \newcommand{\nParamsA}{N_{\paramsA}} \newcommand{\nParamsF}{N_{\paramsF}} \newcommand{\paramsMat}{ \mathcal{M} } \newcommand{\stateRef}{\state_\text{ref}(\paramsA, \paramsF)} \newcommand{\stateRefNoPars}{\state_\text{ref}} %% ROM \newcommand{\romBasis}{\boldsymbol{\Phi}} \newcommand{\romState}{\hat{\state}} \newcommand{\romStateTensor}{\hat{\stateTensor}} \newcommand{\stateTensorRef}{{\stateTensor}_{\text{ref}}} \newcommand{\stateTensorRefWithPars}{{\stateTensor}_{\text{ref}}(\paramsA, \paramsMat)} \newcommand{\stateTensor}{\boldsymbol X} \newcommand{\approxStateTensor}{\tilde{\boldsymbol X}} \newcommand{\romSystemMat}{\hat{\systemMat}} \newcommand{\romDim}{K} \newcommand{\nRuns}{M} \newcommand{\forcingTensor}{\boldsymbol F} \newcommand{\approxState}{\tilde{\state}} %\newcommand{\stateRef}{{\state}_{\text{ref}}} \newcommand{\romBasisCol}{\boldsymbol \phi} \newcommand{\trialSubspace}{\mathcal{V}} \newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \frac{d \romState_{\stresses}(t)}{dt} = \romBasis_{\stresses}^T \systemMat_{\stresses} \romBasis_{\vp} \romState_{\vp}(t)

Galerkin ROM formulation

For simplicity, assume same # of modes = K for velocity/stresses

 

Typical form for LTI systems, e.g.:
P. Benner, S. Gugercin, K. Willcox, SIAM, 2015

%% FOM \newcommand{\state}{\boldsymbol x} \newcommand{\systemMat}{\mathbf{A}} \newcommand{\forcing}{\boldsymbol f} \newcommand{\fomDim}{N} \newcommand{\tfinal}{t_{final}} \newcommand{\paramsA}{\boldsymbol \eta} \newcommand{\paramsF}{\boldsymbol \mu} \newcommand{\paramDomainA}{\mathcal D_{\eta}} \newcommand{\paramDomainF}{\mathcal D_{\mu}} \newcommand{\nParamsA}{N_{\paramsA}} \newcommand{\nParamsF}{N_{\paramsF}} \newcommand{\paramsMat}{ \mathcal{M} } \newcommand{\stateRef}{\state_\text{ref}(\paramsA, \paramsF)} \newcommand{\stateRefNoPars}{\state_\text{ref}} %% ROM \newcommand{\romBasis}{\boldsymbol{\Phi}} \newcommand{\romState}{\hat{\state}} \newcommand{\romStateTensor}{\hat{\stateTensor}} \newcommand{\stateTensorRef}{{\stateTensor}_{\text{ref}}} \newcommand{\stateTensorRefWithPars}{{\stateTensor}_{\text{ref}}(\paramsA, \paramsMat)} \newcommand{\stateTensor}{\boldsymbol X} \newcommand{\approxStateTensor}{\tilde{\boldsymbol X}} \newcommand{\romSystemMat}{\hat{\systemMat}} \newcommand{\romDim}{K} \newcommand{\nRuns}{M} \newcommand{\forcingTensor}{\boldsymbol F} \newcommand{\approxState}{\tilde{\state}} %\newcommand{\stateRef}{{\state}_{\text{ref}}} \newcommand{\romBasisCol}{\boldsymbol \phi} \newcommand{\trialSubspace}{\mathcal{V}} \newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \state_{\vp}(t) \approx \romBasis_{\vp} \hat{\state}_{\vp}(t),
%% FOM \newcommand{\state}{\boldsymbol x} \newcommand{\systemMat}{\mathbf{A}} \newcommand{\forcing}{\boldsymbol f} \newcommand{\fomDim}{N} \newcommand{\tfinal}{t_{final}} \newcommand{\paramsA}{\boldsymbol \eta} \newcommand{\paramsF}{\boldsymbol \mu} \newcommand{\paramDomainA}{\mathcal D_{\eta}} \newcommand{\paramDomainF}{\mathcal D_{\mu}} \newcommand{\nParamsA}{N_{\paramsA}} \newcommand{\nParamsF}{N_{\paramsF}} \newcommand{\paramsMat}{ \mathcal{M} } \newcommand{\stateRef}{\state_\text{ref}(\paramsA, \paramsF)} \newcommand{\stateRefNoPars}{\state_\text{ref}} %% ROM \newcommand{\romBasis}{\boldsymbol{\Phi}} \newcommand{\romState}{\hat{\state}} \newcommand{\romStateTensor}{\hat{\stateTensor}} \newcommand{\stateTensorRef}{{\stateTensor}_{\text{ref}}} \newcommand{\stateTensorRefWithPars}{{\stateTensor}_{\text{ref}}(\paramsA, \paramsMat)} \newcommand{\stateTensor}{\boldsymbol X} \newcommand{\approxStateTensor}{\tilde{\boldsymbol X}} \newcommand{\romSystemMat}{\hat{\systemMat}} \newcommand{\romDim}{K} \newcommand{\nRuns}{M} \newcommand{\forcingTensor}{\boldsymbol F} \newcommand{\approxState}{\tilde{\state}} %\newcommand{\stateRef}{{\state}_{\text{ref}}} \newcommand{\romBasisCol}{\boldsymbol \phi} \newcommand{\trialSubspace}{\mathcal{V}} \newcommand{\spherCoords}{r, \theta, \phi} \newcommand{\axiCoords}{r, \theta} \newcommand{\coords}{r, \theta} \newcommand{\depVars}{t, \coords} \newcommand{\dens}{\rho} \newcommand{\shMod}{G} \newcommand{\stresses}{\sigma} \newcommand{\vp}{v} \newcommand{\dvpdt}{ \frac{\partial \vp}{\partial t} } \newcommand{\dvpdr}{ \frac{\partial \vp}{\partial r} } \newcommand{\dvpdtheta}{ \frac{\partial \vp}{\partial \theta} } \newcommand{\spp}{\sigma_{\phi\phi}} \newcommand{\dsppdp}{ \frac{\partial \spp}{\partial \phi} } \newcommand{\srp}{\sigma_{r\phi}} \newcommand{\dsrpdt}{ \frac{\partial \srp}{\partial t} } \newcommand{\dsrpdr}{ \frac{\partial \srp}{\partial r} } \newcommand{\stp}{\sigma_{\theta\phi}} \newcommand{\dstpdt}{ \frac{\partial \stp}{\partial t} } \newcommand{\dstpdtheta}{ \frac{\partial \stp}{\partial \theta} } \newcommand{\force}{f} \state_{\stresses}(t) \approx \romBasis_{\stresses} \hat{\state}_{\stresses}(t)

Approximation:

Predictive ROM: quantitative

  • Earth, total simulation 2000 (sec)
  • Forcing: Ricker wavelet: depth = 150 Km, period T ∈ [31, 69] (sec) 
  • Relative L2 norm of the error over domain at final time

Many modes!
Not surprising!

Velocity field at final time = 2000 secs computed for the T=69 (extrapolation point)

ROM using 436 modes for velocity and 417 modes for stresses

Predictive ROM: qualitative

Observations

  • Key point: we need many modes or ROM is useless
  • Speedup: a single ROM is 10X faster
  • O(1000) reduction in dofs, O(10) in runtime... something wrong?
    • You might say: "limited" speedup due to 850 modes.
  • Is 850 too much or is it too much for this formulation?
  • High-fidelity: ~3,150,000 spatial dofs
  • ROM: ~850 modes (velocity+stresses)
    Many modes needed, this is expected.

~3700X reduction in # dofs!

Galerkin ROM: Analysis

  • Dense system
     
  • General matrix-vector product (gemv) kernel
     
  • BLAS level-2 kernel

How do we evaluate the efficiency of this kernel?

Arithmetic Intensity (ai)

\text{arithmetic intensity} = \frac{\text{minimal FLOP-count}}{\text{minimal data movement}}

  • Assuming a square system of size     for simplicity

  • Data movement:       (read     ) +       (read       ) +     (write    )

  • FLOPS:

  • ai ~1/4      (hint: this is small)

  • gemv kernel:

y = A x + b
K^2
2K
K
K
2K^2
A
x,b
y

Roofline Model

  • Visual performance model obtained by plotting:

    • ​performance (in GFLOPs/s) against their arithmetic intensity

    • two lines to indicate the theoretically attainable performance depending on the arithmetic intensity​

  • Evaluate resource efficiency by relating its algorithm’s arithmetic intensity relative to the hardware’s peak main-memory bandwidth and floating-point performance

  • Hardware limitations for a given kernel, prioritize optimizations​

1/4

Memory bandwidth bound!

​theoretically attainable performance depending on the arithmetic intensity​

  • BLAS level-1 traces back to 1979
  • Level-2 kernels matrix-vector in 1988
  • Level-3 matrix-matrix in 1990

Memory BW-bound: so what? 

Hardware has changed from the 80's!

Strive for compute-bound

  • MBB kernels: not ideal for modern many cores arch

  • Best when:

    • ​cores are kept busy, data is local

    • access patterns are optimal for the targeted arch

Standard Galerkin ROM 

Can we improve Galerkin ROMs for LTI?

Galerkin Rank-2 formulation

  • This is useful when we need many solves

  • Let's consider M trajectories simultaneously
    e.g. different forcing evaluations

Arithmetic Intensity

\text{arithmetic intensity} = \frac{\text{minimal FLOP-count}}{\text{minimal data movement}}

~ K / 16

Rank-2: why is it good?

Standard formulation

Rank-2 formulation

Results: performance

Kokkos implementation with OpenMP backend;
workstation with two 18-core Intel(R) Xeon(R) Gold 6154 CPU @ 3.00 GHz (24.75MB L3 cache, 125GB total mem)*

  • M=1: very limited

  • M>1: increasing # of threads helps

  • A large K and M is an advantage!
    Allows to fully exploit the machine!

  • M = # of simultaneous trajectories

M = 1  : standard Galerkin ROM

M >= 2: rank-2 ROM formulation

* F.Rizzi et al, CMAME, 2021

When to prefer the Rank-2ROM?

What combination of thread count (n) and number of trajectories, M, would be the most efficient to obtain those P samples while satisfying the given constraints?

  • Suppose: we have a single compute node with 36 cores
  • Want P (large) trajectories (e.g. from forcing realizations)
  • 1024 is the maximum feasible # of trajectories that can simultaneously be processed on the node
    • e.g. due to memory constraints (collecting data, etc)

Launch 36 single-thread ROM runs each using M=1

and repeat until all my P samples are done

CPU 0

CPU 1

Launch 18 two-threaded ROM runs each using M=1

and repeat until all my P samples are done

If M changes...

If we increase # of modes (K), things improve!

# of modes (K) = 512

# of modes (K) = 2048

Compute speedup wrt standard Galerkin ROM (M=1)

Greener is better

Take home message

For K = 256  : rank2-ROM is 13X more efficient than rank-1

For K = 512  : 19X

For K = 1024: 23X

For K = 2048: 26X 

The larger the number of modes, the more efficient it is to evaluate an ensemble of trajectories!

Applicable to other LTI systems

  • Aeroelasticity:
    d    eforming structures modeled as linear, with a nonlinear load

  • Acoustic waves:
    modeled with a linear PDE, but can have a number of nonlinear sources (turbulent shear layers from wakes)

  • Neutral particle (neutron, photon, etc.) transport

  • Linear circuit models

Challenges/outlook

  • What about if the matrix A changes?

  • What about nonlinear problems?

  • Tensors are getting more attention due to DeepLearning

    • ROMs for LTI can benefit from them

    • Leverage hardware evolution: CUDA tensor cores

Rank-3 Galerkin ROM?

Conclusions/references

  • (Sometimes) many modes are necessary for accuracy and provide a computational advantage
  • pROMs can benefit from more work, not just theoretically
    • Previous formulations can (should) be revisited
  • pROMs to efficiently generate data to build other surrogates (ML)
  • A compute-bound formulation of Galerkin model reduction for linear time-invariant dynamical systems, F.Rizzi, E.Parish, P.Blonigan, J.Tencer, CMAME, 2021
  • V. Pereyra et al., Electron. Trans. Numer. Anal. Vol. 30 (2008) 406–41
  • E. Phipps et al, ArXiv, 2015
  • M. Gunzburger et al., SIAM J. Numer. Anal. 55 (1) (2017) 286–304
  • M. Gunzburger et al., IMA J. Numer. Anal. 39 (3) (2018) 1180–1205

Thanks!

 

Wave code is open-source: https://github.com/Pressio/SHAW

(within our Pressio ecosystem for ROMs)

(and part of the Exascale Computing Project [ECP] Proxy apps suite)

 

Questions? 

 

Happy to talk offline too:
francesco.rizzi@ng-analytics.com

fnrizzi@sandia.gov

Slides at: fnrizzi.github.io/ramses-12-2021/

Backup

Why surrogates?

  • If you are here today, likely you use and/or study and/or believe in surrogate modeling. So I could spend minutes on this but...
     

  • Computing/hardware progresses and changes quickly

    • Exascale is already here: China has two machines

    • How does this impact surrogates (if at all)?

    • Can we/how to leverage this for surrogate modeling?

      • ​"It allows me to run my same old surrogate faster": not ideal!
         

  • More synergistic development of surrogates and computing?

Monte Carlo study

MC study: 512 trajectories sampling the forcing period T

Rank-2 ROM is 950 times faster than FOM

Rank-2 ROM is 15 times faster than rank-1 ROM

Source: https://www.alcf.anl.gov/files/DMello-Nguyen-ALCF-CP-Workshop-MKL-2019-05-01-2019.pdf