romtools.hyper_reduction.ecsw

  1#
  2# ************************************************************************
  3#
  4#                         ROM Tools and Workflows
  5# Copyright 2019 National Technology & Engineering Solutions of Sandia,LLC
  6#                              (NTESS)
  7#
  8# Under the terms of Contract DE-NA0003525 with NTESS, the
  9# U.S. Government retains certain rights in this software.
 10#
 11# ROM Tools and Workflows is licensed under BSD-3-Clause terms of use:
 12#
 13# Redistribution and use in source and binary forms, with or without
 14# modification, are permitted provided that the following conditions
 15# are met:
 16#
 17# 1. Redistributions of source code must retain the above copyright
 18# notice, this list of conditions and the following disclaimer.
 19#
 20# 2. Redistributions in binary form must reproduce the above copyright
 21# notice, this list of conditions and the following disclaimer in the
 22# documentation and/or other materials provided with the distribution.
 23#
 24# 3. Neither the name of the copyright holder nor the names of its
 25# contributors may be used to endorse or promote products derived
 26# from this software without specific prior written permission.
 27#
 28# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 29# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 30# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 31# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 32# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 33# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 34# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 35# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 36# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 37# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
 38# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 39# POSSIBILITY OF SUCH DAMAGE.
 40#
 41# Questions? Contact Eric Parish (ejparis@sandia.gov)
 42#
 43# ************************************************************************
 44#
 45
 46'''
 47Implementation of Energy-conserving sampling and weighting (ECSW) hyper-reduction
 48
 49Energy-conserving sampling and weighting (ECSW) is a hyper-reduction approach
 50originally developed specifically for solid mechanics problems, but it has
 51since been generalized. It is a project-then-approximate hyper-reduction
 52approach, similar in spirit and implementation to empirical cubature
 53approaches. The name comes from the energy conservation properties the method
 54has for solid mechanics applications; note that this property is not provable
 55for general systems.
 56
 57Given a set of residual snapshots ECSW computes sampling indices $GID_i$ and
 58weights $\\xi_i$. Specifically, the residual snapshots must be computed for
 59reconstructed full-order model snapshots,
 60
 61$$\\boldsymbol r_i = \\boldsymbol r( \\boldsymbol \\Phi \\boldsymbol \\Phi^T (\\mathbf{u}_i - \\mathbf{u}_0))$$
 62
 63where $\\boldsymbol r_i$ is the residual $i$th residual snapshot,
 64$\\mathbf{u}_i$ is the $i$th state snapshot, and $\\boldsymbol \\Phi$ is the
 65trial basis.
 66
 67The goal of ECSW is to find a sparse set of weights to approximate the reduced
 68residual with a subset of local test-basis/residual products
 69
 70$$\\sum_{e \\in \\mathcal{E}} \\xi_e \\boldsymbol \\Psi_e^T \\boldsymbol r_e \\approx \\Psi^T \\boldsymbol r$$
 71
 72For more details, consult Chapman et al. 2016 DOI: 10.1002/nme.5332.
 73
 74The ECSW class contains the methods needed to compute sampling indices and
 75weights given a set of residual snapshot and trial basis data.
 76'''
 77
 78import sys
 79import abc
 80from typing import Tuple
 81import numpy as np
 82import romtools.linalg.linalg as la
 83
 84
 85class ECSWsolver(abc.ABC):
 86    '''
 87    Abstract base class for ECSW solvers
 88
 89    ECSW solvers should take in a linear system constructed from projected residual vector
 90    snapshots and the contributions at each mesh degree of freedom to the projected snapshot.
 91    The solvers should return arrays with sample mesh indices and weights.
 92
 93    Methods:
 94    '''
 95
 96    @abc.abstractmethod
 97    def __init__(self, solver_param_dict: dict = None):
 98        '''
 99        Set solver parameters to non-default values
100
101        Args:
102            (optional) solver_param_dict: dictionary, with some of the following keys:
103            max_iters: int, maximum overall iterations
104            max_non_neg_iters: int, maximum inner iterations to enforce non-negativity
105            max_iters_res_unchanged: int, maximum number of iterations without any change in the residual norm before terminating
106            zero_tol: int, tolerance used to check if weights or residual norm changes are near zero
107        '''
108        pass
109
110    @abc.abstractmethod
111    def __call__(
112        self, full_mesh_lhs: np.ndarray, full_mesh_rhs: np.array, tolerance: np.double
113    ) -> Tuple[np.ndarray, np.ndarray]:
114        '''
115        Compute the sample mesh DoF indices and corresponding weights
116
117        Args:
118            full_mesh_lhs: (n_snap*n_rom, n_dof) numpy ndarray, where n_snap is the number of residual snapshots,
119                                                                n_rom is the ROM dimension,
120            and n_dof is the number of mesh degrees of freedom (DoFs) (nodes, volumes, or elements)
121            full_mesh_rhs: (n_snap*n_rom,) numpy array
122            tolerance: Double, the ECSW tolerance parameter. Lower values of tolerance will result in more mesh DoF samples
123
124        Returns:
125            Tuple of numpy ndarrays.
126            First array: (n_dof_sample_mesh,) numpy ndarray of ints, the mesh indices in the sample mesh.
127            Second array: (n_dof_sample_mesh,) numpy ndarray of doubles, the corresponding sample mesh weights.
128
129        '''
130        pass
131
132
133# ecsw non-negative least-squares class
134# TODO: scipy loss and enhancement method for testing? Won't work in place of
135#       this function because tau can't be specified
136
137
138class ECSWsolverNNLS(ECSWsolver):
139    '''
140    Given a linear system with left-hand side full_mesh_lhs and right-hand side full_mesh_rhs,
141    compute sample mesh indices and weights for ECSW using the non-negative least squares algorithm
142    from Chapman et al. 2016 DOI: 10.1002/nme.5332.
143    '''
144
145    def __init__(self, solver_param_dict: dict = None):
146        # Set default solver parameter values
147        self.max_iters = 10000  # maximum overall iterations
148        self.max_non_neg_iters = (
149            100  # maximum inner iterations to enforce non-negativity
150        )
151        self.max_iters_res_unchanged = 10  # maximum number of iterations without any change in residual before terminating
152        self.zero_tol = 1e-12  # tolerance used to check if weights or residual norm changes are near zero
153
154        if solver_param_dict is not None:
155            if "max_iters" in solver_param_dict.keys():
156                self.max_iters = solver_param_dict["max_iters"]
157            if "max_non_neg_iters" in solver_param_dict.keys():
158                self.max_non_neg_iters = solver_param_dict["max_non_neg_iters"]
159            if "max_iters_res_unchanged" in solver_param_dict.keys():
160                self.max_iters_res_unchanged = solver_param_dict[
161                    "max_iters_res_unchanged"
162                ]
163            if "zero_tol" in solver_param_dict.keys():
164                self.zero_tol = solver_param_dict["zero_tol"]
165
166    def __call__(
167        self, full_mesh_lhs: np.ndarray, full_mesh_rhs: np.array, tolerance: np.double
168    ):
169        '''
170        Compute the sample mesh DoF indices and corresponding weights using the non-negative least squares algorithm
171        from Chapman et al. 2016 DOI: 10.1002/nme.5332.
172
173        min || full_mesh_lhs*full_mesh_weights-full_mesh_rhs ||_2^2 s.t. full_mesh_weights >=0, ||
174        full_mesh_lhs*full_mesh_weights-full_mesh_rhs ||_2 < tolerance ||full_mesh_rhs||_2
175
176        Args:
177            full_mesh_lhs: (n_snap*n_rom, n_dof) numpy ndarray, where
178                                n_snap is the number of residual snapshots,
179                                n_rom is the ROM dimension, and
180                                n_dof is the number of mesh degrees of freedom (DoFs) (nodes, volumes, or elements)
181            full_mesh_rhs: (n_snap*n_rom,) numpy array
182            tolerance: Double, the ECSW tolerance parameter. Lower values of tolerance will result in more mesh DoF samples
183
184        Returns:
185            Tuple of numpy ndarrays.
186            First array: (n_dof_sample_mesh,) numpy ndarray of ints, the mesh indices in the sample mesh.
187            Second array: (n_dof_sample_mesh,) numpy ndarray of doubles, the corresponding sample mesh weights.
188
189        '''
190
191        n_dof = full_mesh_lhs.shape[1]
192
193        # initialize
194        iters = 0
195        sample_mesh_indicies = []
196        residual = full_mesh_rhs.copy()
197        residual_norm = np.linalg.norm(residual)
198        target_norm = tolerance * residual_norm
199
200        full_mesh_weights = np.zeros(n_dof)
201        full_mesh_candidate_weights = np.zeros(n_dof)
202
203        residual_norm_unchanged = 0
204
205        # add nodes to sample mesh until tolerance is met
206        while (
207            (residual_norm > target_norm)
208            and (residual_norm_unchanged < self.max_iters_res_unchanged)
209            and (iters < self.max_iters)
210        ):
211            # determine new node to add to sample mesh
212            weighted_residual = np.dot(full_mesh_lhs.T, residual)
213
214            # make sure mesh entity hasn't already been selected
215            still_searching = True
216            if len(sample_mesh_indicies) == n_dof:
217                still_searching = False
218            while still_searching:
219                mesh_index = la.argmax(weighted_residual)
220                still_searching = False
221                if mesh_index in sample_mesh_indicies:
222                    still_searching = True
223                    weighted_residual[mesh_index] = la.min(weighted_residual)
224                    if np.all(
225                        (weighted_residual - weighted_residual[mesh_index]) < 1e-15
226                    ):
227                        # All elements are the same, select random element outside of sample mesh
228                        unselected_indicies = np.setdiff1d(
229                            np.arange(n_dof), sample_mesh_indicies, assume_unique=True
230                        )
231                        mesh_index = unselected_indicies[0]
232                        still_searching = False
233
234            # add new mesh entity index
235            if len(sample_mesh_indicies) < n_dof:
236                sample_mesh_indicies.append(mesh_index)
237
238            print(
239                f"iteration={iters}  sample mesh size={len(sample_mesh_indicies)}  "
240                f"residual norm={residual_norm:.8e}  ratio to target={residual_norm / target_norm:.8e}"
241            )
242
243            iters += 1
244
245            # compute corresponding weights
246
247            for j in range(self.max_non_neg_iters):
248                sample_mesh_lhs = full_mesh_lhs[:, sample_mesh_indicies]
249                sample_mesh_candidate_weights = np.dot(
250                    np.linalg.pinv(sample_mesh_lhs), full_mesh_rhs
251                )
252                full_mesh_candidate_weights *= 0
253                full_mesh_candidate_weights[sample_mesh_indicies] = (
254                    sample_mesh_candidate_weights
255                )
256                if np.all(sample_mesh_candidate_weights > 0):
257                    full_mesh_weights = full_mesh_candidate_weights.copy()
258                    break
259                # line search to enforce non-negativity
260                max_step = self.__max_feasible_step(
261                    full_mesh_weights[sample_mesh_indicies],
262                    sample_mesh_candidate_weights,
263                )
264                full_mesh_weights_new = full_mesh_weights + max_step * (
265                    full_mesh_candidate_weights - full_mesh_weights
266                )
267
268                # remove zero valued indices
269                near_zero_inds = np.nonzero(
270                    full_mesh_weights_new[sample_mesh_indicies] < self.zero_tol
271                )[0]
272                samp_inds_for_removal = [
273                    sample_mesh_indicies[i] for i in near_zero_inds
274                ]
275                for samp_ind in samp_inds_for_removal:
276                    sample_mesh_indicies.remove(samp_ind)
277
278                # increment iteration count
279                print(
280                    f"iteration={iters}  sample mesh size={len(sample_mesh_indicies)}  "
281                    f"residual norm={residual_norm:.8e}  ratio to target={residual_norm / target_norm:.8e}"
282                )
283                iters += 1
284                full_mesh_weights = 1 * full_mesh_weights_new
285
286            if j == self.max_non_neg_iters - 1:
287                sys.exit("Error: NNLS algorithm failed to compute weights")
288
289            # update least-squares residual
290            sample_mesh_weights = full_mesh_weights[sample_mesh_indicies]
291            residual = full_mesh_rhs - np.dot(sample_mesh_lhs, sample_mesh_weights)
292            residul_old_norm = 1 * residual_norm
293            residual_norm = np.linalg.norm(residual)
294
295            if np.abs(residual_norm - residul_old_norm) < self.zero_tol:
296                residual_norm_unchanged += 1
297            else:
298                residual_norm_unchanged = 0
299
300            if residual_norm_unchanged >= self.max_iters_res_unchanged:
301                print(
302                    f"WARNING: Norm has not changed more than {self.zero_tol} "
303                    f"in {self.max_iters_res_unchanged} steps, exiting NNLS."
304                )
305
306        print("NNLS complete! Final stats:")
307        print(
308            f"iteration={iters}  sample mesh size={len(sample_mesh_indicies)}  "
309            f"residual norm={residual_norm:.8e}  ratio to target={residual_norm / target_norm:.8e}"
310        )
311
312        return np.array(sample_mesh_indicies, dtype=int), sample_mesh_weights
313
314    def __max_feasible_step(self, weights, candidate_weights):
315        '''
316        determine maximum update step size such that:
317        weights + step_size * (candidate_weights-weights) >=0
318
319        Args:
320            weights: (n,) array
321            candidate_weights: (n, array)
322
323        Returns:
324            step_size: double
325        '''
326        inds = np.argwhere(candidate_weights <= 0)
327        step_size = 1.0
328        for i in inds:
329            if weights[i] == 0.0:
330                step_size = 0
331            else:
332                step_size_i = weights[i] / (weights[i] - candidate_weights[i])
333                step_size = min([step_size, step_size_i])
334        return step_size
335
336
337# ESCW helper functions for specific test basis types
338
339
340def _construct_linear_system(
341    residual_snapshots: np.ndarray, test_basis: np.ndarray, n_var: int
342):
343    '''
344    Construct the linear system required for ECSW with a fixed test basis, such as POD-Galerkin projection.
345
346    Args:
347        residual_snapshots: (n_var, n_dof, n_snap) numpy ndarray, where
348                                n_dof is the number of mesh degrees of freedom (DoFs) (nodes, volumes, or elements),
349                                n_var is the number of residual variables, and
350                                n_snap is the number of snapshots
351        test_basis: (n_var, n_dof, n_mode) numpy ndarray, where n_mode is the number of modes in the basis.
352        n_var: int, the number of residual variables (e.g. for fluid flow, residual variable could be mass,
353                                                           x-momentum, y-momentum, z-momentum, and energy)
354
355    Returns:
356        full_mesh_lhs: (n_snap*n_mode, n_dof) numpy ndarray, the left-hand side of the linear system required by the ECSW solver
357        full_mesh_rhs: (n_snap*n_rom,) numpy array, the right-hand side of the linear system required by the ECSW solver
358    '''
359
360    (n_var, n_dof, n_snap) = residual_snapshots.shape
361    (n_var_tb, n_dof_tb, n_mode) = test_basis.shape
362    assert n_var == n_var_tb
363    assert n_dof == n_dof_tb
364    # construct ECSW system
365    full_mesh_lhs = np.zeros((n_snap * n_mode, n_dof))
366
367    # left-hand side
368    for i in range(n_dof):
369        # should be projection of all variables for a given mesh DoF
370        phi_block = test_basis[:, i, :]  # n_var x n_mode
371        res_snaps_block = residual_snapshots[:, i, :]  # n_var x n_snap
372        full_mesh_lhs_block = np.dot(
373            phi_block.T, res_snaps_block
374        )  # n_modes x n_snaps matrix
375        full_mesh_lhs[:, i] = np.ravel(full_mesh_lhs_block, order="F")
376
377    # right-hand-side
378    full_mesh_rhs = np.sum(full_mesh_lhs, axis=1)
379
380    return full_mesh_lhs, full_mesh_rhs
381
382
383def ecsw_fixed_test_basis(
384    ecsw_solver: ECSWsolver,
385    residual_snapshots: np.ndarray,
386    test_basis: np.ndarray,
387    n_var: int,
388    tolerance: np.double,
389):
390    '''
391    ECSW implementation for a fixed test basis, such as POD-Galerkin projection
392
393    Args:
394        ecsw_solver: ECSWsolver object corresponding to a child class with concrete implementations such as ECSWsolverNNLS.
395        residual_snapshots: (n_var, n_dof, n_snap) numpy ndarray, where
396                                n_dof is the number of mesh degrees of freedom (DoFs) (nodes, volumes, or elements),
397                                n_var is the number of residual variables, and
398                                n_snap is the number of snapshots
399        test_basis: (n_var, n_dof, n_mode) numpy ndarray, where n_mode is the number of modes in the basis.
400        n_var: int, the number of residual variables (e.g. for fluid flow, residual variable could be mass,
401                                                      x-momentum, y-momentum, z-momentum, and energy)
402        tolerance: Double, the ECSW tolerance parameter. Lower values of tolerance will result in more mesh DoF samples
403
404    Returns:
405        Tuple of numpy ndarrays.
406        First array: (n_dof_sample_mesh,) numpy ndarray of ints, the mesh indices in the sample mesh.
407        Second array: (n_dof_sample_mesh,) numpy ndarray of doubles, the corresponding sample mesh weights.
408    '''
409
410    # TODO need to incorporate residual scales here too, perhaps using scaler.py
411    full_mesh_lhs, full_mesh_rhs = _construct_linear_system(
412        residual_snapshots, test_basis, n_var
413    )
414
415    return ecsw_solver(full_mesh_lhs, full_mesh_rhs, tolerance)
416
417
418def ecsw_varying_test_basis(
419    ecsw_solver: ECSWsolver,
420    residual_snapshots: np.ndarray,
421    test_basis: np.ndarray,
422    n_var: int,
423    tolerance: np.double,
424):
425    '''
426    ECSW implementation for a varying test basis, such as Least-Squares Petrov-Galerkin projection
427
428    Args:
429        ecsw_solver: ECSWsolver object corresponding to a child class with concrete implementations such as ECSWsolverNNLS.
430        residual_snapshots: (n_var, n_dof, n_snap) numpy ndarray, where
431                                n_dof is the number of mesh degrees of freedom (DoFs) (nodes, volumes, or elements),
432                                n_var is the number of residual variables,
433                                and n_snap is the number of snapshots
434        test_basis: (n_snap, n_var, n_dof, n_mode) numpy ndarray, where
435                                n_snap is the number of test basis snapshots and
436                                n_mode is the number of modes in the basis.
437        n_var: int, the number of residual variables (e.g. for fluid flow, residual variable could be mass,
438                                                      x-momentum, y-momentum, z-momentum, and energy)
439        tolerance: Double, the ECSW tolerance parameter. Lower values of tolerance will result in more mesh DoF samples
440
441    Returns:
442        Tuple of numpy ndarrays.
443        First array: (n_dof_sample_mesh,) numpy ndarray of ints, the mesh indices in the sample mesh.
444        Second array: (n_dof_sample_mesh,) numpy ndarray of doubles, the corresponding sample mesh weights.
445    '''
446
447    (n_var, n_dof, n_snap) = residual_snapshots.shape
448    (n_snap_tb, n_var_tb, n_dof_tb, _) = test_basis.shape
449    assert n_var == n_var_tb
450    assert n_dof == n_dof_tb
451    assert n_snap == n_snap_tb
452
453    # Construct sequence of fixed test basis matrices, then stack them
454    lhs_list = []
455    rhs_list = []
456    for i in range(n_snap):
457        # TODO need to incorporate residual scales here too, perhaps using scaler.py
458        test_basis_i = test_basis[i]
459        residual_snapshot_i = residual_snapshots[:,:,i]
460        residual_snapshot_i = residual_snapshot_i[:, :, np.newaxis]
461        full_mesh_lhs, full_mesh_rhs = _construct_linear_system(
462            residual_snapshot_i, test_basis_i, n_var
463        )
464        lhs_list.append(full_mesh_lhs)
465        rhs_list.append(full_mesh_rhs)
466
467    full_mesh_lhs = np.concatenate(lhs_list,axis=0)
468    full_mesh_rhs = np.concatenate(rhs_list,axis=0)
469
470    return ecsw_solver(full_mesh_lhs, full_mesh_rhs, tolerance)
471
472
473def ecsw_lspg_zero_residual(
474    ecsw_solver: ECSWsolver,
475    test_basis: np.ndarray,
476    n_var: int,
477    tolerance: np.double,
478):
479    '''
480    ECSW implementation for Least-Squares Petrov-Galerkin projection of a discrete system with near-zero residual snapshots
481    Uses the test basis as a snapshot instead, as proposed in section 3.2 of
482    https://www.sciencedirect.com/science/article/pii/S0045782522004558.
483
484    Args:
485        ecsw_solver: ECSWsolver object corresponding to a child class with concrete implementations such as ECSWsolverNNLS.
486        test_basis: (n_snap, n_var, n_dof, n_mode) numpy ndarray, where
487                        n_snap is the number of test basis snapshots and
488                        n_mode is the number of modes in the basis.
489        n_var: int, the number of residual variables (e.g. for fluid flow, residual variable could be mass,
490                                                      x-momentum, y-momentum, z-momentum, and energy)
491        tolerance: Double, the ECSW tolerance parameter. Lower values of tolerance will result in more mesh DoF samples
492
493    Returns:
494        Tuple of numpy ndarrays.
495        First array: (n_dof_sample_mesh,) numpy ndarray of ints, the mesh indices in the sample mesh.
496        Second array: (n_dof_sample_mesh,) numpy ndarray of doubles, the corresponding sample mesh weights.
497    '''
498
499    (n_snap, n_var, _, _) = test_basis.shape
500
501    # Construct sequence of fixed test basis matrices, then stack them
502    lhs_list = []
503    rhs_list = []
504    for i in range(n_snap):
505        # TODO need to incorporate residual scales here too, perhaps using scaler.py
506        test_basis_i = test_basis[i]
507        full_mesh_lhs, full_mesh_rhs = _construct_linear_system(
508            test_basis_i, test_basis_i, n_var
509        )
510        lhs_list.append(full_mesh_lhs)
511        rhs_list.append(full_mesh_rhs)
512
513    full_mesh_lhs = np.concatenate(lhs_list,axis=0)
514    full_mesh_rhs = np.concatenate(rhs_list,axis=0)
515
516    return ecsw_solver(full_mesh_lhs, full_mesh_rhs, tolerance)