Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F62801869
boundary_value_problem_1D.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
Wed, May 15, 18:25
Size
4 KB
Mime Type
text/x-python
Expires
Fri, May 17, 18:25 (2 d)
Engine
blob
Format
Raw Data
Handle
17690956
Attached To
R6746 RationalROMPy
boundary_value_problem_1D.py
View Options
import
warnings
import
numpy
as
np
import
matplotlib.pyplot
as
plt
from
boundary_value_problem_1D_engine
import
BVP1DEngine
as
engine
,
BVP1DPoles
from
rrompy.reduction_methods
import
(
RationalInterpolant
as
RI
,
RationalInterpolantGreedy
as
RIG
)
from
rrompy.parameter.parameter_sampling
import
(
QuadratureBoxSampler
as
QBS
,
QuadratureSampler
as
QS
)
from
rrompy.utilities.numerical
import
potential
ks
,
shift
=
[
0.
,
5e3
],
200.j
ksRatio
=
2.
*
np
.
abs
(
shift
)
/
np
.
abs
(
ks
[
1
]
-
ks
[
0
])
k0
,
gamma
,
n
=
np
.
mean
(
ks
),
.
1
,
10000
solver
=
engine
(
gamma
,
n
,
shift
)
kEffMult
=
.
1
kEff
=
np
.
real
([
ks
[
0
]
-
kEffMult
*
(
ks
[
1
]
-
ks
[
0
]),
ks
[
1
]
+
kEffMult
*
(
ks
[
1
]
-
ks
[
0
])])
methods
=
[
"25"
,
"30"
,
"35"
,
"10+"
,
"15+"
]
for
kind
in
[
"SEGMENT"
,
"CLOUD"
]:
print
(
"Testing sampling kind {}..."
.
format
(
kind
))
errAbs
,
errRel
=
[],
[]
for
method
in
methods
:
print
(
"S = {}"
.
format
(
method
))
if
len
(
method
)
==
2
:
if
kind
==
"SEGMENT"
:
sampler
=
QS
(
ks
,
"CHEBYSHEV"
)
else
:
#if kind == "CLOUD":
sampler
=
QBS
(
ks
,
"CHEBYSHEV"
,
ksRatio
)
params
=
{
'S'
:
int
(
method
),
'POD'
:
True
,
'polybasis'
:
"CHEBYSHEV"
,
'sampler'
:
sampler
}
algo
=
RI
if
method
[
-
1
]
==
"+"
:
if
kind
==
"SEGMENT"
:
sampler
=
QS
(
ks
,
"UNIFORM"
)
samplerTrainSet
=
QS
(
ks
,
"CHEBYSHEV"
)
else
:
#if kind == "CLOUD":
sampler
=
QBS
(
ks
,
"UNIFORM"
,
ksRatio
)
samplerTrainSet
=
QBS
(
ks
,
"CHEBYSHEV"
,
ksRatio
)
params
=
{
'S'
:
int
(
method
[:
-
1
]),
'POD'
:
True
,
'polybasis'
:
"LEGENDRE"
,
'greedyTol'
:
1e-3
,
'sampler'
:
sampler
,
'errorEstimatorKind'
:
"LOOK_AHEAD"
,
'nTestPoints'
:
1000
,
'samplerTrainSet'
:
samplerTrainSet
}
algo
=
RIG
approx
=
algo
(
solver
,
mu0
=
k0
,
approxParameters
=
params
,
verbosity
=
0
)
approx
.
setupApprox
()
poles
=
approx
.
getPoles
()
polesEff
=
poles
[(
np
.
real
(
poles
)
>=
kEff
[
0
])
*
(
np
.
real
(
poles
)
<=
kEff
[
1
])]
polesEx
=
BVP1DPoles
(
gamma
,
kEff
[
0
],
kEff
[
1
],
shift
)
polesExTop
=
np
.
sort
(
polesEx
[(
np
.
real
(
polesEx
)
>=
np
.
real
(
ks
[
0
]))
*
(
np
.
real
(
polesEx
)
<=
np
.
real
(
ks
[
1
]))
*
(
np
.
imag
(
polesEx
)
>
-
np
.
imag
(
shift
))])
polesExBot
=
np
.
sort
(
polesEx
[(
np
.
real
(
polesEx
)
>=
np
.
real
(
ks
[
0
]))
*
(
np
.
real
(
polesEx
)
<=
np
.
real
(
ks
[
1
]))
*
(
np
.
imag
(
polesEx
)
<
-
np
.
imag
(
shift
))])
polesEx
=
np
.
append
(
np
.
append
(
polesExTop
,
np
.
inf
),
polesExBot
)
polesExErr
=
np
.
abs
(
np
.
tile
(
polesEx
.
reshape
(
-
1
,
1
),
[
1
,
len
(
poles
)])
-
poles
.
reshape
(
1
,
-
1
))
fig
=
plt
.
figure
()
ax
=
fig
.
add_subplot
(
1
,
1
,
1
)
if
method
[
-
1
]
==
"+"
:
print
(
"S_eff = {}"
.
format
(
approx
.
S
))
ax
.
plot
(
approx
.
muTest
.
re
.
data
.
flatten
(),
approx
.
muTest
.
im
.
data
.
flatten
(),
'k,'
,
alpha
=
0.25
)
ax
.
plot
(
approx
.
mus
.
re
.
data
.
flatten
(),
approx
.
mus
.
im
.
data
.
flatten
(),
'k.'
)
ax
.
plot
(
np
.
real
(
polesEff
),
np
.
imag
(
polesEff
),
'r+'
)
ax
.
plot
(
np
.
real
(
polesEx
),
np
.
imag
(
polesEx
),
'bx'
)
ax
.
set_xlim
(
kEff
)
ax
.
grid
()
plt
.
tight_layout
()
plt
.
show
()
errAbs
+=
[
np
.
min
(
polesExErr
,
axis
=
1
)]
with
warnings
.
catch_warnings
():
warnings
.
simplefilter
(
"ignore"
)
potEx
=
potential
(
approx
.
trainedModel
.
centerNormalize
(
polesEx
)(
0
),
sampler
.
normalFoci
())
potEx
[
potEx
<
1.
]
=
1.
potEx
[
len
(
potEx
)
//
2
]
=
np
.
nan
errRel
+=
[
2.
/
np
.
abs
(
ks
[
1
]
-
ks
[
0
])
*
errAbs
[
-
1
]
/
potEx
]
symbols
=
'.x+dso^v<>'
fig
=
plt
.
figure
(
figsize
=
plt
.
figaspect
(
.
5
))
ax1
=
fig
.
add_subplot
(
1
,
2
,
1
)
ax2
=
fig
.
add_subplot
(
1
,
2
,
2
)
for
j
in
range
(
len
(
methods
)):
ax1
.
semilogy
(
np
.
real
(
polesEx
),
errAbs
[
j
],
'{}-'
.
format
(
symbols
[
j
]))
ax2
.
semilogy
(
np
.
real
(
polesEx
),
errRel
[
j
],
'{}-'
.
format
(
symbols
[
j
]))
ax1
.
set_title
(
"Absolute"
)
ax1
.
grid
()
ax2
.
set_title
(
"Relative"
)
ax2
.
legend
([
"S = {}"
.
format
(
m
)
for
m
in
methods
])
ax2
.
grid
()
plt
.
tight_layout
()
plt
.
show
()
print
(
"
\n
"
)
Event Timeline
Log In to Comment