Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61532726
damped_mass_chain.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, May 7, 06:28
Size
7 KB
Mime Type
text/x-python
Expires
Thu, May 9, 06:28 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17524589
Attached To
R6746 RationalROMPy
damped_mass_chain.py
View Options
### example from Lohmann, Eid. Efficient Order Reduction of Parametric and
### Nonlinear Models by Superposition of Locally Reduced Models.
from
copy
import
deepcopy
as
copy
import
numpy
as
np
import
matplotlib.pyplot
as
plt
from
rrompy.reduction_methods
import
(
NearestNeighbor
as
NN
,
RationalInterpolant
as
RI
,
RationalInterpolantGreedy
as
RIG
,
RationalInterpolantPivoted
as
RIP
,
RationalInterpolantGreedyPivoted
as
RIGP
)
from
rrompy.parameter.parameter_sampling
import
(
QuadratureSampler
as
QS
,
EmptySampler
as
ES
)
from
damped_mass_chain_engine
import
(
bode
as
bode0
,
bodeLog
,
MassChainEngine
,
MassChainEngineLog
,
AugmentedMassChainEngine
,
AugmentedMassChainEngineLog
)
from
rrompy.utilities.base.decorators
import
addWhiteNoise
##########################
fullModelOrder
=
1
#+ 1
SMarginal
=
1
#* 2
state
=
0
#+ 1
noise_level
=
0
#+ 1e-5
LS
=
0
#+ 6
##########################
modelSign
=
"Surrogate modeling for frequency response of "
if
fullModelOrder
==
1
:
modelSign
+=
"augmented "
modelSign
+=
"damper-mass-spring model"
if
SMarginal
>
1
:
modelSign
+=
" with 1 design parameter"
modelSign
+=
".
\n
Output is "
if
state
:
modelSign
+=
"vector of mass displacements. "
else
:
modelSign
+=
"displacement of last mass. "
if
LS
:
modelSign
+=
"Least-squares: S - N - 1 = {}. "
.
format
(
LS
)
else
:
modelSign
+=
"Interpolatory: S = N + 1. "
modelSign
+=
"Noise level: {}.
\n
"
.
format
(
noise_level
)
print
(
modelSign
)
M
=
[
np
.
array
([
1.
,
5.
,
25.
,
125.
])]
N
=
len
(
M
[
0
])
D
=
[
np
.
zeros
((
N
,
N
))]
D
[
0
][
0
,
1
],
D
[
0
][
1
,
2
],
D
[
0
][
2
,
3
],
D
[
0
][
3
,
3
]
=
.
1
,
.
4
,
1.6
,
0.
D
[
0
]
=
D
[
0
]
+
D
[
0
]
.
T
K
=
[
np
.
zeros
((
N
,
N
))]
K
[
0
][
0
,
1
],
K
[
0
][
1
,
2
],
K
[
0
][
2
,
3
],
K
[
0
][
0
,
3
],
K
[
0
][
3
,
3
]
=
9.
,
3.
,
1.
,
1.
,
2.
K
[
0
]
=
K
[
0
]
+
K
[
0
]
.
T
B
=
np
.
append
(
27.
,
np
.
zeros
(
N
-
1
))
.
reshape
(
-
1
,
1
)
if
SMarginal
>
1
:
M
+=
[
np
.
zeros
(
N
)]
D
+=
[
np
.
zeros
((
N
,
N
))]
D
[
1
][
3
,
3
]
=
1.
D
[
1
]
=
D
[
1
]
+
D
[
1
]
.
T
K
+=
[
np
.
zeros
((
N
,
N
))]
K
[
1
][
0
,
3
],
K
[
1
][
3
,
3
]
=
2.
,
2.
K
[
1
]
=
K
[
1
]
+
K
[
1
]
.
T
if
state
:
C
=
np
.
eye
(
4
)
else
:
C
=
np
.
append
(
np
.
zeros
(
N
-
1
),
1.
)
.
reshape
(
1
,
-
1
)
for
logspace
in
range
(
2
):
print
(
"Approximation in l{}space"
.
format
(
"og"
*
logspace
+
"in"
*
(
not
logspace
)))
if
logspace
:
bode
=
bodeLog
if
fullModelOrder
==
1
:
engine
=
AugmentedMassChainEngineLog
else
:
engine
=
MassChainEngineLog
else
:
bode
=
bode0
if
fullModelOrder
==
1
:
engine
=
AugmentedMassChainEngine
else
:
engine
=
MassChainEngine
solver
=
addWhiteNoise
(
noise_level
)(
engine
)(
M
,
D
,
K
,
B
,
C
)
ss
,
mu
=
[
1e-2
,
1e1
],
[]
s0
=
10.
**
np
.
mean
(
np
.
log10
(
ss
))
freq
=
np
.
logspace
(
np
.
log10
(
ss
[
0
]),
np
.
log10
(
ss
[
1
]),
100
)
if
logspace
:
ss
,
freq
=
[
np
.
log10
(
ss
[
0
]),
np
.
log10
(
ss
[
1
])],
np
.
log10
(
freq
)
s0
,
parameterMap
=
np
.
log10
(
s0
),
1.
else
:
parameterMap
=
{
"F"
:
[(
"log10"
,
"x"
)],
"B"
:
[(
10.
,
"**"
,
"x"
)]}
krange
=
[[
ss
[
0
]],
[
ss
[
-
1
]]]
k0
,
srange
=
[
s0
],
copy
(
krange
)
if
SMarginal
>
1
:
ms
=
[
0.
,
1.
]
m0
,
mrange
=
np
.
mean
(
ms
),
[[
ms
[
0
]],
[
ms
[
-
1
]]]
krange
[
0
]
+=
mrange
[
0
]
krange
[
1
]
+=
mrange
[
1
]
k0
+=
[
m0
]
mu
=
[
.
5
*
(
ms
[
1
]
-
ms
[
0
])
/
(
SMarginal
-
1
)]
if
not
logspace
:
parameterMap
[
"F"
]
+=
[(
"x"
)]
parameterMap
[
"B"
]
+=
[(
"x"
)]
for
method
in
[
"RI"
,
"RI_GREEDY"
]:
print
(
"Testing {} method"
.
format
(
method
))
if
method
==
"RI"
:
params
=
{
'S'
:
15
,
'POD'
:
True
,
'polybasis'
:
"CHEBYSHEV"
}
if
LS
:
params
[
"N"
]
=
params
[
"S"
]
-
1
-
LS
if
SMarginal
>
1
:
algo
=
RIP
else
:
params
[
'sampler'
]
=
QS
(
srange
,
"CHEBYSHEV"
,
parameterMap
)
algo
=
RI
if
method
==
"RI_GREEDY"
:
params
=
{
'S'
:
5
,
'POD'
:
True
,
'polybasis'
:
"LEGENDRE"
,
'greedyTol'
:
1e-2
,
'errorEstimatorKind'
:
"DISCREPANCY"
,
'trainSetGenerator'
:
QS
(
srange
,
"CHEBYSHEV"
,
parameterMap
)}
if
SMarginal
>
1
:
algo
=
RIGP
else
:
params
[
'sampler'
]
=
QS
(
srange
,
"UNIFORM"
,
parameterMap
)
algo
=
RIG
if
SMarginal
>
1
:
params
[
"paramsMarginal"
]
=
{
"MMarginal"
:
SMarginal
-
1
}
params
[
'SMarginal'
]
=
SMarginal
params
[
'polybasisMarginal'
]
=
"MONOMIAL"
params
[
'radialDirectionalWeightsMarginal'
]
=
[
2.
/
(
ms
[
1
]
-
ms
[
0
])]
params
[
'matchingWeight'
]
=
1.
params
[
'samplerPivot'
]
=
QS
(
srange
,
"UNIFORM"
,
parameterMap
)
params
[
'samplerMarginal'
]
=
QS
(
mrange
,
"UNIFORM"
)
approx
=
algo
([
0
],
solver
,
mu0
=
k0
,
approx_state
=
True
,
approxParameters
=
params
,
verbosity
=
5
,
storeAllSamples
=
True
)
else
:
approx
=
algo
(
solver
,
mu0
=
k0
,
approx_state
=
True
,
approxParameters
=
params
,
verbosity
=
5
)
if
"GREEDY"
in
method
:
approx
.
setupApprox
(
"LAST"
)
else
:
approx
.
setupApprox
()
approxNN
=
NN
(
solver
,
mu0
=
k0
,
approx_state
=
True
,
verbosity
=
5
,
approxParameters
=
{
'S'
:
len
(
approx
.
mus
),
'POD'
:
params
[
'POD'
],
'sampler'
:
ES
()})
if
SMarginal
>
1
:
approxNN
.
setSamples
(
approx
.
storedSamplesFilenames
)
approx
.
purgeStoredSamples
()
for
m
in
approx
.
musMarginal
:
bode
(
freq
,
m
[
0
],
[
approx
.
getHF
,
approx
.
getApprox
,
approxNN
.
getApprox
])
else
:
approxNN
.
setSamples
(
approx
.
samplingEngine
)
bode
(
freq
,
mu
,
[
approx
.
getHF
,
approx
.
getApprox
,
approxNN
.
getApprox
])
if
SMarginal
>
1
:
bode
(
freq
,
[
1.5
*
ms
[
1
]],
[
approx
.
getHF
,
approx
.
getApprox
,
approxNN
.
getApprox
])
bode
(
freq
,
[
2.
*
ms
[
1
]],
[
approx
.
getHF
,
approx
.
getApprox
,
approxNN
.
getApprox
])
verb
=
approx
.
verbosity
approx
.
verbosity
=
0
mspace
=
np
.
linspace
(
ms
[
0
],
ms
[
-
1
],
10
)
for
j
,
t
in
enumerate
(
mspace
):
pls
=
approx
.
getPoles
([
None
,
t
])
if
j
==
0
:
poles
=
np
.
empty
((
len
(
mspace
),
len
(
pls
)),
dtype
=
np
.
complex
)
poles
[
j
]
=
pls
for
j
,
t
in
enumerate
(
approx
.
musMarginal
):
pls
=
approx
.
getPoles
([
None
,
t
[
0
][
0
]])
if
j
==
0
:
polesE
=
np
.
empty
((
SMarginal
,
len
(
pls
)),
dtype
=
np
.
complex
)
polesE
[
j
]
=
pls
approx
.
verbosity
=
verb
fig
=
plt
.
figure
(
figsize
=
(
10
,
6
))
ax
=
fig
.
add_subplot
(
1
,
1
,
1
)
ax
.
plot
(
np
.
real
(
poles
),
np
.
imag
(
poles
),
'--'
)
ax
.
plot
(
np
.
real
(
polesE
),
np
.
imag
(
polesE
),
'ko'
,
markersize
=
4
)
ax
.
set_xlabel
(
'Real'
)
ax
.
set_ylabel
(
'Imag'
)
ax
.
grid
()
plt
.
show
()
else
:
poles
=
approx
.
getPoles
()
print
(
"Poles:
\n
{}"
.
format
(
poles
))
print
(
"
\n
"
)
Event Timeline
Log In to Comment