Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78408038
ROMApproximantLagrangeRB.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 20, 16:36
Size
12 KB
Mime Type
text/x-python
Expires
Thu, Aug 22, 16:36 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20041457
Attached To
R6746 RationalROMPy
ROMApproximantLagrangeRB.py
View Options
#!/usr/bin/python
from
copy
import
copy
import
warnings
import
numpy
as
np
import
scipy
as
sp
import
utilities
from
ROMApproximantLagrange
import
ROMApproximantLagrange
PI
=
np
.
pi
class
ROMApproximantLagrangeRB
(
ROMApproximantLagrange
):
"""
ROM RB approximant computation for parametric problems.
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 k);
- setProblemData: set HF problem data and k to prescribed values;
- getLSBlocks: get algebraic system blocks;
- liftDirichletData: perform lifting of Dirichlet BC as numpy
complex 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.
ks(optional): Array of snapshot parameters. Defaults to np.array([0]).
ws(optional): Array of snapshots weights. Defaults to uniform = 1.
w(optional): Weight for norm computation. If None, set to
Re(np.mean(ks)). Defaults to None.
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;
- 'R': rank for Galerkin projection; defaults to S;
- 'nodesType': sampling node type; available values are:
- 'CHEBYSHEV': Chebyshev nodes;
- 'GAUSSLEGENDRE': Gauss-Legendre nodes;
- 'MANUAL': manual selection;
defaults to 'MANUAL'.
Defaults to empty dict.
plotSnap(optional): What to plot of snapshots of the Helmholtz
solution map. See plot method in HSEngine. Defaults to [].
Attributes:
HFEngine: HF problem solver.
HSEngine: Hilbert space general purpose engine.
solSnapshots: Array whose columns are FE dofs of snapshots of solution
map.
k0: Default parameter.
ks: Array of snapshot parameters.
ws: Array of snapshots weights.
w: Weight for norm computation.
approxRadius: Dummy radius of approximant (i.e. distance from k0 to
farthest sample point).
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of samples current approximant relies upon;
- 'R': rank for Galerkin projection;
- 'nodesType': sampling node type.
S: Number of solution snapshots over which current approximant is
based upon.
R: Rank for Galerkin projection.
nodesType: Sampling node type.
plotSnap: What to plot of snapshots of the Helmholtz solution map.
POD: Whether to compute POD of snapshots.
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.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
extraApproxParameters
=
[
"R"
,
"nodesType"
]
nodesTypeParameters
=
[
"CHEBYSHEV"
,
"GAUSSLEGENDRE"
,
"MANUAL"
]
def
__init__
(
self
,
HFEngine
:
'HF solver'
,
HSEngine
:
'HS engine'
,
ks
:
"Numpy 1D array"
=
np
.
array
([
0
]),
ws
:
"Numpy 1D array"
=
None
,
w
:
float
=
None
,
approxParameters
:
dict
=
{},
plotSnap
:
list
=
[]):
ROMApproximantLagrange
.
__init__
(
self
,
HFEngine
,
HSEngine
,
ks
,
w
,
approxParameters
,
plotSnap
)
if
ws
==
None
:
ws
=
np
.
ones
(
np
.
size
(
self
.
ks
))
self
.
ws
=
ws
self
.
solLifting
=
self
.
HFEngine
.
liftDirichletData
()
@property
def
k0
(
self
):
"""Dummy center of approximant (i.e. k0)."""
self
.
k0
=
np
.
mean
(
self
.
ks
)
return
self
.
_k0
@k0.setter
def
k0
(
self
,
k0
:
bool
):
self
.
_k0
=
k0
def
resetSamples
(
self
):
"""Reset samples."""
ROMApproximantLagrange
.
resetSamples
(
self
)
self
.
projMat
=
None
def
parameterList
(
self
)
->
list
:
"""
List containing keys which are allowed in approxParameters.
Returns:
List of strings.
"""
return
(
ROMApproximantLagrange
.
parameterList
(
self
)
+
ROMApproximantLagrangeRB
.
extraApproxParameters
)
@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
):
if
not
hasattr
(
self
,
"approxParameters"
):
self
.
_approxParameters
=
{}
approxParameters
=
utilities
.
purgeDict
(
approxParams
,
ROMApproximantLagrangeRB
.
parameterList
(
self
),
dictname
=
self
.
name
()
+
".approxParameters"
)
keyList
=
list
(
approxParameters
.
keys
())
approxParametersCopy
=
utilities
.
purgeDict
(
approxParameters
,
ROMApproximantLagrangeRB
.
extraApproxParameters
,
True
,
True
)
ROMApproximantLagrange
.
approxParameters
.
fset
(
self
,
approxParametersCopy
)
if
"R"
in
keyList
:
self
.
R
=
approxParameters
[
"R"
]
elif
hasattr
(
self
,
"R"
):
self
.
R
=
self
.
R
else
:
self
.
R
=
self
.
S
if
"nodesType"
in
keyList
:
self
.
nodesType
=
approxParameters
[
"nodesType"
]
elif
hasattr
(
self
,
"nodesType"
):
self
.
nodesType
=
self
.
nodesType
else
:
self
.
nodesType
=
"MANUAL"
@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
:
warnings
.
warn
(
"Prescribed S is too small. Updating S to R."
)
self
.
S
=
self
.
R
@property
def
nodesType
(
self
):
"""Value of nodesType."""
return
self
.
_nodesType
@nodesType.setter
def
nodesType
(
self
,
nodesType
):
if
hasattr
(
self
,
"nodesType"
):
nodesTypeold
=
self
.
nodesType
else
:
nodesTypeold
=
-
1
try
:
nodesType
=
nodesType
.
upper
()
.
strip
()
.
replace
(
" "
,
""
)
if
nodesType
not
in
self
.
nodesTypeParameters
:
raise
except
:
warnings
.
warn
((
"Prescribed nodesType not recognized. Overriding to"
" 'MANUAL'."
))
nodesType
=
"MANUAL"
self
.
_nodesType
=
nodesType
self
.
_approxParameters
[
"nodesType"
]
=
self
.
nodesType
if
nodesTypeold
!=
self
.
nodesType
:
self
.
resetSamples
()
def
manageSnapshots
(
self
,
u
:
"2-tuple of Fenics function"
,
pos
:
int
):
"""
Post-process snapshots of solution map.
Any specialization should include something like
self.solSnapshots[:, pos] = (np.array(u[0].vector()[:])
+ 1.j * np.array(u[1].vector()[:]))
Args:
u: 2-tuple containing real and imaginary parts of FE dofs of
snapshot.
pos: Derivative index.
"""
return
ROMApproximantLagrange
.
manageSnapshots
(
self
,
u
-
self
.
solLifting
,
pos
)
def
setupApprox
(
self
):
"""
Compute RB projection matrix.
See ``Householder triangularization of a quasimatrix'', L.Trefethen,
2008 for QR algorithm.
"""
if
self
.
solSnapshots
is
None
:
if
self
.
nodesType
==
"MANUAL"
and
len
(
self
.
ks
)
!=
self
.
S
:
warnings
.
warn
((
"Number of prescribed nodes different from S. "
"Overriding S to len(ks)"
))
self
.
S
=
len
(
self
.
ks
)
self
.
approxRadius
=
np
.
max
(
np
.
abs
(
self
.
k0
-
self
.
ks
))
if
self
.
nodesType
!=
"MANUAL"
:
if
self
.
nodesType
==
"CHEBYSHEV"
:
nodes
,
weights
=
np
.
polynomial
.
chebyshev
.
chebgauss
(
self
.
S
)
self
.
ws
=
weights
*
2.
/
PI
elif
self
.
nodesType
==
"GAUSSLEGENDRE"
:
nodes
,
weights
=
np
.
polynomial
.
legendre
.
leggauss
(
self
.
S
)
self
.
ws
=
weights
self
.
ks
=
nodes
*
self
.
approxRadius
+
self
.
k0
self
.
ws
=
self
.
ws
[:,
None
]
self
.
computeSnapshots
()
if
not
self
.
checkComputedApprox
():
if
self
.
POD
:
U
,
_
,
_
=
np
.
linalg
.
svd
(
self
.
RPOD
,
full_matrices
=
False
)
self
.
projMat
=
self
.
solSnapshots
.
dot
(
U
[:,
:
self
.
R
])
else
:
self
.
projMat
=
self
.
solSnapshots
[:,
:
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
]
*
max
(
len
(
self
.
As
),
len
(
self
.
bs
))
for
j
in
range
(
len
(
self
.
As
)):
self
.
ARBs
[
j
]
=
projMatH
.
dot
(
self
.
As
[
j
]
.
dot
(
self
.
projMat
))
if
j
<
len
(
self
.
bs
):
self
.
bRBs
[
j
]
=
projMatH
.
dot
(
self
.
bs
[
j
]
-
self
.
As
[
j
]
.
dot
(
self
.
solLifting
))
else
:
self
.
bRBs
[
j
]
=
-
projMatH
.
dot
(
self
.
As
[
j
]
.
dot
(
self
.
solLifting
))
for
j
in
range
(
len
(
self
.
As
),
len
(
self
.
bs
)):
self
.
bRBs
[
j
]
=
projMatH
.
dot
(
self
.
bs
[
j
])
def
solveReducedSystem
(
self
,
k
:
complex
)
->
"Numpy 1D array"
:
"""
Solve RB linear system.
Args:
k: Target wavenumber.
Returns:
Solution of RB linear system.
"""
self
.
setupApprox
()
ARBk
=
self
.
ARBs
[
0
][:
self
.
R
,:
self
.
R
]
bRBk
=
self
.
bRBs
[
0
][:
self
.
R
]
for
j
in
range
(
1
,
len
(
self
.
ARBs
)):
ARBk
=
ARBk
+
np
.
power
(
k
,
j
)
*
self
.
ARBs
[
j
][:
self
.
R
,
:
self
.
R
]
for
j
in
range
(
1
,
len
(
self
.
bRBs
)):
bRBk
=
bRBk
+
np
.
power
(
k
,
j
)
*
self
.
bRBs
[
j
][:
self
.
R
]
return
self
.
projMat
[:,
:
self
.
R
]
.
dot
(
np
.
linalg
.
solve
(
ARBk
,
bRBk
))
def
evalApprox
(
self
,
k
:
complex
)
->
(
"Fenics function"
,
"Fenics function"
):
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
k: Target wavenumber.
Returns:
Real and imaginary parts of approximant.
"""
self
.
setupApprox
()
self
.
uApp
=
self
.
solLifting
+
self
.
solveReducedSystem
(
k
)
return
self
.
uApp
def
getPoles
(
self
,
centered
:
bool
=
False
)
->
"numpy 1D array"
:
"""
Obtain approximant poles.
Args:
centered(optional): Whether to return pole positions relative to
approximation center. Defaults to False.
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
)
-
self
.
k0
*
(
not
centered
)
except
:
warnings
.
warn
(
"Generalized eig algorithm did not converge."
)
x
=
np
.
empty
(
A
.
shape
[
0
])
x
[:]
=
np
.
NaN
return
x
Event Timeline
Log In to Comment