Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F60472979
approximant_lagrange_greedy_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, Apr 30, 12:02
Size
15 KB
Mime Type
text/x-python
Expires
Thu, May 2, 12:02 (2 d)
Engine
blob
Format
Raw Data
Handle
17356475
Attached To
R6746 RationalROMPy
approximant_lagrange_greedy_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/>.
#
import
numpy
as
np
from
copy
import
copy
from
.generic_approximant_lagrange_greedy
import
(
GenericApproximantLagrangeGreedy
)
from
rrompy.reduction_methods.lagrange
import
ApproximantLagrangeRB
from
rrompy.utilities.base.types
import
DictAny
,
HFEng
,
List
from
rrompy.utilities.base
import
verbosityDepth
from
rrompy.utilities.warning_manager
import
warn
__all__
=
[
'ApproximantLagrangeRBGreedy'
]
class
ApproximantLagrangeRBGreedy
(
GenericApproximantLagrangeGreedy
,
ApproximantLagrangeRB
):
"""
ROM greedy 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;
- 'muBounds': list of bounds for parameter values; defaults to
[[0], [1]];
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'refinementRatio': ratio of training points to be exhausted
before training set refinement; defaults to 0.2;
- 'nTrainingPoints': number of training points; defaults to
maxIter / refinementRatio;
- 'nTestPoints': number of starting test points; defaults to 1;
- 'trainingSetGenerator': training sample points generator;
defaults to uniform sampler within muBounds.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
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;
- 'muBounds': list of bounds for parameter values;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of training points to be exhausted
before training set refinement;
- 'nTrainingPoints': number of training points;
- 'nTestPoints': number of starting test points;
- 'trainingSetGenerator': training sample points generator;
- 'robustTol': tolerance for robust Pade' denominator management.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
POD: whether to compute POD of snapshots.
muBounds: list of bounds for parameter values.
greedyTol: uniform error tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTrainingPoints: number of training points.
nTestPoints: number of starting test points.
trainingSetGenerator: training sample points generator.
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
()
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
()
@property
def
S
(
self
):
"""Value of S."""
return
self
.
_S
@S.setter
def
S
(
self
,
S
):
self
.
_S
=
S
@property
def
R
(
self
):
"""Value of R."""
return
self
.
_S
@R.setter
def
R
(
self
,
R
):
warn
((
"R is used just to simplify inheritance, and its value cannot "
"be changed from that of S."
))
def
resetSamples
(
self
):
"""Reset samples."""
super
()
.
resetSamples
()
self
.
projMat
=
None
self
.
resbb
=
None
self
.
resAb
=
None
self
.
resAA
=
None
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
errorEstimator
(
self
,
mus
:
List
[
np
.
complex
],
homogeneized
:
bool
=
None
)
->
List
[
np
.
complex
]:
"""
Standard residual-based error estimator. Unreliable for unstable
problems (inf-sup constant is missing).
"""
if
homogeneized
is
None
:
homogeneized
=
self
.
homogeneized
nmus
=
len
(
mus
)
err
=
[
0.
]
*
nmus
self
.
assembleReducedSystem
()
self
.
assembleReducedResidualBlocks
(
homogeneized
)
assert
homogeneized
==
self
.
homogeneized
nAs
=
len
(
self
.
As
)
nbs
=
len
(
self
.
bs
)
for
j
in
range
(
nmus
):
mu
=
mus
[
j
]
uAppReduced
=
self
.
getAppReduced
(
mu
)
prodbb
=
0.
prodAb
=
0.
prodAA
=
0.
for
i1
in
range
(
nbs
):
rhobi1
=
self
.
thetabs
(
mu
,
i1
)
for
i2
in
range
(
nbs
):
rhobi2
=
self
.
thetabs
(
mu
,
i2
)
.
conj
()
prodbb
+=
(
rhobi1
*
rhobi2
*
self
.
resbb
[
homogeneized
][
i1
][
i2
])
for
i1
in
range
(
nAs
):
rhoAi1
=
self
.
thetaAs
(
mu
,
i1
)
for
i2
in
range
(
nbs
):
rhobi2
=
self
.
thetabs
(
mu
,
i2
)
.
conj
()
prodAb
+=
(
rhoAi1
*
rhobi2
*
self
.
resAb
[
homogeneized
][
i1
][
i2
])
for
i1
in
range
(
nAs
):
rhoAi1
=
self
.
thetaAs
(
mu
,
i1
)
for
i2
in
range
(
nAs
):
rhoAi2
=
self
.
thetaAs
(
mu
,
i2
)
.
conj
()
prodAA
+=
(
rhoAi1
*
rhoAi2
*
self
.
resAA
[
homogeneized
][
i1
][
i2
])
err
[
j
]
=
np
.
abs
(((
uAppReduced
.
T
.
conj
()
.
dot
(
prodAA
.
dot
(
uAppReduced
))
-
2.
*
np
.
real
(
prodAb
.
dot
(
uAppReduced
)))
/
prodbb
+
1.
)[
0
])
**
.
5
return
err
def
setupApprox
(
self
):
"""Compute RB projection matrix."""
if
not
self
.
checkComputedApprox
():
if
self
.
verbosity
>=
5
:
verbosityDepth
(
"INIT"
,
"Setting up {}."
.
format
(
self
.
name
()))
self
.
greedy
()
self
.
S
=
len
(
self
.
mus
)
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
"Computing projection matrix."
,
end
=
""
)
self
.
projMat
=
self
.
samplingEngine
.
samples
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
assembleReducedResidualBlocks
(
self
,
homogeneized
:
bool
):
"""Build affine blocks of RB linear system through projections."""
self
.
initMassMatrixInverse
()
self
.
assembleReducedSystem
()
computeResbb
=
(
self
.
resbb
is
None
or
homogeneized
not
in
self
.
resbb
.
keys
())
computeResAb
=
(
self
.
resAb
is
None
or
homogeneized
not
in
self
.
resAb
.
keys
()
or
self
.
resAb
[
homogeneized
][
0
][
0
]
.
shape
[
0
]
!=
self
.
S
)
computeResAA
=
(
self
.
resAA
is
None
or
homogeneized
not
in
self
.
resAA
.
keys
()
or
self
.
resAA
[
homogeneized
][
0
][
0
]
.
shape
[
0
]
!=
self
.
S
)
if
computeResbb
or
computeResAb
or
computeResAA
:
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"INIT"
,
"Projecting affine terms of residual."
,
end
=
""
)
nAs
=
len
(
self
.
As
)
nbs
=
len
(
self
.
bs
)
if
computeResbb
:
if
self
.
resbb
is
None
:
self
.
resbb
=
{}
self
.
resbb
[
homogeneized
]
=
[]
for
i
in
range
(
nbs
):
resbbi
=
[]
Mbi
=
self
.
massInv
.
solve
(
self
.
bs
[
i
])
for
j
in
range
(
i
):
Mbj
=
self
.
massInv
.
solve
(
self
.
bs
[
j
])
resbbi
+=
[
self
.
HFEngine
.
innerProduct
(
Mbi
,
Mbj
)]
resbbi
+=
[
self
.
HFEngine
.
innerProduct
(
Mbi
,
Mbi
)]
self
.
resbb
[
homogeneized
]
+=
[
resbbi
]
for
i
in
range
(
nbs
):
for
j
in
range
(
i
+
1
,
nbs
):
self
.
resbb
[
homogeneized
][
i
]
+=
[
self
.
resbb
[
homogeneized
][
j
][
i
]
.
conj
()]
if
computeResAb
:
if
self
.
resAb
is
None
:
self
.
resAb
=
{}
if
homogeneized
not
in
self
.
resAb
.
keys
():
self
.
resAb
[
homogeneized
]
=
[]
for
i
in
range
(
nAs
):
resAbi
=
[]
MAi
=
self
.
massInv
.
solve
(
self
.
As
[
i
]
.
dot
(
self
.
projMat
))
for
j
in
range
(
nbs
):
Mbj
=
self
.
massInv
.
solve
(
self
.
bs
[
j
])
resAbi
+=
[
self
.
HFEngine
.
innerProduct
(
MAi
,
Mbj
)]
self
.
resAb
[
homogeneized
]
+=
[
resAbi
]
else
:
Sold
=
self
.
resAb
[
homogeneized
][
0
][
0
]
.
shape
[
0
]
for
i
in
range
(
nAs
):
MAi
=
self
.
massInv
.
solve
(
self
.
As
[
i
]
.
dot
(
self
.
projMat
[:,
Sold
:]))
for
j
in
range
(
nbs
):
Mbj
=
self
.
massInv
.
solve
(
self
.
bs
[
j
])
self
.
resAb
[
homogeneized
][
i
][
j
]
=
np
.
append
(
self
.
resAb
[
homogeneized
][
i
][
j
],
self
.
HFEngine
.
innerProduct
(
MAi
,
Mbj
))
if
computeResAA
:
if
self
.
resAA
is
None
:
self
.
resAA
=
{}
if
homogeneized
not
in
self
.
resAA
.
keys
():
self
.
resAA
[
homogeneized
]
=
[]
for
i
in
range
(
nAs
):
resAAi
=
[]
MAi
=
self
.
massInv
.
solve
(
self
.
As
[
i
]
.
dot
(
self
.
projMat
))
for
j
in
range
(
i
):
MAj
=
self
.
massInv
.
solve
(
self
.
As
[
j
]
.
dot
(
self
.
projMat
))
resAAi
+=
[
self
.
HFEngine
.
innerProduct
(
MAi
,
MAj
)]
resAAi
+=
[
self
.
HFEngine
.
innerProduct
(
MAi
,
MAi
)]
self
.
resAA
[
homogeneized
]
+=
[
resAAi
]
for
i
in
range
(
nAs
):
for
j
in
range
(
i
+
1
,
nAs
):
self
.
resAA
[
homogeneized
][
i
]
+=
[
self
.
resAA
[
homogeneized
][
j
][
i
]
.
conj
()]
else
:
Sold
=
self
.
resAA
[
homogeneized
][
0
][
0
]
.
shape
[
0
]
for
i
in
range
(
nAs
):
MAi
=
self
.
massInv
.
solve
(
self
.
As
[
i
]
.
dot
(
self
.
projMat
))
for
j
in
range
(
i
):
MAj
=
self
.
massInv
.
solve
(
self
.
As
[
j
]
.
dot
(
self
.
projMat
))
resAAcross1
=
self
.
HFEngine
.
innerProduct
(
MAi
[:,
Sold
:],
MAj
[:,
:
Sold
])
resAAcross2
=
self
.
HFEngine
.
innerProduct
(
MAi
[:,
:
Sold
],
MAj
[:,
Sold
:])
resAAdiag
=
self
.
HFEngine
.
innerProduct
(
MAi
[:,
Sold
:],
MAj
[:,
Sold
:])
self
.
resAA
[
homogeneized
][
i
][
j
]
=
np
.
block
(
[[
self
.
resAA
[
homogeneized
][
i
][
j
],
resAAcross1
],
[
resAAcross2
,
resAAdiag
]])
resAAcross
=
self
.
HFEngine
.
innerProduct
(
MAi
[:,
Sold
:],
MAi
[:,
:
Sold
])
resAAdiag
=
self
.
HFEngine
.
innerProduct
(
MAi
[:,
Sold
:],
MAi
[:,
Sold
:])
self
.
resAA
[
homogeneized
][
i
][
i
]
=
np
.
block
(
[[
self
.
resAA
[
homogeneized
][
i
][
i
],
resAAcross
],
[
resAAcross
.
T
.
conj
(),
resAAdiag
]])
for
i
in
range
(
nAs
):
for
j
in
range
(
i
+
1
,
nAs
):
self
.
resAA
[
homogeneized
][
i
][
j
]
=
(
self
.
resAA
[
homogeneized
][
j
][
i
]
.
conj
())
if
self
.
verbosity
>=
7
:
verbosityDepth
(
"DEL"
,
(
"Done setting up affine "
"decomposition of residual."
))
Event Timeline
Log In to Comment