# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE', which is part of this source code package.
#
#
import os
try:
from mpi4py import MPI
except ImportError:# pragma: no cover
print("MPI capabilities not available")
import numpy as np
from scipy import sparse
from Mordicus.BasicAlgorithms import SVD as SVD
[docs]
def CompressData(
collectionProblemData, solutionName, tolerance = None, snapshotCorrelationOperator = None, snapshots = None, compressSolutions = False, nbModes = None
):
"""
Computes a reducedOrderBasis using the SnapshotPOD algorithm, from the
snapshots contained in the iterator snapshotsIterator, which a correlation
operator between the snapshots defined by the matrix
snapshotCorrelationOperator, with tolerance as target accuracy of the data
compression. If a reducedOrderBasis prexists, it is updated using the
snapshots from the solutions.
Parameters
----------
collectionProblemData : CollectionProblemData
input collectionProblemData containing the data
solutionName : str
name of the solutions from which snapshots are taken
tolerance : float, cannot be provided with nbModes
target accuracy of the data compression
snapshotCorrelationOperator : scipy.sparse.csr_matrix, optional
correlation operator between the snapshots
snapshots : np.ndarray, optional
of size (nbSnapshots, numberOfDofs): snapshots of the solutions to
compress
compressSolutions : bool, optional
True to compresse solutions using the constructed reducedOrderBasis
nbModes : int, cannot be provided with tolerance
the number of keps snapshot POD modes
"""
assert isinstance(solutionName, str)
if tolerance == None and nbModes == None:# pragma: no cover
raise("must specify epsilon or nbModes")
if tolerance != None and nbModes != None:# pragma: no cover
raise("cannot specify both epsilon and nbModes")
if snapshots is None:
snapshots = collectionProblemData.GetSnapshots(solutionName)
if snapshotCorrelationOperator is None:
snapshotCorrelationOperator = sparse.eye(snapshots.shape[1])
numberOfSnapshots = snapshots.shape[0]
previousReducedOrderBasis = collectionProblemData.GetReducedOrderBasis(solutionName)
correlationMatrix = np.zeros((numberOfSnapshots,numberOfSnapshots))
for i, snapshot1 in enumerate(snapshots):
matVecProduct = snapshotCorrelationOperator.dot(snapshot1)
for j, snapshot2 in enumerate(snapshots):
if j <= i and j < numberOfSnapshots:
correlationMatrix[i, j] = np.dot(matVecProduct, snapshot2)
mpiReducedCorrelationMatrix = np.zeros((numberOfSnapshots, numberOfSnapshots))
MPI.COMM_WORLD.Allreduce([correlationMatrix, MPI.DOUBLE], [mpiReducedCorrelationMatrix, MPI.DOUBLE])
if tolerance != None:
eigenValuesRed, eigenVectorsRed = SVD.TruncatedSVDSymLower(mpiReducedCorrelationMatrix, tolerance)
else:
eigenValuesRed, eigenVectorsRed = SVD.TruncatedSVDSymLower(mpiReducedCorrelationMatrix, nbModes = nbModes)
nbePODModes = eigenValuesRed.shape[0]
changeOfBasisMatrix = np.zeros((nbePODModes,numberOfSnapshots))
for j in range(nbePODModes):
changeOfBasisMatrix[j,:] = eigenVectorsRed[:,j]/np.sqrt(eigenValuesRed[j])
reducedOrderBasis = np.dot(changeOfBasisMatrix,snapshots)
if previousReducedOrderBasis is None:
collectionProblemData.AddReducedOrderBasis(solutionName, reducedOrderBasis)
else:
print("detecting existing POD basis")
snapshots = np.append(previousReducedOrderBasis, reducedOrderBasis, axis=0)
numberOfSnapshots = snapshots.shape[0]
correlationMatrix = np.zeros((numberOfSnapshots,numberOfSnapshots))
for i, snapshot1 in enumerate(snapshots):
matVecProduct = snapshotCorrelationOperator.dot(snapshot1)
for j, snapshot2 in enumerate(snapshots):
if j <= i and j < numberOfSnapshots:
correlationMatrix[i, j] = np.dot(matVecProduct, snapshot2)
mpiReducedCorrelationMatrix = np.zeros((numberOfSnapshots, numberOfSnapshots))
MPI.COMM_WORLD.Allreduce([correlationMatrix, MPI.DOUBLE], [mpiReducedCorrelationMatrix, MPI.DOUBLE])
if tolerance != None:
eigenValuesRed, eigenVectorsRed = SVD.TruncatedSVDSymLower(mpiReducedCorrelationMatrix, tolerance)
else:
eigenValuesRed, eigenVectorsRed = SVD.TruncatedSVDSymLower(mpiReducedCorrelationMatrix, nbModes = nbModes)
nbePODModes = eigenValuesRed.shape[0]
changeOfBasisMatrix = np.zeros((nbePODModes,numberOfSnapshots))
for j in range(nbePODModes):
changeOfBasisMatrix[j,:] = eigenVectorsRed[:,j]/np.sqrt(eigenValuesRed[j])
reducedOrderBasis = np.dot(changeOfBasisMatrix,snapshots)
collectionProblemData.AddReducedOrderBasis(solutionName, reducedOrderBasis)
if compressSolutions == True:
collectionProblemData.CompressSolutions(solutionName, snapshotCorrelationOperator)
print("POD for "+solutionName+" : number of snapshots = "+str(numberOfSnapshots)+", number of modes = "+str(nbePODModes))
if __name__ == "__main__":# pragma: no cover
from genericROM import RunTestFile
RunTestFile(__file__)