Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/rrompy/hfengines/scipy_engines/oscillator_problem_engine.py
6c4eee04e0dfmaster

Last Change
Raw File

oscillator_problem_engine.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
import scipy.sparse as scsp
from numbers import Number
from rrompy.hfengines.base.linear_affine_engine import LinearAffineEngine
from rrompy.hfengines.base.scipy_engine_base import ScipyEngineBaseTensorized
from rrompy.utilities.base.types import List, Np1D, Np2D
from rrompy.utilities.exception_manager import RROMPyAssert
__all__ = ['OscillatorProblemEngine', 'AugmentedOscillatorProblemEngine']
class OscillatorProblemEngine(LinearAffineEngine, ScipyEngineBaseTensorized):
"""
Solver for damped mass oscillator problems.
Each i represents the factor in front of mu_i, except for i=0, which is a
constant. For each i:
(mass_i)_j = mass of j-th mass
(damping_i)_jk = damping between j-th and k-th mass
(stiffness_i)_jk = stiffness between j-th and k-th mass
Others:
input_ij = stiffness between i-th input and j-th mass
output_ij = contribution of j-th mass's position to i-th output
"""
def __init__(self, mass:List[Np1D], damping:List[Np2D],
stiffness:List[Np2D], input:Np2D, output:Np2D,
inputPosition : bool = True, verbosity : int = 10,
timestamp : bool = True):
super().__init__(verbosity = verbosity, timestamp = timestamp)
self._affinePoly = True
self.npar = len(mass)
RROMPyAssert(len(damping), self.npar, "Number of parameters")
RROMPyAssert(len(stiffness), self.npar, "Number of parameters")
input = np.array(input)
if input.ndim < 2: input = input.reshape(-1, 1)
self.nports = input.shape[1]
self.nAs, self.As, self.thAs = 3 * self.npar, [], []
self.nbs, self.bs = 1, [input.flatten()]
self._C = output
for j in range(self.npar):
deg = [0] * (self.npar - 1)
if j != 0: deg[j - 1] = 1
K = - stiffness[j]
Dj = - 1.j * damping[j]
Kcorr = - np.array(np.sum(K, axis = 1)).flatten()
if j == 0 and inputPosition:
Kcorr += np.array(np.sum(input, axis = 1)).flatten()
Djcorr = - np.array(np.sum(Dj, axis = 1)).flatten()
if isinstance(Dj, (np.ndarray,)):
Mm = np.diag(- mass[j])
K = K + np.diag(Kcorr - np.diag(K))
Dj = Dj + np.diag(Djcorr - np.diag(Dj))
else:
Mm = scsp.diags(- mass[j], 0, format = Dj.format)
K = K + scsp.diags([Kcorr - K.diagonal()], [0])
Dj = Dj + scsp.diags([Djcorr - Dj.diagonal()], [0])
self.As += [K, Dj, Mm]
self.thAs += [self.getMonomialSingleWeight([0] + deg),
self.getMonomialSingleWeight([1] + deg),
self.getMonomialSingleWeight([2] + deg)]
if isinstance(self.As[0], (np.ndarray,)):
self.setSolver("SOLVE")
else:
self.setSolver("SPSOLVE", {"use_umfpack" : False})
class AugmentedOscillatorProblemEngine(OscillatorProblemEngine):
"""
Solver for damped mass oscillator problems.
Each i represents the factor in front of mu_i, except for i=0, which is a
constant. For each i:
(mass_i)_j = mass of j-th mass
(damping_i)_jk = damping between j-th and k-th mass
(stiffness_i)_jk = stiffness between j-th and k-th mass
Others:
input_ij = stiffness between i-th input and j-th mass
output_ij = contribution of j-th mass's position to i-th output
"""
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
N, self._nAs = self.As[0].shape[1], 2 * self.npar
self.bs = [np.pad(self.bs[0].reshape(-1, self.nports),
[(N, 0), (0, 0)], 'constant').flatten()]
if (not isinstance(self._C, (Number,))
and self._C.shape[1] != 2 * N):
RROMPyAssert(self._C.shape[1], N, "Shape of output")
self._C = np.pad(self._C, [(0, 0), (0, N)],
'constant')
As, thAs = [], []
for j in range(self.npar):
K, Dj, Mm = self.As[3 * j], self.As[3 * j + 1], self.As[3 * j + 2]
if isinstance(Dj, (np.ndarray,)):
I, Z = 1.j * np.eye(N), np.zeros((N, N))
Deff = np.block([[Z, - I], [K, Dj]])
Meff = np.block([[I, Z], [Z, Mm]])
else:
I = 1.j * scsp.eye(N)
Deff = scsp.bmat([[None, - I], [K, Dj]], format = Mm.format)
Meff = scsp.block_diag([I, Mm], format = Mm.format)
As += [Deff, Meff]
thAs += [self.thAs[3 * j], self.thAs[3 * j + 1]]
self.As, self.thAs = As, thAs
if isinstance(self.As[0], (np.ndarray,)):
self.setSolver("SOLVE")
else:
self.setSolver("SPSOLVE", {"use_umfpack" : False})
c4science · Help