Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/rrompy/sampling/scipy/sampling_engine_arnoldi.py
51cd9a0a331afeature_greedy_coarseni…

Last Change
Raw File

sampling_engine_arnoldi.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/>.
#
from copy import copy
import numpy as np
from rrompy.sampling.base.pod_engine import PODEngine
from .sampling_engine_krylov import SamplingEngineKrylov
from rrompy.utilities.base.types import Np1D
from rrompy.utilities.base import verbosityDepth
__all__ = ['SamplingEngineArnoldi']
class SamplingEngineArnoldi(SamplingEngineKrylov):
"""HERE"""
def resetHistory(self):
super().resetHistory()
self.HArnoldi = None
self.RArnoldi = None
self.RHSs = None
self.samplesAug = None
def popSample(self):
if hasattr(self, "nsamples") and self.nsamples > 1:
self.HArnoldi = self.HArnoldi[: -1, : -1]
self.RArnoldi = self.RArnoldi[: -1, : -1]
if self.nsamples > 2:
self.RHSs = self.RHSs[:, : -1]
else:
self.RHSs = None
self.samplesAug = self.RHSs[self.HFEngine.V.dim() :, : -1]
super().popSample()
@property
def HFEngine(self):
"""Value of HFEngine. Its assignment resets history."""
return self._HFEngine
@HFEngine.setter
def HFEngine(self, HFEngine):
self._HFEngine = HFEngine
self.resetHistory()
self.PODEngine = PODEngine(self._HFEngine)
def preprocesssamples(self):
ns = self.nsamples
if ns <= 0: return
return self.samplesAug[:, ns - 1].reshape((-1,self.HFEngine.V.dim())).T
def preprocessb(self, mu:complex, overwrite : bool = False,
homogeneized : bool = False):
ns = self.nsamples
r = super().preprocessb(mu, overwrite, homogeneized)
if ns == 0:
return r
elif ns == 1:
r = r / self.RArnoldi[0, 0]
else:
r = ((r - self.RHSs[:, :ns-1].dot(self.RArnoldi[:ns-1, ns-1]))
/ self.RArnoldi[ns-1, ns-1])
if overwrite:
self.RHSs[:, ns - 1] = r
else:
if ns == 1:
self.RHSs = r.reshape((- 1, 1))
else:
self.RHSs = np.hstack((self.RHSs, r[:, None]))
return r
def postprocessu(self, u:Np1D, overwrite : bool = False):
if self.verbosity >= 10:
verbosityDepth("INIT", "Starting orthogonalization.")
ns = self.nsamples
nsAug = (ns + 1) * self.HFEngine.V.dim()
if ns == 0:
u, h, _ = self.PODEngine.GS(u, np.empty((0, 0)))
r = h[0]
uAug = copy(u)
else:
uAug = np.concatenate((self.samplesAug[self.HFEngine.V.dim()
- nsAug :, ns - 1],
u), axis = None)
u, h, uAug = self.PODEngine.GS(u, self.samples[:, : ns], ns, uAug,
self.samplesAug[- nsAug :, : ns])
if overwrite:
self.HArnoldi[: ns + 1, ns] = h
if ns > 0:
r = self.HArnoldi[: ns + 1, 1 : ns + 1].dot(
self.RArnoldi[: ns, ns - 1])
self.RArnoldi[: ns + 1, ns] = r
self.samplesAug[- nsAug :, ns] = uAug
else:
if ns == 0:
self.HArnoldi = h.reshape((1, 1))
self.RArnoldi = r.reshape((1, 1))
self.samplesAug = uAug.reshape((-1, 1))
else:
self.HArnoldi=np.block([[ self.HArnoldi, h[:-1, None]],
[np.zeros((1, ns)), h[-1]]])
if ns > 0:
r = self.HArnoldi[: ns + 1, 1 : ns + 1].dot(
self.RArnoldi[: ns, ns - 1])
self.RArnoldi=np.block([[ self.RArnoldi, r[:-1, None]],
[np.zeros((1, ns)), r[-1]]])
self.samplesAug=np.vstack((np.zeros((self.HFEngine.V.dim(),
ns)),
self.samplesAug))
self.samplesAug = np.hstack((self.samplesAug, uAug[:, None]))
if self.verbosity >= 10:
verbosityDepth("DEL", "Done orthogonalizing.")
return u
def preallocateSamples(self, u:Np1D, n:int):
super().preallocateSamples(u, n)
h = self.HArnoldi
r = self.RArnoldi
saug = self.samplesAug
self.HArnoldi = np.zeros((n, n), dtype = u.dtype)
self.HArnoldi[0, 0] = h[0, 0]
self.RArnoldi = np.zeros((n, n), dtype = u.dtype)
self.RArnoldi[0, 0] = r[0, 0]
self.RHSs = np.empty((u.size, n - 1), dtype = u.dtype)
self.samplesAug = np.zeros((self.HFEngine.V.dim() * (n + 1), n),
dtype = u.dtype)
self.samplesAug[- self.HFEngine.V.dim() :, 0] = saug[:, 0]
c4science · Help