Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88668625
greedy_internalBox.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, Oct 20, 02:01
Size
3 KB
Mime Type
text/x-python
Expires
Tue, Oct 22, 02:01 (2 d)
Engine
blob
Format
Raw Data
Handle
21805327
Attached To
R6746 RationalROMPy
greedy_internalBox.py
View Options
import
numpy
as
np
import
fenics
as
fen
from
rrompy.hfengines.linear_problem
import
HelmholtzProblemEngine
as
HPE
from
rrompy.reduction_methods.distributed_greedy
import
\
RationalInterpolantGreedy
as
Pade
from
rrompy.reduction_methods.distributed_greedy
import
\
RBDistributedGreedy
as
RB
from
rrompy.solver.fenics
import
L2NormMatrix
dim
=
3
verb
=
2
timed
=
False
algo
=
"Pade"
algo
=
"RB"
polyBasis
=
"LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
k0s
=
np
.
power
(
np
.
linspace
(
500
**
2.
,
2250
**
2.
,
200
),
.
5
)
k0
=
np
.
mean
(
np
.
power
(
k0s
,
2.
))
**
.
5
kl
,
kr
=
min
(
k0s
),
max
(
k0s
)
params
=
{
'muBounds'
:[
kl
,
kr
],
'nTestPoints'
:
500
,
'Delta'
:
0
,
'greedyTol'
:
1e-2
,
'S'
:
2
,
'basis'
:
polyBasis
}
if
dim
==
2
:
mesh
=
fen
.
RectangleMesh
(
fen
.
Point
(
0.
,
0.
),
fen
.
Point
(
.
1
,
.
15
),
10
,
15
)
x
,
y
=
fen
.
SpatialCoordinate
(
mesh
)[:]
f
=
fen
.
exp
(
-
1e2
*
(
x
+
y
))
else
:
#if dim == 3:
mesh
=
fen
.
BoxMesh
(
fen
.
Point
(
0.
,
0.
,
0.
),
fen
.
Point
(
.
1
,
.
15
,
.
25
),
4
,
6
,
10
)
x
,
y
,
z
=
fen
.
SpatialCoordinate
(
mesh
)[:]
f
=
fen
.
exp
(
-
1e2
*
(
x
+
y
+
z
))
solver
=
HPE
(
np
.
real
(
k0
),
verbosity
=
verb
)
solver
.
V
=
fen
.
FunctionSpace
(
mesh
,
"P"
,
3
)
solver
.
refractionIndex
=
fen
.
Constant
(
1.
/
54.6
)
solver
.
forcingTerm
=
f
solver
.
NeumannBoundary
=
"ALL"
#########
if
algo
==
"Pade"
:
approx
=
Pade
(
solver
,
mu0
=
k0
,
approxParameters
=
params
,
verbosity
=
verb
)
else
:
approx
=
RB
(
solver
,
mu0
=
k0
,
approxParameters
=
params
,
verbosity
=
verb
)
approx
.
initEstimatorNormEngine
(
L2NormMatrix
(
solver
.
V
))
if
timed
:
from
time
import
clock
start_time
=
clock
()
approx
.
greedy
()
print
(
"Time: "
,
clock
()
-
start_time
)
else
:
approx
.
greedy
(
True
)
approx
.
samplingEngine
.
verbosity
=
0
approx
.
verbosity
=
0
kl
,
kr
=
np
.
real
(
kl
),
np
.
real
(
kr
)
from
matplotlib
import
pyplot
as
plt
normApp
=
np
.
zeros
(
len
(
k0s
))
norm
=
np
.
zeros_like
(
normApp
)
res
=
np
.
zeros_like
(
normApp
)
err
=
np
.
zeros_like
(
normApp
)
for
j
in
range
(
len
(
k0s
)):
normApp
[
j
]
=
approx
.
normApprox
(
k0s
[
j
])
norm
[
j
]
=
approx
.
normHF
(
k0s
[
j
])
res
[
j
]
=
(
approx
.
estimatorNormEngine
.
norm
(
approx
.
getRes
(
k0s
[
j
]))
/
approx
.
estimatorNormEngine
.
norm
(
approx
.
getRHS
(
k0s
[
j
])))
err
[
j
]
=
approx
.
normErr
(
k0s
[
j
])
/
approx
.
normHF
(
k0s
[
j
])
resApp
=
approx
.
errorEstimator
(
k0s
)
plt
.
figure
()
plt
.
plot
(
k0s
,
norm
)
plt
.
plot
(
k0s
,
normApp
,
'--'
)
plt
.
plot
(
np
.
real
(
approx
.
mus
(
0
)),
1.05
*
np
.
max
(
norm
)
*
np
.
ones
(
approx
.
mus
.
shape
,
dtype
=
float
),
'rx'
)
plt
.
xlim
([
kl
,
kr
])
plt
.
grid
()
plt
.
show
()
plt
.
close
()
plt
.
figure
()
plt
.
semilogy
(
k0s
,
res
)
plt
.
semilogy
(
k0s
,
resApp
,
'--'
)
plt
.
semilogy
(
np
.
real
(
approx
.
mus
(
0
)),
4.
*
np
.
max
(
resApp
)
*
np
.
ones
(
approx
.
mus
.
shape
,
dtype
=
float
),
'rx'
)
plt
.
xlim
([
kl
,
kr
])
plt
.
grid
()
plt
.
show
()
plt
.
close
()
plt
.
figure
()
plt
.
semilogy
(
k0s
,
err
)
plt
.
xlim
([
kl
,
kr
])
plt
.
grid
()
plt
.
show
()
plt
.
close
()
polesApp
=
approx
.
getPoles
()
mask
=
(
np
.
real
(
polesApp
)
<
kl
)
|
(
np
.
real
(
polesApp
)
>
kr
)
print
(
"Outliers:"
,
polesApp
[
mask
])
polesApp
=
polesApp
[
~
mask
]
plt
.
figure
()
plt
.
plot
(
np
.
real
(
polesApp
),
np
.
imag
(
polesApp
),
'kx'
)
plt
.
axis
(
'equal'
)
plt
.
grid
()
plt
.
show
()
plt
.
close
()
Event Timeline
Log In to Comment