Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F60515539
generic_approximant_lagrange.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
Tue, Apr 30, 19:11
Size
11 KB
Mime Type
text/x-python
Expires
Thu, May 2, 19:11 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17365624
Attached To
R6746 RationalROMPy
generic_approximant_lagrange.py
View Options
#!/usr/bin/python
# 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/>.
#
import
numpy
as
np
import
warnings
from
rrompy.reduction_methods.base.generic_approximant
import
GenericApproximant
from
rrompy.utilities.types
import
Np1D
,
ListAny
,
DictAny
,
HFEng
,
HSEng
__all__
=
[
'GenericApproximantLagrange'
]
class
GenericApproximantLagrange
(
GenericApproximant
):
"""
ROM Lagrange interpolant computation for parametric problems (ABSTRACT).
Args:
HFEngine: HF problem solver. Should have members:
- 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 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.
mus(optional): Array of snapshot parameters. Defaults to np.array([0]).
ws(optional): Array of snapshot weigths (possibly unused). Defaults to
np.ones_like(self.mus).
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 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: Array of snapshot parameters.
ws: Array of snapshot weigths (possibly unused).
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 snapshots current approximant relies upon.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
autoNodeTypes: List of possible values for autoNode.
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.
Possibly unused.
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.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
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.
"""
autoNodeTypes
=
[
"UNIFORM"
,
"CHEBYSHEV"
,
"GAUSSLEGENDRE"
]
def
__init__
(
self
,
HFEngine
:
HFEng
,
HSEngine
:
HSEng
,
mus
:
Np1D
=
np
.
array
([
0
]),
ws
:
Np1D
=
None
,
approxParameters
:
DictAny
=
{},
plotSnap
:
ListAny
=
[],
plotSnapSpecs
:
DictAny
=
{}):
self
.
mus
=
mus
if
ws
is
None
:
ws
=
np
.
ones_like
(
self
.
mus
)
self
.
ws
=
ws
super
()
.
__init__
(
self
,
HFEngine
=
HFEngine
,
HSEngine
=
HSEngine
,
mu0
=
self
.
mu0
,
approxParameters
=
approxParameters
)
self
.
parameterList
+=
[
"S"
]
self
.
plotSnap
=
plotSnap
self
.
plotSnapSpecs
=
plotSnapSpecs
self
.
resetSamples
()
@property
def
mu0
(
self
):
"""Dummy center of approximant (i.e. mu0)."""
self
.
mu0
=
np
.
mean
(
self
.
mus
)
return
self
.
_mu0
@mu0.setter
def
mu0
(
self
,
mu0
:
bool
):
self
.
_mu0
=
mu0
@property
def
mus
(
self
):
"""Value of mus. Its assignment may reset snapshots."""
return
self
.
_mus
@mus.setter
def
mus
(
self
,
mus
):
musOld
=
self
.
mus
if
hasattr
(
self
,
'mus'
)
else
None
self
.
_mus
=
np
.
array
(
mus
)
_
,
musCounts
=
np
.
unique
(
self
.
_mus
,
return_counts
=
True
)
if
len
(
np
.
where
(
musCounts
>
1
)[
0
])
>
0
:
raise
Exception
(
"Repeated sample points not allowed."
)
if
(
musOld
is
None
or
len
(
self
.
mus
)
!=
len
(
musOld
)
or
not
np
.
allclose
(
self
.
mus
,
musOld
,
1e-14
)):
self
.
resetSamples
()
self
.
autoNode
=
None
@property
def
ws
(
self
):
"""Value of ws."""
return
self
.
_mus
@ws.setter
def
ws
(
self
,
ws
):
if
hasattr
(
self
,
'ws'
):
wsOld
=
self
.
ws
else
:
wsOld
=
None
ws
=
np
.
array
(
ws
)
self
.
_ws
=
np
.
abs
(
ws
)
if
not
np
.
array_equal
(
self
.
_ws
,
ws
):
warnings
.
warn
((
"Negative weights not allowed. Overriding to "
"absolute values."
),
stacklevel
=
2
)
if
(
wsOld
is
None
or
len
(
self
.
ws
)
!=
len
(
wsOld
)
or
not
np
.
allclose
(
self
.
ws
,
wsOld
,
1e-14
)):
self
.
autoNode
=
None
def
checkNodesWeightsConsistency
(
self
)
->
bool
:
"""Check if mus and ws have consistent lengths."""
return
self
.
mus
.
shape
==
self
.
ws
.
shape
@property
def
approxParameters
(
self
):
"""Value of approximant parameters. Its assignment may change S."""
return
self
.
_approxParameters
@approxParameters.setter
def
approxParameters
(
self
,
approxParams
):
approxParameters
=
super
()
.
approxParameters
.
fset
(
self
,
approxParams
)
keyList
=
list
(
approxParameters
.
keys
())
if
"S"
in
keyList
:
self
.
S
=
approxParameters
[
"S"
]
elif
hasattr
(
self
,
"S"
):
self
.
S
=
self
.
S
else
:
self
.
S
=
2
return
approxParameters
@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
self
.
_S
=
S
self
.
_approxParameters
[
"S"
]
=
self
.
S
if
Sold
!=
self
.
S
:
self
.
resetSamples
()
def
resetSamples
(
self
):
"""Reset samples."""
self
.
solSnapshots
=
None
self
.
RPOD
=
None
def
setNodesWeights
(
self
,
xs
:
Np1D
=
np
.
array
([
0
,
1
]),
nodeTypes
:
str
=
"UNIFORM"
):
"""
Assign sampling nodes and weights based on a quadrature formula.
Args:
xs: extremes of the sampling interval;
nodeTypes: type of nodes/weights; allowed values are in
self.autoNodeTypes.
"""
if
len
(
xs
)
!=
2
:
raise
Exception
((
"Value of xs not recognized. Should be numpy 1D "
"vector or list with length 2."
))
try
:
nodeTypes
=
nodeTypes
.
upper
()
.
strip
()
.
replace
(
" "
,
""
)
if
nodeTypes
not
in
self
.
autoNodeTypes
:
raise
except
:
warnings
.
warn
((
"Prescribed nodeTypes not recognized. Overriding "
"to 'UNIFORM'."
),
stacklevel
=
2
)
nodeTypes
=
"UNIFORM"
if
nodeTypes
==
"UNIFORM"
:
nodes
=
np
.
linspace
(
1.
,
-
1.
,
self
.
S
)
weights
=
np
.
ones
((
self
.
S
,
1
))
elif
nodeTypes
==
"CHEBYSHEV"
:
nodes
,
weights
=
np
.
polynomial
.
chebyshev
.
chebgauss
(
self
.
S
)
weights
=
weights
[:,
None
]
*
2.
/
np
.
pi
else
:
#if nodeTypes == "GAUSSLEGENDRE":
nodes
,
weights
=
np
.
polynomial
.
legendre
.
leggauss
(
self
.
S
)
weights
=
weights
[:,
None
]
mu0
=
np
.
mean
(
xs
)
r
=
(
xs
[
0
]
-
xs
[
1
])
/
2.
self
.
mus
=
r
*
nodes
+
mu0
self
.
ws
=
np
.
abs
(
r
)
*
weights
self
.
autoNode
=
nodeTypes
def
computeSnapshots
(
self
):
"""
Compute snapshots of solution map.
"""
if
self
.
solSnapshots
is
None
:
if
len
(
self
.
mus
)
!=
self
.
S
:
warnings
.
warn
((
"Number of prescribed nodes different from "
"S. Overriding S to len(mus)."
),
stacklevel
=
2
)
self
.
S
=
len
(
self
.
mus
)
for
j
,
mu
in
enumerate
(
self
.
mus
):
self
.
solveHF
(
mu
)
self
.
HSEngine
.
plot
(
self
.
uHF
,
name
=
"u({:.4f})"
.
format
(
mu
),
what
=
self
.
plotSnap
,
**
self
.
plotSnapSpecs
)
self
.
manageSnapshots
(
self
.
uHF
,
j
)
def
manageSnapshots
(
self
,
u
:
Np1D
,
pos
:
int
):
"""
Store snapshots of solution map.
Args:
u: solution derivative as numpy complex vector;
pos: Derivative index.
"""
if
pos
==
0
:
self
.
solSnapshots
=
np
.
empty
((
u
.
shape
[
0
],
self
.
S
),
dtype
=
np
.
complex
)
self
.
solSnapshots
[:,
pos
]
=
u
if
self
.
POD
:
if
pos
==
0
:
self
.
RPOD
=
np
.
eye
(
self
.
S
,
dtype
=
np
.
complex
)
beta
=
1
for
j
in
range
(
2
):
#twice is enough!
nu
=
self
.
solSnapshots
[:,
pos
]
.
conj
()
.
dot
(
self
.
energyNormMatrix
.
dot
(
self
.
solSnapshots
[:,
:
pos
])
)
.
conj
()
self
.
RPOD
[:
pos
,
pos
]
=
self
.
RPOD
[:
pos
,
pos
]
+
beta
*
nu
eta
=
(
self
.
solSnapshots
[:,
pos
]
-
self
.
solSnapshots
[:,
:
pos
]
.
dot
(
nu
))
beta
=
eta
.
conj
()
.
dot
(
self
.
energyNormMatrix
.
dot
(
eta
))
**.
5
self
.
solSnapshots
[:,
pos
]
=
eta
/
beta
self
.
RPOD
[
pos
,
pos
]
=
beta
*
self
.
RPOD
[
pos
,
pos
]
def
checkComputedApprox
(
self
)
->
bool
:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return
(
self
.
solSnapshots
is
not
None
and
super
()
.
checkComputedApprox
(
self
))
Event Timeline
Log In to Comment