Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85271935
mhd.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
Fri, Sep 27, 22:32
Size
3 KB
Mime Type
text/x-python
Expires
Sun, Sep 29, 22:32 (2 d)
Engine
blob
Format
Raw Data
Handle
21152689
Attached To
R6746 RationalROMPy
mhd.py
View Options
import
numpy
as
np
import
matplotlib.pyplot
as
plt
from
mhd_engine
import
MHDEngine
as
engine
from
rrompy.reduction_methods
import
(
RationalInterpolant
as
RI
,
RationalInterpolantGreedy
as
RIG
)
from
rrompy.parameter.parameter_sampling
import
(
FFTSampler
as
FFTS
,
QuadratureCircleSampler
as
QCS
,
QuadratureBoxSampler
as
QBS
)
ks
=
[
-.
35
+
.
5j
,
0.
+
.
5j
]
k0
=
np
.
mean
(
ks
)
solver
=
engine
(
5
)
kEffDelta
=
.
1
*
(
ks
[
1
]
-
ks
[
0
])
kEff
=
np
.
real
([
ks
[
0
]
-
kEffDelta
,
ks
[
1
]
+
kEffDelta
])
iEff
=
kEff
-
.
5
*
np
.
sum
(
np
.
real
(
ks
))
+
np
.
imag
(
ks
[
0
])
nPoles
=
50
polesEx
=
solver
.
getPolesExact
(
nPoles
,
k0
)
for
corrector
in
[
False
,
True
]:
for
method
in
[
"FFT"
,
"BOX"
,
"GREEDY"
]:
print
(
"Testing {} method with{} corrector"
.
format
(
method
,
"out"
*
(
not
corrector
)))
if
method
==
"FFT"
:
params
=
{
'S'
:
64
,
'POD'
:
True
,
'polybasis'
:
"MONOMIAL"
,
'sampler'
:
FFTS
(
ks
),
'correctorTol'
:
1e-5
}
algo
=
RI
if
method
==
"BOX"
:
params
=
{
'S'
:
64
,
'POD'
:
True
,
'polybasis'
:
"MONOMIAL"
,
'sampler'
:
QBS
(
ks
),
'correctorTol'
:
1e-5
}
algo
=
RI
if
method
==
"GREEDY"
:
params
=
{
'S'
:
30
,
'POD'
:
True
,
'greedyTol'
:
1e-2
,
'polybasis'
:
"MONOMIAL"
,
'sampler'
:
QCS
(
ks
),
'errorEstimatorKind'
:
"LOOK_AHEAD"
,
'nTestPoints'
:
10000
,
'trainSetGenerator'
:
FFTS
(
ks
),
'correctorTol'
:
1e-5
}
algo
=
RIG
params
[
'correctorForce'
]
=
corrector
approx
=
algo
(
solver
,
mu0
=
k0
,
approx_state
=
True
,
approxParameters
=
params
,
verbosity
=
10
)
approx
.
setupApprox
()
poles
,
residues
=
approx
.
getResidues
()
inRange
=
np
.
logical_and
(
np
.
logical_and
(
np
.
real
(
poles
)
>=
kEff
[
0
],
np
.
real
(
poles
)
<=
kEff
[
1
]),
np
.
logical_and
(
np
.
imag
(
poles
)
>=
iEff
[
0
],
np
.
imag
(
poles
)
<=
iEff
[
1
]))
polesEff
=
poles
[
inRange
]
resNormEff
=
np
.
linalg
.
norm
(
residues
,
axis
=
1
)[
inRange
]
rLm
=
np
.
min
(
np
.
log
(
resNormEff
))
rLmM
=
np
.
max
(
np
.
log
(
resNormEff
))
-
rLm
fig
=
plt
.
figure
(
figsize
=
(
10
,
10
))
ax
=
fig
.
add_subplot
(
1
,
1
,
1
)
if
method
==
"GREEDY"
:
ax
.
plot
(
approx
.
muTest
.
re
.
data
.
flatten
(),
approx
.
muTest
.
im
.
data
.
flatten
(),
'k,'
,
alpha
=
0.25
)
for
pl
,
rN
in
zip
(
polesEff
,
resNormEff
):
if
corrector
:
alpha
=
.
35
+
.
4
*
(
np
.
log
(
rN
)
-
rLm
)
/
rLmM
else
:
alpha
=
.
1
+
.
65
*
(
np
.
log
(
rN
)
-
rLm
)
/
rLmM
ax
.
annotate
(
"{0:.0e}"
.
format
(
rN
),
(
np
.
real
(
pl
),
np
.
imag
(
pl
)),
alpha
=
alpha
)
ax
.
plot
(
np
.
real
(
pl
),
np
.
imag
(
pl
),
'r+'
,
alpha
=
alpha
+
.
25
)
ax
.
plot
(
approx
.
mus
.
re
.
data
.
flatten
(),
approx
.
mus
.
im
.
data
.
flatten
(),
'k.'
)
ax
.
plot
(
np
.
real
(
polesEx
),
np
.
imag
(
polesEx
),
'bx'
)
ax
.
set_xlim
(
kEff
)
ax
.
set_ylim
(
iEff
)
ax
.
grid
()
plt
.
tight_layout
()
plt
.
show
()
print
(
"
\n
"
)
Event Timeline
Log In to Comment