Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85450678
ROMApproximantLagrangePade.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
Sun, Sep 29, 08:01
Size
19 KB
Mime Type
text/x-python
Expires
Tue, Oct 1, 08:01 (2 d)
Engine
blob
Format
Raw Data
Handle
21184030
Attached To
R6746 RationalROMPy
ROMApproximantLagrangePade.py
View Options
#!/usr/bin/python
from
copy
import
copy
import
warnings
import
numpy
as
np
import
utilities
from
ROMApproximantLagrange
import
ROMApproximantLagrange
from
RROMPyTypes
import
Np1D
,
ListAny
,
DictAny
,
HFEng
,
HSEng
PI
=
np
.
pi
class
ROMApproximantLagrangePade
(
ROMApproximantLagrange
):
"""
ROM Lagrange Pade' interpolant 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;
- getLSBlocks: get algebraic system blocks.
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.
ks(optional): Array of snapshot parameters. Defaults to np.array([0]).
ws(optional): Array of snapshots weights. Defaults to uniform = 1.
w(optional): Weight for norm computation. If None, set to
Re(np.mean(ks)). 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;
- 'S': total number of samples current approximant relies upon;
defaults to 2;
- 'M': degree of Pade' interpolant numerator; defaults to 0;
- 'N': degree of Pade' interpolant denominator; defaults to 0.
- 'polyBasis': label of polynomial basis for LS problem; available
values are:
- 'CHEBYSHEV': Chebyshev nodes and weights;
- 'GAUSSLEGENDRE': Gauss-Legendre nodes and weights;
defaults to 'CHEBYSHEV'.
Defaults to empty dict.
plotSnap(optional): What to plot of snapshots of the Helmholtz
solution map. See plot method in HSEngine. Defaults to [].
plotSnapSpecs(optional): How to plot snapshots of the Helmholtz
solution map. See plot method in HSEngine. Defaults to {}.
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.
solSnapshots: Array whose columns are FE dofs of snapshots of solution
map.
k0: Default parameter.
ks: Array of snapshot parameters.
ws: Array of snapshots weights.
w: Weight for norm computation.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of samples current approximant relies upon;
- 'M': degree of Pade' interpolant numerator;
- 'N': degree of Pade' interpolant denominator;
- 'polyBasis': label of polynomial basis for LS problem.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
S: Number of solution snapshots over which current approximant is
based upon.
polyBasis: Polynomial basis for LS problem.
plotSnap: What to plot of snapshots of the Helmholtz solution map.
plotSnapSpecs: How to plot snapshots of the Helmholtz solution map.
POD: Whether to compute POD of snapshots.
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.
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"
,
"polyBasis"
]
polyBasisParameters
=
[
"CHEBYSHEV"
,
"GAUSSLEGENDRE"
]
def
__init__
(
self
,
HFEngine
:
HFEng
,
HSEngine
:
HSEng
,
ks
:
Np1D
=
np
.
array
([
0
]),
ws
:
Np1D
=
None
,
w
:
float
=
None
,
approxParameters
:
DictAny
=
{},
plotSnap
:
ListAny
=
[],
plotSnapSpecs
:
DictAny
=
{},
verboseRobust
:
int
=
1
,
cleanupParameters
:
DictAny
=
{}):
ROMApproximantLagrange
.
__init__
(
self
,
HFEngine
,
HSEngine
,
ks
,
w
,
approxParameters
,
plotSnap
,
plotSnapSpecs
)
if
ws
==
None
:
ws
=
np
.
ones
(
np
.
size
(
self
.
ks
))
self
.
ws
=
ws
self
.
verboseRobust
=
verboseRobust
self
.
cleanupParameters
=
cleanupParameters
def
parameterList
(
self
)
->
ListAny
:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return
(
ROMApproximantLagrange
.
parameterList
(
self
)
+
ROMApproximantLagrangePade
.
extraApproxParameters
)
@property
def
k0
(
self
):
"""Dummy center of approximant (i.e. k0)."""
self
.
k0
=
np
.
mean
(
self
.
ks
)
return
self
.
_k0
@k0.setter
def
k0
(
self
,
k0
:
bool
):
self
.
_k0
=
k0
@property
def
approxParameters
(
self
):
"""
Value of approximant parameters. Its assignment may change M, N and S.
"""
return
self
.
_approxParameters
@approxParameters.setter
def
approxParameters
(
self
,
approxParams
):
if
not
hasattr
(
self
,
"approxParameters"
):
self
.
_approxParameters
=
{}
approxParameters
=
utilities
.
purgeDict
(
approxParams
,
ROMApproximantLagrangePade
.
parameterList
(
self
),
dictname
=
self
.
name
()
+
".approxParameters"
)
keyList
=
list
(
approxParameters
.
keys
())
if
"M"
in
keyList
:
self
.
M
=
0
#to avoid warnings if M and S are changed simultaneously
if
"N"
in
keyList
:
self
.
N
=
0
#to avoid warnings if N and S are changed simultaneously
approxParametersCopy
=
utilities
.
purgeDict
(
approxParameters
,
ROMApproximantLagrangePade
.
extraApproxParameters
,
True
,
True
)
ROMApproximantLagrange
.
approxParameters
.
fset
(
self
,
approxParametersCopy
)
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
"polyBasis"
in
keyList
:
self
.
polyBasis
=
approxParameters
[
"polyBasis"
]
elif
hasattr
(
self
,
"polyBasis"
):
self
.
polyBasis
=
self
.
polyBasis
else
:
self
.
polyBasis
=
"CHEBYSHEV"
@property
def
M
(
self
):
"""Value of M. Its assignment may change S."""
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
,
"S"
)
and
self
.
S
<
self
.
M
+
1
:
warnings
.
warn
(
"Prescribed S is too small. Updating S to M + 1."
)
self
.
S
=
self
.
M
+
1
@property
def
N
(
self
):
"""Value of N. Its assignment may change S."""
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
,
"S"
)
and
self
.
S
<
self
.
N
+
1
:
warnings
.
warn
(
"Prescribed S is too small. Updating S to N + 1."
)
self
.
S
=
self
.
N
+
1
@property
def
S
(
self
):
"""Value of S."""
return
self
.
_S
@S.setter
def
S
(
self
,
S
):
if
S
<=
0
:
raise
ArithmeticError
(
"S must be positive."
)
if
hasattr
(
self
,
"S"
):
Sold
=
self
.
S
else
:
Sold
=
-
1
vals
,
label
=
[
0
]
*
2
,
{
0
:
"M"
,
1
:
"N"
}
if
hasattr
(
self
,
"M"
):
vals
[
0
]
=
self
.
M
if
hasattr
(
self
,
"N"
):
vals
[
1
]
=
self
.
N
idxmax
=
np
.
argmax
(
vals
)
if
vals
[
idxmax
]
+
1
>
S
:
warnings
.
warn
(
"Prescribed S is too small. Updating S to {} + 1."
\
.
format
(
label
[
idxmax
]))
self
.
Emax
=
vals
[
idxmax
]
+
1
else
:
self
.
_S
=
S
self
.
_approxParameters
[
"S"
]
=
self
.
S
if
Sold
!=
self
.
S
:
self
.
resetSamples
()
@property
def
polyBasis
(
self
):
"""Value of polyBasis."""
return
self
.
_polyBasis
@polyBasis.setter
def
polyBasis
(
self
,
polyBasis
):
if
hasattr
(
self
,
"polyBasis"
):
polyBasisold
=
self
.
polyBasis
else
:
polyBasisold
=
-
1
try
:
polyBasis
=
polyBasis
.
upper
()
.
strip
()
.
replace
(
" "
,
""
)
if
polyBasis
not
in
self
.
polyBasisParameters
:
raise
except
:
warnings
.
warn
((
"Prescribed polyBasis not recognized. Overriding to"
" 'CHEBYSHEV'."
))
polyBasis
=
"CHEBYSHEV"
self
.
_polyBasis
=
polyBasis
self
.
_approxParameters
[
"polyBasis"
]
=
self
.
polyBasis
if
polyBasisold
!=
self
.
polyBasis
:
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
}
def
setupApprox
(
self
):
"""
Compute Pade' interpolant.
SVD-based robust eigenvalue management.
"""
S0
=
self
.
S
M1
=
self
.
M
+
1
N1
=
self
.
N
+
1
if
self
.
solSnapshots
is
None
:
self
.
approxRadius
=
np
.
max
(
np
.
abs
(
self
.
k0
-
self
.
ks
))
if
self
.
polyBasis
==
"CHEBYSHEV"
:
nodes
,
weights
=
np
.
polynomial
.
chebyshev
.
chebgauss
(
S0
)
self
.
ws
=
weights
*
2.
/
PI
elif
self
.
polyBasis
==
"GAUSSLEGENDRE"
:
nodes
,
weights
=
np
.
polynomial
.
legendre
.
leggauss
(
S0
)
self
.
ws
=
weights
self
.
ks
=
nodes
*
self
.
approxRadius
+
self
.
k0
self
.
ws
=
self
.
ws
[:,
None
]
self
.
computeSnapshots
()
if
not
self
.
checkComputedApprox
():
Phi
=
np
.
zeros
((
S0
,
max
(
M1
,
N1
)),
dtype
=
np
.
complex
)
Phi
[:,
0
]
=
np
.
ones
((
S0
,))
/
2
**.
5
for
j
in
range
(
1
,
max
(
M1
,
N1
)):
c
=
np
.
zeros
((
j
+
1
,))
c
[
-
1
]
=
1.
if
self
.
polyBasis
==
"CHEBYSHEV"
:
Phi
[:,
j
]
=
np
.
polynomial
.
chebyshev
.
chebval
(
self
.
radiusPade
(
self
.
ks
),
c
)
elif
self
.
polyBasis
==
"GAUSSLEGENDRE"
:
Phi
[:,
j
]
=
(
j
+
.
5
)
**
.
5
*
np
.
polynomial
.
legendre
.
legval
(
self
.
radiusPade
(
self
.
ks
),
c
)
if
self
.
POD
:
II
=
np
.
array
(
np
.
arange
(
0
,
S0
**
3
,
S0
)
+
np
.
kron
(
np
.
arange
(
S0
),
np
.
ones
(
S0
)),
dtype
=
np
.
int
)
Rtall
=
np
.
zeros
(
S0
**
3
,
dtype
=
np
.
complex
)
Rtall
[
II
]
=
np
.
reshape
(
self
.
RPOD
.
T
,
(
S0
**
2
,))
Rtall
=
np
.
reshape
(
Rtall
,
(
S0
,
S0
**
2
))
.
T
B
=
Rtall
.
dot
(
Phi
[:,
:
N1
])
Z
=
copy
(
B
)
Z
=
np
.
reshape
(
Z
.
T
,
(
S0
*
(
N1
),
S0
))
.
T
for
j
in
range
(
2
):
Z
=
Z
-
Phi
[:,
:
M1
]
.
dot
(
Phi
[:,
:
M1
]
.
conj
()
.
T
.
dot
(
np
.
multiply
(
self
.
ws
,
Z
)))
Z
=
np
.
reshape
(
Z
.
T
,
(
N1
,
S0
**
2
))
.
T
else
:
ker
=
self
.
solSnapshots
.
conj
()
.
T
.
dot
(
self
.
energyNormMat
.
dot
(
self
.
solSnapshots
))
WPhi
=
np
.
reshape
(
np
.
multiply
(
self
.
ws
,
Phi
[:,
:
M1
])
.
T
,
(
S0
*
M1
))
.
conj
()[:,
None
]
Y
=
np
.
multiply
(
WPhi
,
np
.
kron
(
np
.
ones
((
M1
,
1
)),
Phi
[:,
:
N1
]))
Ylarge
=
np
.
reshape
(
Y
.
T
,
(
M1
*
N1
,
S0
))
.
T
B
=
np
.
reshape
(
ker
.
dot
(
Ylarge
)
.
T
,
(
N1
,
M1
*
S0
))
.
T
D
=
np
.
multiply
(
self
.
ws
,
np
.
diag
(
ker
)[:,
None
])
D
=
Phi
[:,
:
N1
]
.
conj
()
.
T
.
dot
(
np
.
multiply
(
D
,
Phi
[:,
:
N1
]))
Z
=
B
.
conj
()
.
T
.
dot
(
Y
)
-
D
_
,
ev
,
eV
=
np
.
linalg
.
svd
(
Z
,
full_matrices
=
False
)
eV
=
eV
.
conj
()
.
T
phi
=
eV
[:,
np
.
argmin
(
ev
)]
if
self
.
POD
:
c
=
Phi
[:,
:
M1
]
.
conj
()
.
T
.
dot
(
np
.
multiply
(
self
.
ws
,
np
.
reshape
(
B
.
dot
(
phi
)
.
T
,
(
S0
,
S0
))
.
T
))
else
:
c
=
np
.
reshape
(
Y
.
dot
(
phi
)
.
T
,
(
M1
,
S0
))
polybase
=
np
.
zeros
((
max
(
M1
,
N1
),
max
(
M1
,
N1
)))
polybase
[
0
,
0
]
=
1
polybase
[
1
,
1
]
=
1
for
j
in
range
(
2
,
max
(
M1
,
N1
)):
if
self
.
polyBasis
==
"CHEBYSHEV"
:
polybase
[
1
:
j
+
1
,
j
]
=
2
*
polybase
[:
j
,
j
-
1
]
polybase
[:
j
+
1
,
j
]
=
(
polybase
[:
j
+
1
,
j
]
-
polybase
[:
j
+
1
,
j
-
2
])
elif
self
.
polyBasis
==
"GAUSSLEGENDRE"
:
polybase
[
1
:
j
+
1
,
j
]
=
((
2
*
j
-
1
)
/
j
*
polybase
[:
j
,
j
-
1
])
polybase
[:
j
+
1
,
j
]
=
(
polybase
[:
j
+
1
,
j
]
-
((
j
-
1
)
/
j
*
polybase
[:
j
+
1
,
j
-
2
]))
if
self
.
polyBasis
==
"CHEBYSHEV"
:
polybase
[
0
,
0
]
=
.
5
**
.
5
elif
self
.
polyBasis
==
"GAUSSLEGENDRE"
:
polybase
=
np
.
multiply
(
np
.
power
(
np
.
arange
(
.
5
,
max
(
M1
,
N1
)),
.
5
),
polybase
)
self
.
P
=
self
.
solSnapshots
.
dot
(
polybase
[:
M1
,
:
M1
]
.
dot
(
c
)
.
T
)
self
.
Q
=
polybase
[:
N1
,
:
N1
]
.
dot
(
phi
)
self
.
approxParameters
=
{
"N"
:
self
.
N
,
"M"
:
self
.
M
,
"S"
:
S0
}
self
.
lastApproxParameters
=
copy
(
self
.
approxParameters
)
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
radiusPade
(
self
,
k
:
complex
):
"""
Compute translated radius to be plugged into Pade' approximant.
Args:
k: Target wavenumber.
Returns:
Translated radius to be plugged into Pade' approximant.
"""
return
(
k
-
self
.
k0
)
/
self
.
approxRadius
def
evalApprox
(
self
,
k
:
complex
)
->
Np1D
:
"""
Evaluate Pade' approximant at arbitrary wavenumber.
Args:
k: Target wavenumber.
Returns:
Real and imaginary parts of approximant.
"""
self
.
setupApprox
()
powerlist
=
np
.
power
(
self
.
radiusPade
(
k
),
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
.
approxRadius
+
self
.
k0
*
centered
Event Timeline
Log In to Comment