Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62268437
approximant_taylor_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
Sun, May 12, 01:53
Size
13 KB
Mime Type
text/x-python
Expires
Tue, May 14, 01:53 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17626558
Attached To
R6746 RationalROMPy
approximant_taylor_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
warnings
import
numpy
as
np
import
scipy
as
sp
from
rrompy.reduction_methods.taylor.generic_approximant_taylor
import
GenericApproximantTaylor
from
rrompy.utilities.types
import
Np1D
,
ListAny
,
DictAny
,
HFEng
,
HSEng
from
rrompy.utilities
import
purgeDict
__all__
=
[
'ApproximantTaylorRB'
]
class
ApproximantTaylorRB
(
GenericApproximantTaylor
):
"""
ROM single-point fast RB approximant computation for parametric problems
with polynomial dependence up to degree 2.
Args:
HFEngine: HF problem solver. Should have members:
- energyNormMatrix: sparse matrix (in CSC format) associated to
w-weighted H10 norm;
- problemData: list of HF problem data (excluding mu);
- setProblemData: set HF problem data and mu to prescribed values;
- getLSBlocks: get blocks of HF linear system as sparse matrices in
CSC format;
- liftDirichletData: perform lifting of Dirichlet BC as numpy
complex vector;
- setupDerivativeComputation: setup HF problem data to compute j-th
solution derivative at mu0;
- solve: get HF solution as complex numpy vector.
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.
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;
- 'R': rank for Galerkin projection; defaults to E + 1;
- 'E': total number of derivatives current approximant relies upon;
defaults to Emax;
- 'Emax': total number of derivatives of solution map to be
computed; defaults to E;
- 'sampleType': label of sampling type; available values are:
- 'ARNOLDI': orthogonalization of solution derivatives through
Arnoldi algorithm;
- 'KRYLOV': standard computation of solution derivatives.
Defaults to 'KRYLOV'.
Defaults to empty dict.
plotDer(optional): What to plot of derivatives of the Helmholtz
solution map. See plot method in HSEngine. Defaults to [].
plotDerSpecs(optional): How to plot derivatives 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.
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;
- 'R': rank for Galerkin projection;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'sampleType': label of sampling type.
solDerivatives: Array whose columns are FE dofs of solution
derivatives.
RArnoldi: Right factor of derivatives QR factorization. If no Arnoldi,
set to None.
HArnoldi: Hessenberg factor of derivatives Arnoldi algorithm. If no
Arnoldi, set to None.
POD: Whether to compute QR factorization of derivatives.
R: Rank for Galerkin projection.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
sampleType: Label of sampling type, i.e. 'KRYLOV'.
plotDer: What to plot of derivatives of the Helmholtz solution map.
plotDerSpecs: How to plot derivatives of the Helmholtz solution map.
energyNormMatrix: Sparse matrix (in CSC format) associated 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.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
solLifting: Numpy complex vector with lifting of real part of Dirichlet
boundary datum.
projMat: Numpy matrix representing projection onto RB space.
projMat: Numpy matrix representing projection onto RB space.
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.
ARBs: List of sparse matrices (in CSC format) representing RB
coefficients of linear system matrix wrt mu.
bRBs: List of numpy vectors representing RB coefficients of linear
system RHS wrt mu.
"""
def
__init__
(
self
,
HFEngine
:
HFEng
,
HSEngine
:
HSEng
,
mu0
:
complex
=
0
,
approxParameters
:
DictAny
=
{},
plotDer
:
ListAny
=
[],
plotDerSpecs
:
DictAny
=
{}):
if
not
hasattr
(
self
,
"parameterList"
):
self
.
parameterList
=
[]
self
.
parameterList
+=
[
"R"
]
super
()
.
__init__
(
HFEngine
=
HFEngine
,
HSEngine
=
HSEngine
,
mu0
=
mu0
,
approxParameters
=
approxParameters
,
plotDer
=
plotDer
,
plotDerSpecs
=
plotDerSpecs
)
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"
)
approxParametersCopy
=
purgeDict
(
approxParameters
,
[
"R"
],
True
,
True
)
GenericApproximantTaylor
.
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
.
E
+
1
@property
def
POD
(
self
):
"""Value of POD."""
return
self
.
_POD
@POD.setter
def
POD
(
self
,
POD
):
GenericApproximantTaylor
.
POD
.
fset
(
self
,
POD
)
if
(
hasattr
(
self
,
"sampleType"
)
and
self
.
sampleType
==
"ARNOLDI"
and
not
self
.
POD
):
warnings
.
warn
((
"Arnoldi sampling implicitly forces POD-type "
"derivative management."
),
stacklevel
=
2
)
@property
def
sampleType
(
self
):
"""Value of sampleType."""
return
self
.
_sampleType
@sampleType.setter
def
sampleType
(
self
,
sampleType
):
GenericApproximantTaylor
.
sampleType
.
fset
(
self
,
sampleType
)
if
(
hasattr
(
self
,
"POD"
)
and
not
self
.
POD
and
self
.
sampleType
==
"ARNOLDI"
):
warnings
.
warn
((
"Arnoldi sampling implicitly forces POD-type "
"derivative management."
),
stacklevel
=
2
)
@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
,
"E"
)
and
self
.
E
+
1
<
self
.
R
:
warnings
.
warn
(
"Prescribed E is too small. Updating E to R - 1."
,
stacklevel
=
2
)
self
.
E
=
self
.
R
-
1
def
setupApprox
(
self
):
"""
Setup RB system. For usage of correlation matrix in POD see Section
6.3.1 in 'Reduced Basis Methods for PDEs. An introduction', A.
Quarteroni, A. Manzoni, F. Negri, Springer, 2016.
"""
if
not
self
.
checkComputedApprox
():
self
.
computeDerivatives
()
if
self
.
POD
and
not
self
.
sampleType
==
"ARNOLDI"
:
A
=
copy
(
self
.
solDerivatives
)
M
=
self
.
energyNormMatrix
V
=
np
.
zeros_like
(
A
,
dtype
=
np
.
complex
)
Q
=
np
.
zeros_like
(
A
,
dtype
=
np
.
complex
)
R
=
np
.
zeros
((
A
.
shape
[
1
],
A
.
shape
[
1
]),
dtype
=
np
.
complex
)
for
k
in
range
(
A
.
shape
[
1
]):
Q
[:,
k
]
=
(
np
.
random
.
randn
(
Q
.
shape
[
0
])
+
1.j
*
np
.
random
.
randn
(
Q
.
shape
[
0
]))
if
k
>
0
:
for
times
in
range
(
2
):
#twice is enough!
Q
[:,
k
]
=
Q
[:,
k
]
-
Q
[:,
:
k
]
.
dot
(
(
Q
[:,
k
]
.
conj
()
.
dot
(
M
.
dot
(
Q
[:,
:
k
])))
.
conj
())
Q
[:,
k
]
=
Q
[:,
k
]
/
(
Q
[:,
k
]
.
conj
()
.
dot
(
M
.
dot
(
Q
[:,
k
])))
**.
5
R
[
k
,
k
]
=
np
.
abs
(
A
[:,
k
]
.
conj
()
.
dot
(
M
.
dot
(
A
[:,
k
])))
**
.
5
alpha
=
Q
[:,
k
]
.
conj
()
.
dot
(
M
.
dot
(
A
[:,
k
]))
if
np
.
isclose
(
np
.
abs
(
alpha
),
0.
):
s
=
1.
else
:
s
=
-
alpha
/
np
.
abs
(
alpha
)
Q
[:,
k
]
=
s
*
Q
[:,
k
]
v
=
R
[
k
,
k
]
*
Q
[:,
k
]
-
A
[:,
k
]
for
times
in
range
(
2
):
#twice is enough!
v
=
v
-
Q
[:,
:
k
]
.
dot
((
v
.
conj
()
.
dot
(
M
.
dot
(
Q
[:,
:
k
]))
)
.
conj
())
sigma
=
np
.
abs
(
v
.
conj
()
.
dot
(
M
.
dot
(
v
)))
**
.
5
if
np
.
isclose
(
sigma
,
0.
):
v
=
Q
[:,
k
]
else
:
v
=
v
/
sigma
V
[:,
k
]
=
v
J
=
np
.
arange
(
k
+
1
,
A
.
shape
[
1
])
vtQ
=
v
.
conj
()
.
dot
(
M
.
dot
(
A
[:,
J
]))
A
[:,
J
]
=
A
[:,
J
]
-
2
*
np
.
outer
(
v
,
vtQ
)
R
[
k
,
J
]
=
Q
[:,
k
]
.
conj
()
.
dot
(
M
.
dot
(
A
[:,
J
]))
A
[:,
J
]
=
A
[:,
J
]
-
np
.
outer
(
Q
[:,
k
],
R
[
k
,
J
])
for
k
in
range
(
A
.
shape
[
1
]
-
1
,
-
1
,
-
1
):
v
=
V
[:,
k
]
J
=
np
.
arange
(
k
,
A
.
shape
[
1
])
vtQ
=
v
.
conj
()
.
dot
(
M
.
dot
(
Q
[:,
J
]))
Q
[:,
J
]
=
Q
[:,
J
]
-
2
*
np
.
outer
(
v
,
vtQ
)
self
.
projMatQ
=
Q
self
.
projMatR
=
R
if
self
.
POD
and
not
self
.
sampleType
==
"ARNOLDI"
:
U
,
_
,
_
=
np
.
linalg
.
svd
(
self
.
projMatR
[:
self
.
R
,
:
self
.
R
])
self
.
projMat
=
self
.
projMatQ
[:,
:
self
.
R
]
.
dot
(
U
)
else
:
self
.
projMat
=
self
.
solDerivatives
[:,
:
self
.
R
]
self
.
assembleReducedSystem
()
self
.
lastApproxParameters
=
copy
(
self
.
approxParameters
)
def
assembleReducedSystem
(
self
):
"""Build affine blocks of RB linear system through projections."""
projMatH
=
self
.
projMat
.
T
.
conjugate
()
self
.
ARBs
=
[
None
]
*
len
(
self
.
As
)
self
.
bRBs
=
[
None
]
*
len
(
self
.
bs
)
for
j
in
range
(
len
(
self
.
As
)):
self
.
ARBs
[
j
]
=
projMatH
.
dot
(
self
.
As
[
j
]
.
dot
(
self
.
projMat
))
for
j
in
range
(
len
(
self
.
bs
)):
self
.
bRBs
[
j
]
=
projMatH
.
dot
(
self
.
bs
[
j
])
def
solveReducedSystem
(
self
,
mu
:
complex
)
->
Np1D
:
"""
Solve RB linear system.
Args:
mu: Target parameter.
Returns:
Solution of RB linear system.
"""
self
.
setupApprox
()
ARBmu
=
self
.
ARBs
[
0
][:
self
.
R
,
:
self
.
R
]
bRBmu
=
self
.
bRBs
[
0
][:
self
.
R
]
for
j
in
range
(
1
,
len
(
self
.
ARBs
)):
ARBmu
=
ARBmu
+
np
.
power
(
mu
,
j
)
*
self
.
ARBs
[
j
][:
self
.
R
,
:
self
.
R
]
for
j
in
range
(
1
,
len
(
self
.
bRBs
)):
bRBmu
=
bRBmu
+
np
.
power
(
mu
,
j
)
*
self
.
bRBs
[
j
][:
self
.
R
]
return
self
.
projMat
[:,
:
self
.
R
]
.
dot
(
np
.
linalg
.
solve
(
ARBmu
,
bRBmu
))
def
evalApprox
(
self
,
mu
:
complex
):
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
"""
self
.
setupApprox
()
self
.
uApp
=
self
.
solveReducedSystem
(
mu
)
def
getPoles
(
self
)
->
Np1D
:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self
.
setupApprox
()
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.
try
:
return
sp
.
linalg
.
eigvals
(
A
,
B
)
except
:
warnings
.
warn
(
"Generalized eig algorithm did not converge."
,
stacklevel
=
2
)
x
=
np
.
empty
(
A
.
shape
[
0
])
x
[:]
=
np
.
NaN
return
x
Event Timeline
Log In to Comment