Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90026510
ROMApproximantLagrangeVF.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, Oct 28, 14:49
Size
13 KB
Mime Type
text/x-python
Expires
Wed, Oct 30, 14:49 (2 d)
Engine
blob
Format
Raw Data
Handle
21995755
Attached To
R6746 RationalROMPy
ROMApproximantLagrangeVF.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
class
ROMApproximantLagrangeVF
(
ROMApproximantLagrange
):
"""
ROM Lagrange VF 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 mu);
- setProblemData: set HF problem data and mu 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.
mus(optional): Boundaries of snapshot positions. Defaults to
np.array([0, 1]).
ws(optional): Array of snapshot weigths. Defaults to
np.ones_like(self.mus).
w(optional): Weight for norm computation. If None, set to
Re(np.mean(self.mus)). 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 VF interpolant numerator; defaults to 0;
- 'N': degree of VF interpolant denominator; defaults to 0.
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 {}.
Attributes:
HFEngine: HF problem solver.
HSEngine: Hilbert space general purpose engine.
mu0: Default parameter.
mus: Boundaries of snapshot positions.
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 VF interpolant numerator;
- 'N': degree of VF interpolant denominator;
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
autoNodeTypes: List of possible values for autoNode.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
S: Number of solution snapshots over which current approximant is
based upon.
plotSnap: What to plot of snapshots of the Helmholtz solution map.
plotSnapSpecs: How to plot snapshots of the Helmholtz solution map.
solSnapshots: Array whose columns are FE dofs of snapshots of solution
map.
RPOD: Right factor of snapshots POD. If no POD, set to None.
POD: Whether to compute POD of snapshots.
autoNode: Type of nodes, if automatically generated. Otherwise None.
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.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix wrt mu.
bs: List of numpy vectors representing coefficients of linear system
RHS wrt mu.
"""
extraApproxParameters
=
[
"M"
,
"N"
]
def
parameterList
(
self
)
->
ListAny
:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return
(
ROMApproximantLagrange
.
parameterList
(
self
)
+
ROMApproximantLagrangeVF
.
extraApproxParameters
)
@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
,
ROMApproximantLagrangeVF
.
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
,
ROMApproximantLagrangeVF
.
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
@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."
,
stacklevel
=
2
)
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."
,
stacklevel
=
2
)
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
]),
stacklevel
=
2
)
self
.
Emax
=
vals
[
idxmax
]
+
1
else
:
self
.
_S
=
S
self
.
_approxParameters
[
"S"
]
=
self
.
S
if
Sold
!=
self
.
S
:
self
.
resetSamples
()
def
approxRadius
(
self
)
->
float
:
"""Sample radius."""
return
np
.
max
(
self
.
mu0
-
self
.
mus
)
def
setupApprox
(
self
):
"""Compute VF interpolant. SVD-based robust eigenvalue management."""
if
not
self
.
checkComputedApprox
():
S0
=
self
.
S
M1
=
self
.
M
+
1
N1
=
self
.
N
+
1
if
self
.
autoNode
not
in
self
.
autoNodeTypes
[
1
:]:
raise
Exception
(
"Manual sample selection not yet implemented."
)
self
.
computeSnapshots
()
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
.
radiusVF
(
self
.
mus
),
c
)
elif
self
.
polyBasis
==
"GAUSSLEGENDRE"
:
Phi
[:,
j
]
=
(
j
+
.
5
)
**
.
5
*
np
.
polynomial
.
legendre
.
legval
(
self
.
radiusVF
(
self
.
mus
),
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
):
#twice is enough!
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
)
phi
=
eV
[
-
1
,
:]
.
conj
()
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
radiusVF
(
self
,
mu
:
complex
)
->
float
:
"""
Compute translated radius to be plugged into VF approximant.
Args:
mu: Target parameter.
Returns:
Translated radius to be plugged into VF approximant.
"""
return
(
mu
-
self
.
mu0
)
/
self
.
approxRadius
()
def
evalApprox
(
self
,
mu
:
complex
)
->
Np1D
:
"""
Evaluate VF approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
Returns:
Real and imaginary parts of approximant.
"""
self
.
setupApprox
()
powerlist
=
np
.
power
(
self
.
radiusVF
(
mu
),
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
)
->
Np1D
:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self
.
setupApprox
()
return
np
.
roots
(
self
.
Q
[::
-
1
])
*
self
.
approxRadius
()
+
self
.
mu0
Event Timeline
Log In to Comment