Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79574008
approximant_lagrange_rb.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, Aug 27, 04:21
Size
12 KB
Mime Type
text/x-python
Expires
Thu, Aug 29, 04:21 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20205132
Attached To
R6746 RationalROMPy
approximant_lagrange_rb.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
import
scipy
as
sp
from
.generic_approximant_lagrange
import
GenericApproximantLagrange
from
rrompy.utilities.base.types
import
Np1D
,
DictAny
,
HFEng
from
rrompy.utilities.base
import
purgeDict
,
verbosityDepth
from
rrompy.utilities.warning_manager
import
warn
__all__
=
[
'ApproximantLagrangeRB'
]
class
ApproximantLagrangeRB
(
GenericApproximantLagrange
):
"""
ROM RB approximant computation for parametric problems.
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];
- 'R': rank for Galerkin projection; defaults to S.
Defaults to empty dict.
homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths (unused).
homogeneized: Whether to homogeneize Dirichlet BCs.
approxRadius: Dummy radius of approximant (i.e. distance from mu0 to
farthest sample point).
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;
- 'R': rank for Galerkin projection.
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.
R: Rank for Galerkin projection.
POD: Whether to compute POD of snapshots.
samplingEngine: Sampling engine.
projMat: Projection matrix for RB system assembly.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
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 theta(mu).
bs: List of numpy vectors representing coefficients of linear system
RHS wrt theta(mu).
thetaAs: List of callables representing coefficients of linear system
matrix wrt mu.
thetabs: List of callables representing coefficients of linear system
RHS wrt mu.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix wrt theta(mu).
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS wrt theta(mu).
"""
def
__init__
(
self
,
HFEngine
:
HFEng
,
mu0
:
complex
=
0.
,
approxParameters
:
DictAny
=
{},
homogeneized
:
bool
=
False
,
verbosity
:
int
=
10
):
self
.
_preInit
()
self
.
_addParametersToList
([
"R"
])
super
()
.
__init__
(
HFEngine
=
HFEngine
,
mu0
=
mu0
,
approxParameters
=
approxParameters
,
homogeneized
=
homogeneized
,
verbosity
=
verbosity
)
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"INIT"
,
"Computing affine blocks of system."
)
self
.
As
,
self
.
thetaAs
=
self
.
HFEngine
.
affineBlocksA
(
self
.
mu0
)
self
.
bs
,
self
.
thetabs
=
self
.
HFEngine
.
affineBlocksb
(
self
.
mu0
,
self
.
homogeneized
)
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"DEL"
,
"Done computing affine blocks."
)
self
.
_postInit
()
def
resetSamples
(
self
):
"""Reset samples."""
super
()
.
resetSamples
()
self
.
projMat
=
None
@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
):
approxParameters
=
purgeDict
(
approxParams
,
self
.
parameterList
,
dictname
=
self
.
name
()
+
".approxParameters"
,
baselevel
=
1
)
approxParametersCopy
=
purgeDict
(
approxParameters
,
[
"R"
],
True
,
True
,
baselevel
=
1
)
GenericApproximantLagrange
.
approxParameters
.
fset
(
self
,
approxParametersCopy
)
keyList
=
list
(
approxParameters
.
keys
())
if
"R"
in
keyList
:
self
.
R
=
approxParameters
[
"R"
]
elif
hasattr
(
self
,
"R"
):
self
.
R
=
self
.
R
else
:
self
.
R
=
self
.
S
@property
def
R
(
self
):
"""Value of R. Its assignment may change S."""
return
self
.
_R
@R.setter
def
R
(
self
,
R
):
if
R
<
0
:
raise
ArithmeticError
(
"R must be non-negative."
)
self
.
_R
=
R
self
.
_approxParameters
[
"R"
]
=
self
.
R
if
hasattr
(
self
,
"S"
)
and
self
.
S
<
self
.
R
:
warn
(
"Prescribed S is too small. Updating S to R."
)
self
.
S
=
self
.
R
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
.
projMat
is
not
None
and
super
()
.
checkComputedApprox
())
def
setupApprox
(
self
):
"""Compute RB projection matrix."""
if
not
self
.
checkComputedApprox
():
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"INIT"
,
"Setting up {}."
.
format
(
self
.
name
()))
self
.
computeSnapshots
()
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
"Computing projection matrix."
,
end
=
""
)
if
self
.
POD
:
U
,
_
,
_
=
np
.
linalg
.
svd
(
self
.
samplingEngine
.
RPOD
,
full_matrices
=
False
)
self
.
projMat
=
self
.
samplingEngine
.
samples
.
dot
(
U
[:,
:
self
.
R
])
else
:
self
.
projMat
=
self
.
samplingEngine
.
samples
[:,
:
self
.
R
]
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"DEL"
,
" Done."
,
inline
=
True
)
self
.
lastApproxParameters
=
copy
(
self
.
approxParameters
)
if
hasattr
(
self
,
"lastSolvedApp"
):
del
self
.
lastSolvedApp
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"DEL"
,
"Done setting up approximant.
\n
"
)
def
assembleReducedSystem
(
self
):
"""Build affine blocks of RB linear system through projections."""
self
.
setupApprox
()
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"INIT"
,
"Projecting affine terms of HF model."
,
end
=
""
)
projMatH
=
self
.
projMat
.
T
.
conj
()
self
.
ARBs
=
[
None
]
*
len
(
self
.
As
)
self
.
bRBs
=
[
None
]
*
len
(
self
.
bs
)
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"MAIN"
,
"."
,
end
=
""
,
inline
=
True
)
for
j
in
range
(
len
(
self
.
As
)):
self
.
ARBs
[
j
]
=
projMatH
.
dot
(
self
.
As
[
j
]
.
dot
(
self
.
projMat
))
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"MAIN"
,
"."
,
end
=
""
,
inline
=
True
)
for
j
in
range
(
len
(
self
.
bs
)):
self
.
bRBs
[
j
]
=
projMatH
.
dot
(
self
.
bs
[
j
])
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"DEL"
,
"Done."
,
inline
=
True
)
def
solveReducedSystem
(
self
,
mu
:
complex
)
->
Np1D
:
"""
Solve RB linear system.
Args:
mu: Target parameter.
Returns:
Solution of RB linear system.
"""
self
.
assembleReducedSystem
()
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"INIT"
,
"Assembling reduced model for mu = {}."
.
format
(
mu
),
end
=
""
)
ARBmu
=
self
.
thetaAs
(
mu
,
0
)
*
self
.
ARBs
[
0
][:
self
.
R
,:
self
.
R
]
bRBmu
=
self
.
thetabs
(
mu
,
0
)
*
self
.
bRBs
[
0
][:
self
.
R
]
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"MAIN"
,
"."
,
end
=
""
,
inline
=
True
)
for
j
in
range
(
1
,
len
(
self
.
ARBs
)):
ARBmu
+=
self
.
thetaAs
(
mu
,
j
)
*
self
.
ARBs
[
j
][:
self
.
R
,
:
self
.
R
]
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"MAIN"
,
"."
,
end
=
""
,
inline
=
True
)
for
j
in
range
(
1
,
len
(
self
.
bRBs
)):
bRBmu
+=
self
.
thetabs
(
mu
,
j
)
*
self
.
bRBs
[
j
][:
self
.
R
]
if
self
.
verbosity
>=
10
:
verbosityDepth
(
"DEL"
,
"Done."
,
inline
=
True
)
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"INIT"
,
"Solving reduced model for mu = {}."
.
format
(
mu
),
end
=
""
)
uRB
=
np
.
linalg
.
solve
(
ARBmu
,
bRBmu
)
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"DEL"
,
" Done."
,
inline
=
True
)
return
uRB
def
evalApproxReduced
(
self
,
mu
:
complex
):
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
"""
self
.
assembleReducedSystem
()
if
(
not
hasattr
(
self
,
"lastSolvedApp"
)
or
not
np
.
isclose
(
self
.
lastSolvedApp
,
mu
)):
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"INIT"
,
"Computing RB solution at mu = {}."
.
format
(
mu
))
self
.
uAppReduced
=
self
.
solveReducedSystem
(
mu
)
self
.
lastSolvedApp
=
mu
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"DEL"
,
"Done computing RB solution."
)
def
evalApprox
(
self
,
mu
:
complex
):
"""
Evaluate approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
self
.
evalApproxReduced
(
mu
)
self
.
uApp
=
self
.
projMat
[:,
:
self
.
R
]
.
dot
(
self
.
uAppReduced
)
def
getPoles
(
self
)
->
Np1D
:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
warn
((
"Impossible to compute poles in general affine parameter "
"dependence. Results subject to interpretation/rescaling, or "
"possibly completely wrong."
))
self
.
assembleReducedSystem
()
if
len
(
self
.
ARBs
)
<
2
:
return
A
=
np
.
eye
(
self
.
ARBs
[
0
]
.
shape
[
0
]
*
(
len
(
self
.
ARBs
)
-
1
),
dtype
=
np
.
complex
)
B
=
np
.
zeros_like
(
A
)
A
[:
self
.
ARBs
[
0
]
.
shape
[
0
],
:
self
.
ARBs
[
0
]
.
shape
[
0
]]
=
-
self
.
ARBs
[
0
]
for
j
in
range
(
len
(
self
.
ARBs
)
-
1
):
Aj
=
self
.
ARBs
[
j
+
1
]
B
[:
Aj
.
shape
[
0
],
j
*
Aj
.
shape
[
0
]
:
(
j
+
1
)
*
Aj
.
shape
[
0
]]
=
Aj
II
=
np
.
arange
(
self
.
ARBs
[
0
]
.
shape
[
0
],
self
.
ARBs
[
0
]
.
shape
[
0
]
*
(
len
(
self
.
ARBs
)
-
1
))
B
[
II
,
II
-
self
.
ARBs
[
0
]
.
shape
[
0
]]
=
1.
return
self
.
HFEngine
.
rescalingInv
(
sp
.
linalg
.
eigvals
(
A
,
B
)
+
self
.
HFEngine
.
rescaling
(
self
.
mu0
))
Event Timeline
Log In to Comment