Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78332758
ROMApproximantTaylorPade.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
Mon, Aug 19, 21:24
Size
26 KB
Mime Type
text/x-python
Expires
Wed, Aug 21, 21:24 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20027531
Attached To
R6746 RationalROMPy
ROMApproximantTaylorPade.py
View Options
#!/usr/bin/python
from
copy
import
copy
import
warnings
import
numpy
as
np
import
utilities
from
ROMApproximantTaylor
import
ROMApproximantTaylor
from
RROMPyTypes
import
Np1D
,
Np2D
,
ListAny
,
DictAny
,
Tuple
,
HFEng
,
HSEng
PI
=
np
.
pi
class
ROMApproximantTaylorPade
(
ROMApproximantTaylor
):
"""
ROM single-point fast Pade' approximant computation for parametric
problems.
Args:
HFEngine: HF problem solver. Should have members:
- Vdim: domain dimension;
- energyNormMatrix: assemble sparse matrix (in CSC format)
associated to weighted H10 norm;
- problemData: list of HF problem data (excluding k);
- setProblemData: set HF problem data and k to prescribed values;
- setupDerivativeComputation: setup HF problem data to compute j-th
solution derivative at k0;
- solve: get HF solution as complex numpy vector.
HSEngine: Hilbert space general purpose engine. Should have members:
- norm: compute Hilbert norm of expression represented as complex
numpy vector;
- plot: plot Hilbert expression represented as complex numpy vector.
k0(optional): Default parameter. Defaults to 0.
w(optional): Weight for norm computation. If None, set to Re(k0).
Defaults to None.
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 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 HSEngine. Defaults to [].
plotDerSpecs(optional): How to plot derivatives of the Helmholtz
solution map. See plot method in HSEngine. Defaults to {}.
equilibration(optional): Whether to apply equilibration to Gram matrix.
Defaults to False.
verboseRobust(optional): Verbosity level for robustness-related
parameter modifications. Defaults to 1.
cleanupParameters(optional): Parameter values for pole cleanup.
Available fields are:
- 'boolCondition'(bool function handle): if evaluation on pole
returns False, pole is removed; defaults to always True;
- 'residueCheck'(bool): whether to apply residue check; defaults to
False;
- 'residueNPoints'(int): number of sample points to be used in
residue estimation; defaults to 2;
- 'residueRadius'(float): sample radius to be used in residue
estimation by Cauchy formula; defaults to 1e-5;
- 'residueTol'(float): target tolerance to be used in residue
estimation; defaults to 1e-4.
Attributes:
HFEngine: HF problem solver.
HSEngine: Hilbert space general purpose engine.
solDerivatives: Array whose columns are FE dofs of solution
derivatives.
k0: Default parameter.
w: Weight for norm computation.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots;
- 'rho': weight for computation of original Pade' approximant;
defaults to inf, i.e. fast 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.
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.
equilibration: Whether to apply equilibration to Gram matrix.
verboseRobust: Verbosity level for robustness-related parameter
modifications.
cleanupParameters; Parameter values for pole cleanup. Available fields
are:
- 'boolCondition': if evaluation on pole returns False, pole is
removed;
- 'residueCheck': whether to apply residue check;
- 'residueNPoints': number of sample points to be used in residue
estimation;
- 'residueRadius': sample radius to be used in residue estimation
by Cauchy formula;
- 'residueTol': target tolerance to be used in residue estimation.
energyNormMatrix: Sparse matrix (in CSC format) assoociated to
w-weighted H10 norm.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
G: Square Numpy 2D vector of size (N+1) corresponding to Pade'
denominator matrix (see paper).
equilPowers: Numpy 1D (column) array containing equilibration powers.
Does not exist if equilibration is disabled.
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.
"""
extraApproxParameters
=
[
"M"
,
"N"
,
"robustTol"
,
"rho"
]
def
__init__
(
self
,
HFEngine
:
HFEng
,
HSEngine
:
HSEng
,
k0
:
complex
=
0
,
w
:
float
=
None
,
approxParameters
:
DictAny
=
{},
plotDer
:
ListAny
=
[],
plotDerSpecs
:
DictAny
=
{},
equilibration
:
bool
=
False
,
verboseRobust
:
int
=
1
,
cleanupParameters
:
DictAny
=
{}):
self
.
equilibration
=
equilibration
self
.
verboseRobust
=
verboseRobust
self
.
cleanupParameters
=
cleanupParameters
ROMApproximantTaylor
.
__init__
(
self
,
HFEngine
,
HSEngine
,
k0
,
w
,
approxParameters
,
plotDer
,
plotDerSpecs
)
def
parameterList
(
self
)
->
ListAny
:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return
(
ROMApproximantTaylor
.
parameterList
(
self
)
+
ROMApproximantTaylorPade
.
extraApproxParameters
)
@property
def
approxParameters
(
self
):
"""
Value of approximant parameters. Its assignment may change M, N, E,
Emax and robustTol.
"""
return
self
.
_approxParameters
@approxParameters.setter
def
approxParameters
(
self
,
approxParams
):
if
not
hasattr
(
self
,
"approxParameters"
):
self
.
_approxParameters
=
{}
approxParameters
=
utilities
.
purgeDict
(
approxParams
,
ROMApproximantTaylorPade
.
parameterList
(
self
),
dictname
=
self
.
name
()
+
".approxParameters"
)
keyList
=
list
(
approxParameters
.
keys
())
if
"M"
in
keyList
:
self
.
M
=
0
#to avoid warnings if M and E are changed simultaneously
if
"N"
in
keyList
:
self
.
N
=
0
#to avoid warnings if N and E are changed simultaneously
approxParametersCopy
=
utilities
.
purgeDict
(
approxParameters
,
ROMApproximantTaylorPade
.
extraApproxParameters
,
True
,
True
)
ROMApproximantTaylor
.
approxParameters
.
fset
(
self
,
approxParametersCopy
)
if
"rho"
in
keyList
:
self
.
rho
=
approxParameters
[
"rho"
]
elif
hasattr
(
self
,
"rho"
):
self
.
rho
=
self
.
rho
else
:
self
.
rho
=
np
.
inf
if
"M"
in
keyList
:
self
.
M
=
approxParameters
[
"M"
]
elif
hasattr
(
self
,
"M"
):
self
.
M
=
self
.
M
else
:
self
.
M
=
0
if
"N"
in
keyList
:
self
.
N
=
approxParameters
[
"N"
]
elif
hasattr
(
self
,
"N"
):
self
.
N
=
self
.
N
else
:
self
.
N
=
0
if
"robustTol"
in
keyList
:
self
.
robustTol
=
approxParameters
[
"robustTol"
]
elif
hasattr
(
self
,
"robustTol"
):
self
.
robustTol
=
self
.
robustTol
else
:
self
.
robustTol
=
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
hasattr
(
self
,
"Emax"
)
and
self
.
Emax
<
self
.
M
:
warnings
.
warn
(
"Prescribed Emax is too small. Updating Emax to M."
)
self
.
Emax
=
self
.
M
if
hasattr
(
self
,
"E"
)
and
self
.
E
<
self
.
M
:
warnings
.
warn
(
"Prescribed E is too small. Updating E to M."
)
self
.
E
=
self
.
M
@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
if
hasattr
(
self
,
"Emax"
)
and
self
.
Emax
<
self
.
N
:
warnings
.
warn
(
"Prescribed Emax is too small. Updating Emax to N."
)
self
.
Emax
=
self
.
N
if
hasattr
(
self
,
"E"
)
and
self
.
E
<
self
.
N
:
warnings
.
warn
(
"Prescribed E is too small. Updating E to N."
)
self
.
E
=
self
.
N
@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."
))
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
if
hasattr
(
self
,
"M"
)
and
self
.
M
>
self
.
E
:
warnings
.
warn
(
"Prescribed E is too small. Updating E to M."
)
self
.
E
=
self
.
M
if
hasattr
(
self
,
"N"
)
and
self
.
N
>
self
.
E
:
warnings
.
warn
(
"Prescribed E is too small. Updating E to N."
)
self
.
E
=
self
.
N
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."
)
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
vals
,
label
=
[
0
]
*
3
,
{
0
:
"E"
,
1
:
"M"
,
2
:
"N"
}
if
hasattr
(
self
,
"E"
):
vals
[
0
]
=
self
.
E
if
hasattr
(
self
,
"M"
):
vals
[
1
]
=
self
.
M
if
hasattr
(
self
,
"N"
):
vals
[
2
]
=
self
.
N
idxmax
=
np
.
argmax
(
vals
)
if
vals
[
idxmax
]
>
Emax
:
warnings
.
warn
(
"Prescribed Emax is too small. Updating Emax to {}."
\
.
format
(
label
[
idxmax
]))
self
.
Emax
=
vals
[
idxmax
]
else
:
self
.
_Emax
=
Emax
self
.
_approxParameters
[
"Emax"
]
=
Emax
if
(
EmaxOld
>
self
.
Emax
and
self
.
solDerivatives
is
not
None
and
hasattr
(
self
,
'HArnoldi'
)
and
hasattr
(
self
,
'RArnoldi'
)
and
self
.
HArnoldi
is
not
None
and
self
.
RArnoldi
is
not
None
):
self
.
solDerivatives
=
self
.
solDerivatives
[:,
:
self
.
Emax
+
1
]
if
self
.
sampleType
==
"ARNOLDI"
:
self
.
HArnoldi
=
self
.
HArnoldi
[:
self
.
Emax
+
1
,
:
self
.
Emax
+
1
]
self
.
RArnoldi
=
self
.
RArnoldi
[:
self
.
Emax
+
1
,
:
self
.
Emax
+
1
]
else
:
self
.
resetSamples
()
@property
def
cleanupParameters
(
self
):
"""Value of cleanupParameters."""
return
self
.
_cleanupParameters
@cleanupParameters.setter
def
cleanupParameters
(
self
,
cleanupParameters
):
allowedCleanupKeys
=
[
'boolCondition'
,
'residueCheck'
,
'residueNPoints'
,
'residueRadius'
,
'residueTol'
]
cleanupKeys
=
cleanupParameters
.
keys
()
if
'boolCondition'
not
in
cleanupKeys
:
cleanupParameters
[
'boolCondition'
]
=
lambda
x
:
True
if
'residueCheck'
not
in
cleanupKeys
:
cleanupParameters
[
'residueCheck'
]
=
False
if
'residueNPoints'
not
in
cleanupKeys
:
cleanupParameters
[
'residueNPoints'
]
=
2
if
'residueRadius'
not
in
cleanupKeys
:
cleanupParameters
[
'residueRadius'
]
=
1e-5
if
'residueTol'
not
in
cleanupKeys
:
cleanupParameters
[
'residueTol'
]
=
1e-4
cleanupParameters
[
'boolCondition'
]
=
np
.
vectorize
(
cleanupParameters
[
'boolCondition'
])
self
.
_cleanupParameters
=
{
key
:
cleanupParameters
[
key
]
for
key
in
allowedCleanupKeys
}
#???????????????????????????????????????????????????????????????????????????????
#from scipy.special import comb
#import scipy.sparse.linalg as spla
# def computeDerivatives(self):
# """Compute derivatives of solution map starting from order 0."""
# if self.solDerivatives is None:
# lAs = len(self.As)
# lbs = len(self.bs)
# uOlds = [np.zeros(self.HFEngine.Vdim())] * (lAs - 1)
# self.HFEngine.setProblemData(self.initialHFData, self.k0)
# bs = self.bs
# _, b = self.buildLS(self.k0, bs = bs)
# del uOlds[0]
# uOlds = uOlds + [self.manageDerivatives(b, 0)]
# for j in range(self.Emax + 1):
# bs = [np.zeros_like(uOlds[0])] * max(lAs - 1, lbs - j)
# for ii in range(lbs - j - 1):
# bs[ii] = self.bs[j + ii]
# for ii in range(lAs - 1):
# for jj in range(1, min(j + 2, lAs - ii)):
# bs[ii] =(bs[ii]
# - comb(jj + ii, jj, exact = False)
# * self.As[jj + ii].dot(uOlds[lAs - jj - 1]))
# A, b = self.buildLS(self.k0, bs = bs)
# uj = spla.spsolve(A, b)
# self.HSEngine.plot(uj, name = "u_{0}".format(j),
# what = self.plotDer, **self.plotDerSpecs)
# del uOlds[0]
# uOlds = uOlds + [self.manageDerivatives(uj, j + 1)]
#???????????????????????????????????????????????????????????????????????????????
def
setupApprox
(
self
):
"""
Compute Pade' approximant.
See Section 4 in Fast Pade' paper.
SVD-based robust eigenvalue management.
"""
self
.
computeDerivatives
()
if
not
self
.
checkComputedApprox
():
while
True
:
if
self
.
POD
:
if
self
.
sampleType
==
"ARNOLDI"
:
ev
,
eV
=
self
.
findeveVGArnoldi
()
else
:
raise
Exception
((
"To be implemented correctly. Change 'POD' "
"to False or 'sampleType' to 'ARNOLDI' if "
"possible."
))
ev
,
eV
=
self
.
findeveVGQR
()
else
:
self
.
buildG
()
ev
,
eV
=
np
.
linalg
.
eigh
(
self
.
G
)
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
)
if
self
.
verboseRobust
>=
1
:
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
]
if
self
.
equilibration
:
self
.
Q
=
np
.
multiply
(
self
.
equilPowers
.
flatten
()[::
-
1
],
self
.
Q
)
self
.
_cleanup
()
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
.
RArnoldi
.
dot
(
QToeplitz
)
self
.
P
=
self
.
solDerivatives
.
dot
(
QToeplitz
)
self
.
lastApproxParameters
=
copy
(
self
.
approxParameters
)
def
buildG
(
self
):
"""Assemble Pade' denominator matrix."""
if
self
.
N
==
0
:
self
.
G
=
np
.
array
([[
1
]],
dtype
=
np
.
complex
)
return
self
.
computeDerivatives
()
if
np
.
isinf
(
self
.
rho
):
lgSize
=
self
.
N
+
1
else
:
lgSize
=
self
.
E
+
self
.
N
-
self
.
M
Gshift
=
self
.
E
-
lgSize
+
1
largeG
=
np
.
zeros
((
lgSize
,
lgSize
),
dtype
=
np
.
complex
)
largeG
=
self
.
solDerivatives
[:,
Gshift
:
self
.
E
+
1
]
.
conj
()
.
T
.
dot
(
self
.
energyNormMatrix
.
dot
(
self
.
solDerivatives
[:,
Gshift
:
self
.
E
+
1
]))
if
np
.
isinf
(
self
.
rho
):
self
.
G
=
largeG
else
:
self
.
G
=
np
.
zeros
((
self
.
N
+
1
,
self
.
N
+
1
),
dtype
=
np
.
complex
)
for
k
in
range
(
self
.
E
-
self
.
M
):
self
.
G
=
(
self
.
G
+
self
.
rho
**
(
2
*
k
)
*
largeG
[
k
:
k
+
self
.
N
+
1
,
k
:
k
+
self
.
N
+
1
])
if
self
.
equilibration
:
Gd
=
np
.
diag
(
self
.
G
)
gamma
=
np
.
average
(
np
.
abs
(
np
.
divide
(
Gd
[
1
:],
Gd
[:
-
1
])),
weights
=
np
.
power
(
1.2
,
np
.
arange
(
self
.
N
)))
**.
5
self
.
equilPowers
=
np
.
power
(
gamma
,
np
.
arange
(
self
.
N
+
1
)[:,
None
])
self
.
G
=
np
.
multiply
(
self
.
equilPowers
,
np
.
multiply
(
self
.
equilPowers
,
self
.
G
)
.
T
)
.
T
def
_cleanup
(
self
):
"""Cleanup Pade' denominator by removing unwanted poles."""
if
self
.
N
==
0
:
return
poles
=
np
.
roots
(
self
.
Q
[::
-
1
])
+
self
.
k0
NpolesOld
=
len
(
poles
)
poles
=
poles
[
self
.
cleanupParameters
[
'boolCondition'
](
poles
)]
if
self
.
cleanupParameters
[
'residueCheck'
]:
resR
=
self
.
cleanupParameters
[
'residueRadius'
]
resTol
=
self
.
cleanupParameters
[
'residueTol'
]
residues
=
np
.
zeros_like
(
poles
)
for
l
,
pole
in
enumerate
(
poles
):
resV
=
np
.
zeros
(
self
.
solDerivatives
.
shape
[
0
],
dtype
=
np
.
complex
)
NPoints
=
self
.
cleanupParameters
[
'residueNPoints'
]
for
theta
in
2
*
PI
*
np
.
arange
(
NPoints
)
/
NPoints
:
deltapole
=
resR
*
np
.
exp
(
1.j
*
theta
)
ksample
=
pole
+
deltapole
self
.
solveHF
(
ksample
)
resV
=
deltapole
/
NPoints
*
(
resV
+
self
.
uHF
)
residues
[
l
]
=
np
.
abs
(
resV
.
dot
(
self
.
energyNormMatrix
.
dot
(
resV
)
.
conj
()))
**
.
5
poles
=
poles
[
residues
>=
resTol
]
NpolesNew
=
len
(
poles
)
if
NpolesOld
>
NpolesNew
:
if
self
.
verboseRobust
>=
1
:
sSing
=
"s"
*
(
NpolesOld
-
NpolesNew
>
1
)
print
((
"Identified {0} pole{1} to be removed.
\n
"
"Reducing N from {2} to {3}."
)
.
format
(
NpolesOld
-
NpolesNew
,
sSing
,
NpolesOld
,
NpolesNew
))
self
.
N
=
NpolesNew
newQ
=
np
.
polyfit
(
np
.
append
(
poles
-
self
.
k0
,
0.
),
np
.
append
(
np
.
zeros
(
NpolesNew
),
1.
),
NpolesNew
)
self
.
Q
=
newQ
[::
-
1
]
/
np
.
linalg
.
norm
(
newQ
)
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
()
A
=
copy
(
self
.
solDerivatives
[:,
(
self
.
E
-
self
.
N
)
:
(
self
.
E
+
1
)])
M
=
self
.
energyNormMatrix
E
=
np
.
zeros_like
(
A
,
dtype
=
np
.
complex
)
R
=
np
.
zeros
((
self
.
N
+
1
,
self
.
N
+
1
),
dtype
=
np
.
complex
)
for
k
in
range
(
self
.
N
+
1
):
E
[:,
k
]
=
(
np
.
random
.
randn
(
E
.
shape
[
0
])
+
1.j
*
np
.
random
.
randn
(
E
.
shape
[
0
]))
if
k
>
0
:
for
times
in
range
(
2
):
E
[:,
k
]
=
E
[:,
k
]
-
E
[:,
:
k
]
.
dot
(
(
E
[:,
k
]
.
conj
()
.
dot
(
M
.
dot
(
E
[:,
:
k
])))
.
conj
())
E
[:,
k
]
=
E
[:,
k
]
/
(
E
[:,
k
]
.
conj
()
.
dot
(
M
.
dot
(
E
[:,
k
])))
**
.
5
R
[
k
,
k
]
=
np
.
abs
(
A
[:,
k
]
.
conj
()
.
dot
(
M
.
dot
(
A
[:,
k
])))
**
.
5
alpha
=
E
[:,
k
]
.
conj
()
.
dot
(
M
.
dot
(
A
[:,
k
]))
if
np
.
isclose
(
np
.
abs
(
alpha
),
0.
):
s
=
1.
else
:
s
=
-
alpha
/
np
.
abs
(
alpha
)
E
[:,
k
]
=
s
*
E
[:,
k
]
v
=
R
[
k
,
k
]
*
E
[:,
k
]
-
A
[:,
k
]
for
times
in
range
(
2
):
v
=
v
-
E
[:,
:
k
]
.
dot
((
v
.
conj
()
.
dot
(
M
.
dot
(
E
[:,
:
k
])))
.
conj
())
sigma
=
np
.
abs
(
v
.
conj
()
.
dot
(
M
.
dot
(
v
)))
**
.
5
if
np
.
isclose
(
sigma
,
0.
):
v
=
E
[:,
k
]
else
:
v
=
v
/
sigma
J
=
np
.
arange
(
k
+
1
,
self
.
N
+
1
)
vtQ
=
v
.
conj
()
.
dot
(
M
.
dot
(
A
[:,
J
]))
A
[:,
J
]
=
A
[:,
J
]
-
2
*
np
.
outer
(
v
,
vtQ
)
R
[
k
,
J
]
=
E
[:,
k
]
.
conj
()
.
dot
(
M
.
dot
(
A
[:,
J
]))
A
[:,
J
]
=
A
[:,
J
]
-
np
.
outer
(
E
[:,
k
],
R
[
k
,
J
])
if
self
.
equilibration
:
Rd
=
np
.
diag
(
R
)
gamma
=
np
.
average
(
np
.
abs
(
np
.
divide
(
Rd
[
1
:],
Rd
[:
-
1
])),
weights
=
np
.
power
(
1.2
,
np
.
arange
(
self
.
N
)))
**.
5
self
.
equilPowers
=
np
.
power
(
gamma
,
np
.
arange
(
self
.
N
+
1
)[:,
None
])
R
=
np
.
multiply
(
self
.
equilPowers
,
R
.
T
)
.
T
_
,
s
,
V
=
np
.
linalg
.
svd
(
R
,
full_matrices
=
False
)
return
s
[::
-
1
],
V
.
conj
()
.
T
[:,
::
-
1
]
def
findeveVGArnoldi
(
self
)
->
Tuple
[
Np1D
,
Np2D
]:
"""
Compute eigenvalues and eigenvectors of Pade' denominator matrix
through SVD of R factor of solution derivatives orthogonalized by
Arnoldi algorithm.
Returns:
Eigenvalues in ascending order and corresponding eigenvector
matrix.
"""
self
.
computeDerivatives
()
if
self
.
sampleType
!=
"ARNOLDI"
:
raise
Exception
((
"Eigensolver different from 'ARNOLDI'."
"Arnoldi eigenproblem solver cannot be called."
))
R
=
self
.
RArnoldi
[:
self
.
E
+
1
,
self
.
E
-
self
.
N
:
self
.
E
+
1
]
if
self
.
equilibration
:
Rd
=
np
.
diag
(
R
)
gamma
=
np
.
average
(
np
.
abs
(
np
.
divide
(
Rd
[
1
:],
Rd
[:
-
1
])),
weights
=
np
.
power
(
1.2
,
np
.
arange
(
self
.
N
)))
**.
5
self
.
equilPowers
=
np
.
power
(
gamma
,
np
.
arange
(
self
.
N
+
1
)[:,
None
])
R
=
np
.
multiply
(
self
.
equilPowers
,
R
.
T
)
.
T
_
,
s
,
V
=
np
.
linalg
.
svd
(
R
,
full_matrices
=
False
)
return
s
[::
-
1
],
V
.
conj
()
.
T
[:,
::
-
1
]
def
evalApprox
(
self
,
k
:
complex
)
->
Np1D
:
"""
Evaluate Pade' approximant at arbitrary wavenumber.
Args:
k: Target wavenumber.
Returns:
Approximant as numpy complex vector.
"""
self
.
setupApprox
()
powerlist
=
np
.
power
(
k
-
self
.
k0
,
range
(
max
(
self
.
M
,
self
.
N
)
+
1
))
self
.
uApp
=
(
self
.
P
.
dot
(
powerlist
[:
self
.
M
+
1
])
/
self
.
Q
.
dot
(
powerlist
[:
self
.
N
+
1
]))
return
self
.
uApp
def
getPoles
(
self
,
centered
:
bool
=
False
)
->
Np1D
:
"""
Obtain approximant poles.
Args:
centered(optional): Whether to return pole positions relative to
approximation center. Defaults to False.
Returns:
Numpy complex vector of poles.
"""
self
.
setupApprox
()
return
np
.
roots
(
self
.
Q
[::
-
1
])
+
self
.
k0
*
centered
Event Timeline
Log In to Comment