Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62853150
approximant_taylor_pade.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, May 16, 02:28
Size
17 KB
Mime Type
text/x-python
Expires
Sat, May 18, 02:28 (2 d)
Engine
blob
Format
Raw Data
Handle
17698696
Attached To
R6746 RationalROMPy
approximant_taylor_pade.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
warnings
import
numpy
as
np
from
rrompy.reduction_methods.taylor.generic_approximant_taylor
import
(
GenericApproximantTaylor
)
from
rrompy.sampling.scipy.pod_engine
import
PODEngine
from
rrompy.utilities.base.types
import
Np1D
,
Np2D
,
ListAny
,
Tuple
,
DictAny
from
rrompy.utilities.base.types
import
HFEng
from
rrompy.utilities.base
import
purgeDict
__all__
=
[
'ApproximantTaylorPade'
]
class
ApproximantTaylorPade
(
GenericApproximantTaylor
):
"""
ROM single-point fast Pade' approximant computation for parametric
problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'rho': weight for computation of original Pade' approximant;
defaults to np.inf, i.e. fast approximant;
- 'M': degree of Pade' approximant numerator; defaults to 0;
- 'N': degree of Pade' approximant denominator; defaults to 0;
- 'E': total number of derivatives current approximant relies upon;
defaults to Emax;
- 'Emax': total number of derivatives of solution map to be
computed; defaults to E;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0;
- 'sampleType': label of sampling type; available values are:
- 'ARNOLDI': orthogonalization of solution derivatives through
Arnoldi algorithm;
- 'KRYLOV': standard computation of solution derivatives.
Defaults to 'KRYLOV'.
Defaults to empty dict.
plotDer(optional): What to plot of derivatives of the Helmholtz
solution map. See plot method in HFEngine. Defaults to [].
plotDerSpecs(optional): How to plot derivatives of the Helmholtz
solution map. See plot method in HFEngine. Defaults to {}.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'rho': weight for computation of original Pade' approximant;
- 'M': degree of Pade' approximant numerator;
- 'N': degree of Pade' approximant denominator;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'robustTol': tolerance for robust Pade' denominator management;
- 'sampleType': label of sampling type.
POD: Whether to compute QR factorization of derivatives.
rho: Weight of approximant.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
robustTol: tolerance for robust Pade' denominator management.
sampleType: Label of sampling type.
plotDer: What to plot of derivatives of the Helmholtz solution map.
plotDerSpecs: How to plot derivatives of the Helmholtz solution map.
initialHFData: HF problem initial data.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
G: Square Numpy 2D vector of size (N+1) corresponding to Pade'
denominator matrix (see paper).
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
def
__init__
(
self
,
HFEngine
:
HFEng
,
mu0
:
complex
=
0
,
approxParameters
:
DictAny
=
{},
plotDer
:
ListAny
=
[],
plotDerSpecs
:
DictAny
=
{}):
if
not
hasattr
(
self
,
"parameterList"
):
self
.
parameterList
=
[]
self
.
parameterList
+=
[
"M"
,
"N"
,
"robustTol"
,
"rho"
]
super
()
.
__init__
(
HFEngine
=
HFEngine
,
mu0
=
mu0
,
approxParameters
=
approxParameters
,
plotDer
=
plotDer
,
plotDerSpecs
=
plotDerSpecs
)
@property
def
approxParameters
(
self
):
"""Value of approximant parameters."""
return
self
.
_approxParameters
@approxParameters.setter
def
approxParameters
(
self
,
approxParams
):
approxParameters
=
purgeDict
(
approxParams
,
self
.
parameterList
,
dictname
=
self
.
name
()
+
".approxParameters"
,
baselevel
=
1
)
approxParametersCopy
=
purgeDict
(
approxParameters
,
[
"M"
,
"N"
,
"robustTol"
,
"rho"
],
True
,
True
,
baselevel
=
1
)
keyList
=
list
(
approxParameters
.
keys
())
if
"rho"
in
keyList
:
self
.
_rho
=
approxParameters
[
"rho"
]
elif
hasattr
(
self
,
"rho"
):
self
.
_rho
=
self
.
rho
else
:
self
.
_rho
=
np
.
inf
GenericApproximantTaylor
.
approxParameters
.
fset
(
self
,
approxParametersCopy
)
self
.
rho
=
self
.
_rho
if
"robustTol"
in
keyList
:
self
.
robustTol
=
approxParameters
[
"robustTol"
]
elif
hasattr
(
self
,
"robustTol"
):
self
.
robustTol
=
self
.
robustTol
else
:
self
.
robustTol
=
0
self
.
_ignoreParWarnings
=
True
if
"M"
in
keyList
:
self
.
M
=
approxParameters
[
"M"
]
elif
hasattr
(
self
,
"M"
):
self
.
M
=
self
.
M
else
:
self
.
M
=
0
del
self
.
_ignoreParWarnings
if
"N"
in
keyList
:
self
.
N
=
approxParameters
[
"N"
]
elif
hasattr
(
self
,
"N"
):
self
.
N
=
self
.
N
else
:
self
.
N
=
0
@property
def
rho
(
self
):
"""Value of rho."""
return
self
.
_rho
@rho.setter
def
rho
(
self
,
rho
):
self
.
_rho
=
np
.
abs
(
rho
)
self
.
_approxParameters
[
"rho"
]
=
self
.
rho
@property
def
M
(
self
):
"""Value of M. Its assignment may change Emax and E."""
return
self
.
_M
@M.setter
def
M
(
self
,
M
):
if
M
<
0
:
raise
ArithmeticError
(
"M must be non-negative."
)
self
.
_M
=
M
self
.
_approxParameters
[
"M"
]
=
self
.
M
if
not
hasattr
(
self
,
"_ignoreParWarnings"
):
self
.
checkMNEEmax
()
@property
def
N
(
self
):
"""Value of N. Its assignment may change Emax."""
return
self
.
_N
@N.setter
def
N
(
self
,
N
):
if
N
<
0
:
raise
ArithmeticError
(
"N must be non-negative."
)
self
.
_N
=
N
self
.
_approxParameters
[
"N"
]
=
self
.
N
self
.
checkMNEEmax
()
def
checkMNEEmax
(
self
):
"""Check consistency of M, N, E, and Emax."""
M
=
self
.
M
if
hasattr
(
self
,
"_M"
)
else
0
N
=
self
.
N
if
hasattr
(
self
,
"_N"
)
else
0
E
=
self
.
E
if
hasattr
(
self
,
"_E"
)
else
M
+
N
Emax
=
self
.
Emax
if
hasattr
(
self
,
"_Emax"
)
else
M
+
N
if
self
.
rho
==
np
.
inf
:
if
Emax
<
max
(
N
,
M
):
warnings
.
warn
((
"Prescribed Emax is too small. Updating Emax "
"to max(M, N)."
),
stacklevel
=
3
)
self
.
Emax
=
max
(
N
,
M
)
if
E
<
max
(
N
,
M
):
warnings
.
warn
((
"Prescribed E is too small. Updating E to "
"max(M, N)."
),
stacklevel
=
3
)
self
.
E
=
max
(
N
,
M
)
else
:
if
Emax
<
N
+
M
:
warnings
.
warn
((
"Prescribed Emax is too small. Updating "
"Emax to M + N."
),
stacklevel
=
3
)
self
.
Emax
=
self
.
N
+
M
if
E
<
N
+
M
:
warnings
.
warn
((
"Prescribed E is too small. Updating E to "
"M + N."
),
stacklevel
=
3
)
self
.
E
=
self
.
N
+
M
@property
def
robustTol
(
self
):
"""Value of tolerance for robust Pade' denominator management."""
return
self
.
_robustTol
@robustTol.setter
def
robustTol
(
self
,
robustTol
):
if
robustTol
<
0.
:
warnings
.
warn
((
"Overriding prescribed negative robustness "
"tolerance to 0."
),
stacklevel
=
2
)
robustTol
=
0.
self
.
_robustTol
=
robustTol
self
.
_approxParameters
[
"robustTol"
]
=
self
.
robustTol
@property
def
E
(
self
):
"""Value of E. Its assignment may change Emax."""
return
self
.
_E
@E.setter
def
E
(
self
,
E
):
if
E
<
0
:
raise
ArithmeticError
(
"E must be non-negative."
)
self
.
_E
=
E
self
.
checkMNEEmax
()
self
.
_approxParameters
[
"E"
]
=
self
.
E
if
hasattr
(
self
,
"Emax"
)
and
self
.
Emax
<
self
.
E
:
warnings
.
warn
((
"Prescribed Emax is too small. Updating Emax "
"to E."
),
stacklevel
=
2
)
self
.
Emax
=
self
.
E
@property
def
Emax
(
self
):
"""Value of Emax. Its assignment may reset computed derivatives."""
return
self
.
_Emax
@Emax.setter
def
Emax
(
self
,
Emax
):
if
Emax
<
0
:
raise
ArithmeticError
(
"Emax must be non-negative."
)
if
hasattr
(
self
,
"Emax"
):
EmaxOld
=
self
.
Emax
else
:
EmaxOld
=
-
1
if
hasattr
(
self
,
"_N"
):
N
=
self
.
N
else
:
N
=
0
if
hasattr
(
self
,
"_M"
):
M
=
self
.
M
else
:
M
=
0
if
hasattr
(
self
,
"_E"
):
E
=
self
.
E
else
:
E
=
0
if
self
.
rho
==
np
.
inf
:
if
max
(
N
,
M
,
E
)
>
Emax
:
warnings
.
warn
((
"Prescribed Emax is too small. Updating Emax "
"to max(N, M, E)."
),
stacklevel
=
2
)
Emax
=
max
(
N
,
M
,
E
)
else
:
if
max
(
N
+
M
,
E
)
>
Emax
:
warnings
.
warn
((
"Prescribed Emax is too small. Updating Emax "
"to max(N + M, E)."
),
stacklevel
=
2
)
Emax
=
max
(
N
+
M
,
E
)
self
.
_Emax
=
Emax
self
.
_approxParameters
[
"Emax"
]
=
Emax
if
EmaxOld
>=
self
.
Emax
and
self
.
samplingEngine
.
samples
is
not
None
:
self
.
samplingEngine
.
samples
=
self
.
samplingEngine
.
samples
[:,
:
self
.
Emax
+
1
]
if
(
self
.
sampleType
==
"ARNOLDI"
and
self
.
samplingEngine
.
HArnoldi
is
not
None
):
self
.
samplingEngine
.
HArnoldi
=
self
.
samplingEngine
.
HArnoldi
[
:
self
.
Emax
+
1
,
:
self
.
Emax
+
1
]
self
.
samplingEngine
.
RArnoldi
=
self
.
samplingEngine
.
RArnoldi
[
:
self
.
Emax
+
1
,
:
self
.
Emax
+
1
]
def
setupApprox
(
self
):
"""
Compute Pade' approximant. SVD-based robust eigenvalue management.
"""
if
not
self
.
checkComputedApprox
():
self
.
computeDerivatives
()
while
True
:
if
self
.
POD
:
ev
,
eV
=
self
.
findeveVGQR
()
else
:
ev
,
eV
=
self
.
findeveVGExplicit
()
ts
=
self
.
robustTol
*
np
.
linalg
.
norm
(
ev
)
Nnew
=
np
.
sum
(
np
.
abs
(
ev
)
>=
ts
)
diff
=
self
.
N
-
Nnew
if
diff
<=
0
:
break
Enew
=
self
.
E
-
diff
Mnew
=
min
(
self
.
M
,
Enew
)
if
Mnew
==
self
.
M
:
strM
=
""
else
:
strM
=
", M from {0} to {1},"
.
format
(
self
.
M
,
Mnew
)
print
((
"Smallest {0} eigenvalues below tolerance.
\n
"
"Reducing N from {1} to {2}{5} and E from {3} to {4}."
)
\
.
format
(
diff
+
1
,
self
.
N
,
Nnew
,
self
.
E
,
Enew
,
strM
))
newParameters
=
{
"N"
:
Nnew
,
"M"
:
Mnew
,
"E"
:
Enew
}
self
.
approxParameters
=
newParameters
self
.
Q
=
eV
[::
-
1
,
0
]
QToeplitz
=
np
.
zeros
((
self
.
Emax
+
1
,
self
.
M
+
1
),
dtype
=
np
.
complex
)
for
i
in
range
(
len
(
self
.
Q
)):
rng
=
np
.
arange
(
self
.
M
+
1
-
i
)
QToeplitz
[
rng
,
rng
+
i
]
=
self
.
Q
[
i
]
if
self
.
sampleType
==
"ARNOLDI"
:
QToeplitz
=
self
.
samplingEngine
.
RArnoldi
.
dot
(
QToeplitz
)
self
.
P
=
self
.
samplingEngine
.
samples
.
dot
(
QToeplitz
)
self
.
lastApproxParameters
=
copy
(
self
.
approxParameters
)
def
buildG
(
self
):
"""Assemble Pade' denominator matrix."""
self
.
computeDerivatives
()
if
self
.
rho
==
np
.
inf
:
Nmin
=
self
.
E
-
self
.
N
else
:
Nmin
=
self
.
M
-
self
.
N
+
1
if
self
.
sampleType
==
"KRYLOV"
:
DerE
=
self
.
samplingEngine
.
samples
[:,
Nmin
:
self
.
E
+
1
]
G
=
self
.
HFEngine
.
innerProduct
(
DerE
,
DerE
)
else
:
RArnE
=
self
.
samplingEngine
.
RArnoldi
[:
self
.
E
+
1
,
Nmin
:
self
.
E
+
1
]
G
=
RArnE
.
conj
()
.
T
.
dot
(
RArnE
)
if
self
.
rho
==
np
.
inf
:
self
.
G
=
G
else
:
Gbig
=
G
self
.
G
=
np
.
zeros
((
self
.
N
+
1
,
self
.
N
+
1
),
dtype
=
np
.
complex
)
for
k
in
range
(
self
.
E
-
self
.
M
):
self
.
G
+=
self
.
rho
**
(
2
*
k
)
*
Gbig
[
k
:
k
+
self
.
N
+
1
,
k
:
k
+
self
.
N
+
1
]
def
findeveVGExplicit
(
self
)
->
Tuple
[
Np1D
,
Np2D
]:
"""
Compute explicitly eigenvalues and eigenvectors of Pade' denominator
matrix.
"""
self
.
buildG
()
ev
,
eV
=
np
.
linalg
.
eigh
(
self
.
G
)
return
ev
,
eV
def
findeveVGQR
(
self
)
->
Tuple
[
Np1D
,
Np2D
]:
"""
Compute eigenvalues and eigenvectors of Pade' denominator matrix
through SVD of R factor. See ``Householder triangularization of a
quasimatrix'', L.Trefethen, 2008 for QR algorithm.
Returns:
Eigenvalues in ascending order and corresponding eigenvector
matrix.
"""
self
.
computeDerivatives
()
if
self
.
rho
==
np
.
inf
:
Nmin
=
self
.
E
-
self
.
N
else
:
Nmin
=
self
.
M
-
self
.
N
+
1
if
self
.
sampleType
==
"KRYLOV"
:
A
=
copy
(
self
.
samplingEngine
.
samples
[:,
Nmin
:
self
.
E
+
1
])
self
.
PODEngine
=
PODEngine
(
self
.
HFEngine
)
R
=
self
.
PODEngine
.
QRHouseholder
(
A
,
only_R
=
True
)
else
:
R
=
self
.
samplingEngine
.
RArnoldi
[:
self
.
E
+
1
,
Nmin
:
self
.
E
+
1
]
if
self
.
rho
==
np
.
inf
:
_
,
s
,
V
=
np
.
linalg
.
svd
(
R
,
full_matrices
=
False
)
else
:
Rtower
=
np
.
zeros
((
R
.
shape
[
0
]
*
(
self
.
E
-
self
.
M
),
self
.
N
+
1
),
dtype
=
np
.
complex
)
for
k
in
range
(
self
.
E
-
self
.
M
):
Rtower
[
k
*
R
.
shape
[
0
]
:
(
k
+
1
)
*
R
.
shape
[
0
],
:]
=
(
self
.
rho
**
k
*
R
[:,
self
.
M
-
self
.
N
+
1
+
k
:
self
.
M
+
1
+
k
])
_
,
s
,
V
=
np
.
linalg
.
svd
(
Rtower
,
full_matrices
=
False
)
return
s
[::
-
1
],
V
.
conj
()
.
T
[:,
::
-
1
]
def
radiusPade
(
self
,
mu
:
Np1D
,
mu0
:
float
=
None
)
->
float
:
"""
Compute translated radius to be plugged into Pade' approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Translated radius to be plugged into Pade' approximant.
"""
if
mu0
is
None
:
mu0
=
self
.
mu0
return
self
.
HFEngine
.
rescaling
(
mu
)
-
self
.
HFEngine
.
rescaling
(
mu0
)
def
evalApprox
(
self
,
mu
:
complex
):
"""
Evaluate Pade' approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
"""
self
.
setupApprox
()
powerlist
=
np
.
power
(
self
.
radiusPade
(
mu
),
range
(
max
(
self
.
M
,
self
.
N
)
+
1
))
self
.
uApp
=
(
self
.
P
.
dot
(
powerlist
[:
self
.
M
+
1
])
/
self
.
Q
.
dot
(
powerlist
[:
self
.
N
+
1
]))
def
getPoles
(
self
)
->
Np1D
:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self
.
setupApprox
()
return
self
.
HFEngine
.
rescalingInv
(
np
.
roots
(
self
.
Q
[::
-
1
])
+
self
.
HFEngine
.
rescaling
(
self
.
mu0
))
Event Timeline
Log In to Comment