Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F56454666
approximant_lagrange_pade_orthogonal.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, Mar 28, 09:20
Size
13 KB
Mime Type
text/x-python
Expires
Sat, Mar 30, 09:20 (2 d)
Engine
blob
Format
Raw Data
Handle
16718958
Attached To
R6746 RationalROMPy
approximant_lagrange_pade_orthogonal.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.reduction_methods.base
import
checkRobustTolerance
from
.approximant_lagrange_pade
import
ApproximantLagrangePade
from
rrompy.utilities.base.types
import
Np1D
,
List
from
rrompy.utilities.base
import
verbosityDepth
from
rrompy.utilities.warning_manager
import
warn
__all__
=
[
'ApproximantLagrangePadeOrthogonal'
]
class
ApproximantLagrangePadeOrthogonal
(
ApproximantLagrangePade
):
"""
ROM Lagrange Pade' interpolant computation for parametric problems with
rational-type interpolation basis functions.
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;
- 'S': total number of samples current approximant relies upon;
defaults to 2;
- 'sampler': sample point generator; defaults to uniform sampler on
[0, 1];
- 'E': coefficient of interpolant to be minimized; defaults to
min(S, M + 1);
- 'M': degree of Pade' interpolant numerator; defaults to 0;
- 'N': degree of Pade' interpolant denominator; defaults to 0;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
defaults to None;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths.
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;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'E': coefficient of interpolant to be minimized;
- 'M': degree of Pade' interpolant numerator;
- 'N': degree of Pade' interpolant denominator;
- 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust Pade' denominator management.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
POD: Whether to compute POD of snapshots.
interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust Pade' denominator management.
samplingEngine: Sampling engine.
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.
"""
@property
def
sampler
(
self
):
"""Value of sampler."""
return
self
.
_sampler
@sampler.setter
def
sampler
(
self
,
sampler
):
if
'generatePoints'
not
in
dir
(
sampler
):
raise
Exception
(
"Sampler type not recognized."
)
if
hasattr
(
self
,
'_sampler'
):
samplerOld
=
self
.
sampler
self
.
_sampler
=
sampler
self
.
_approxParameters
[
"sampler"
]
=
self
.
sampler
if
not
'samplerOld'
in
locals
()
or
samplerOld
!=
self
.
sampler
:
if
hasattr
(
self
.
sampler
,
"kind"
):
if
self
.
sampler
.
kind
==
"CHEBYSHEV"
:
self
.
_val
=
np
.
polynomial
.
chebyshev
.
chebval
self
.
_vander
=
np
.
polynomial
.
chebyshev
.
chebvander
self
.
_fit
=
np
.
polynomial
.
chebyshev
.
chebfit
self
.
_roots
=
np
.
polynomial
.
chebyshev
.
chebroots
self
.
_domcoeff
=
lambda
n
:
2.
**
max
(
0
,
n
-
1
)
else
:
#if self.sampler.kind in ["UNIFORM", "GAUSSLEGENDRE"]:
self
.
_val
=
np
.
polynomial
.
legendre
.
legval
self
.
_vander
=
np
.
polynomial
.
legendre
.
legvander
self
.
_fit
=
np
.
polynomial
.
legendre
.
legfit
self
.
_roots
=
np
.
polynomial
.
legendre
.
legroots
from
scipy.special
import
binom
self
.
_domcoeff
=
lambda
n
:
binom
(
2
*
n
,
n
)
else
:
self
.
_val
=
np
.
polynomial
.
polynomial
.
polyval
self
.
_vander
=
np
.
polynomial
.
polynomial
.
polyvander
self
.
_fit
=
np
.
polynomial
.
polynomial
.
polyfit
self
.
_roots
=
np
.
polynomial
.
polynomial
.
polyroots
self
.
_domcoeff
=
lambda
n
:
1.
self
.
resetSamples
()
def
setupApprox
(
self
):
"""
Compute Pade' interpolant.
SVD-based robust eigenvalue management.
"""
if
not
self
.
checkComputedApprox
():
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"INIT"
,
"Setting up {}."
.
format
(
self
.
name
()))
self
.
computeRescaleParameter
()
self
.
computeSnapshots
()
if
self
.
N
>
0
:
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
(
"Starting computation of "
"denominator."
))
TN
=
self
.
_vander
(
self
.
radiusPade
(
self
.
mus
),
self
.
N
)
while
self
.
N
>
0
:
TN
=
TN
[:,
:
self
.
N
+
1
]
if
self
.
POD
:
data
=
self
.
samplingEngine
.
RPOD
else
:
data
=
self
.
samplingEngine
.
samples
RHSFull
=
np
.
empty
((
self
.
S
,
data
.
shape
[
0
]
*
(
self
.
N
+
1
)),
dtype
=
np
.
complex
)
for
j
in
range
(
self
.
S
):
RHSFull
[
j
,
:]
=
np
.
kron
(
data
[:,
j
],
TN
[
j
,
:])
fitOut
=
self
.
_fit
(
self
.
radiusPade
(
self
.
mus
),
RHSFull
,
self
.
E
,
w
=
self
.
ws
,
full
=
True
,
rcond
=
self
.
interpRcond
)
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"MAIN"
,
(
"Fitting {} samples with "
"degree {} through {}... "
"Conditioning of LS system: "
"{:.4e}."
)
.
format
(
self
.
S
,
self
.
E
,
self
.
_fit
.
__name__
,
fitOut
[
1
][
2
][
0
]
/
fitOut
[
1
][
2
][
-
1
]))
if
fitOut
[
1
][
1
]
<
self
.
E
+
1
:
Enew
=
fitOut
[
1
][
1
]
-
1
Nnew
=
min
(
self
.
N
,
Enew
)
Mnew
=
min
(
self
.
M
,
Enew
)
if
Nnew
==
self
.
N
:
strN
=
""
else
:
strN
=
"N from {} to {} and "
.
format
(
self
.
N
,
Nnew
)
if
Mnew
==
self
.
M
:
strM
=
""
else
:
strM
=
"M from {} to {} and "
.
format
(
self
.
M
,
Mnew
)
warn
((
"Polyfit is poorly conditioned.
\n
Reducing {}{}E "
"from {} to {}."
)
.
format
(
strN
,
strM
,
self
.
E
,
Enew
))
newParameters
=
{
"N"
:
Nnew
,
"M"
:
Mnew
,
"E"
:
Enew
}
self
.
approxParameters
=
newParameters
continue
G
=
fitOut
[
0
][
-
1
,
:]
.
reshape
((
self
.
N
+
1
,
data
.
shape
[
0
]))
if
self
.
POD
:
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
(
"Solving svd for square "
"root of gramian matrix."
),
end
=
""
)
_
,
ev
,
eV
=
np
.
linalg
.
svd
(
G
,
full_matrices
=
False
)
ev
=
ev
[::
-
1
]
eV
=
eV
[::
-
1
,
:]
.
conj
()
.
T
else
:
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"INIT"
,
"Building gramian matrix."
,
end
=
""
)
G2
=
self
.
HFEngine
.
innerProduct
(
G
,
G
)
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"DEL"
,
"Done building gramian."
,
inline
=
True
)
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
(
"Solving eigenvalue "
"problem for gramian "
"matrix."
),
end
=
""
)
ev
,
eV
=
np
.
linalg
.
eigh
(
G2
)
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"DEL"
,
" Done."
,
inline
=
True
)
newParameters
=
checkRobustTolerance
(
ev
,
self
.
E
,
self
.
robustTol
)
if
not
newParameters
:
break
self
.
approxParameters
=
newParameters
if
self
.
N
<=
0
:
eV
=
np
.
ones
((
1
,
1
))
self
.
Q
=
eV
[:,
0
]
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"DEL"
,
"Done computing denominator."
)
else
:
self
.
Q
=
np
.
ones
(
1
,
dtype
=
np
.
complex
)
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
"Starting computation of numerator."
)
self
.
lastApproxParameters
=
copy
(
self
.
approxParameters
)
Qevaldiag
=
np
.
diag
(
self
.
getQVal
(
self
.
mus
))
while
self
.
M
>=
0
:
fitOut
=
self
.
_fit
(
self
.
radiusPade
(
self
.
mus
),
Qevaldiag
,
self
.
M
,
w
=
self
.
ws
,
full
=
True
,
rcond
=
self
.
interpRcond
)
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"MAIN"
,
(
"Fitting {} samples with degree "
"{} through {}... Conditioning of "
"LS system: {:.4e}."
)
.
format
(
self
.
S
,
self
.
M
,
self
.
_fit
.
__name__
,
fitOut
[
1
][
2
][
0
]
/
fitOut
[
1
][
2
][
-
1
]))
if
fitOut
[
1
][
1
]
==
self
.
M
+
1
:
P
=
fitOut
[
0
]
.
T
break
warn
((
"Polyfit is poorly conditioned. Reducing M from {} to "
"{}. Exact snapshot interpolation not guaranteed."
)
\
.
format
(
self
.
M
,
fitOut
[
1
][
1
]
-
1
))
self
.
M
=
fitOut
[
1
][
1
]
-
1
self
.
P
=
np
.
atleast_2d
(
P
)
if
self
.
POD
:
self
.
P
=
self
.
samplingEngine
.
RPOD
.
dot
(
self
.
P
)
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"DEL"
,
"Done computing numerator."
)
self
.
lastApproxParameters
=
copy
(
self
.
approxParameters
)
if
hasattr
(
self
,
"lastSolvedApp"
):
del
self
.
lastSolvedApp
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"DEL"
,
"Done setting up approximant.
\n
"
)
def
getPVal
(
self
,
mu
:
List
[
complex
]):
"""
Evaluate Pade' numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
self
.
setupApprox
()
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"INIT"
,
"Evaluating numerator at mu = {}."
.
format
(
mu
),
end
=
""
)
try
:
len
(
mu
)
except
:
mu
=
[
mu
]
p
=
self
.
_val
(
self
.
radiusPade
(
mu
),
self
.
P
.
T
)
if
len
(
mu
)
==
1
:
p
=
p
.
flatten
()
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"DEL"
,
" Done."
,
inline
=
True
)
return
p
def
getQVal
(
self
,
mu
:
List
[
complex
]):
"""
Evaluate Pade' denominator at arbitrary parameter.
Args:
mu: Target parameter.
"""
self
.
setupApprox
()
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"INIT"
,
"Evaluating denominator at mu = {}."
.
format
(
mu
),
end
=
""
)
q
=
self
.
_val
(
self
.
radiusPade
(
mu
),
self
.
Q
)
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"DEL"
,
" Done."
,
inline
=
True
)
return
q
def
getPoles
(
self
)
->
Np1D
:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self
.
setupApprox
()
return
self
.
HFEngine
.
rescalingInv
(
self
.
scaleFactor
*
self
.
_roots
(
self
.
Q
)
+
self
.
HFEngine
.
rescaling
(
self
.
mu0
))
Event Timeline
Log In to Comment