Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61142439
greedy.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
Sat, May 4, 19:38
Size
5 KB
Mime Type
text/x-python
Expires
Mon, May 6, 19:38 (2 d)
Engine
blob
Format
Raw Data
Handle
17471396
Attached To
R6746 RationalROMPy
greedy.py
View Options
import
numpy
as
np
import
fenics
as
fen
import
ufl
from
rrompy.hfengines.vector_linear_problem
import
\
LinearElasticityHelmholtzProblemEngine
as
LEHPE
from
rrompy.hfengines.vector_linear_problem
import
\
LinearElasticityHelmholtzProblemEngineDamped
as
LEHPED
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
verb
=
2
timed
=
False
algo
=
"Pade"
#algo = "RB"
polyBasis
=
"LEGENDRE"
polyBasis
=
"CHEBYSHEV"
#polyBasis = "MONOMIAL"
if
timed
:
verb
=
0
dampingEta
=
0
*
1e4
/
2.
/
np
.
pi
k0s
=
np
.
linspace
(
2.5e2
,
7.5e3
,
100
)
k0s
=
np
.
linspace
(
2.5e3
,
1.5e4
,
100
)
k0s
=
np
.
linspace
(
5.0e4
,
1.0e5
,
100
)
k0s
=
np
.
linspace
(
2.0e5
,
2.5e5
,
100
)
k0
=
np
.
mean
(
np
.
power
(
k0s
,
2.
))
**
.
5
# [Hz]
kl
,
kr
=
min
(
k0s
),
max
(
k0s
)
params
=
{
'muBounds'
:[
kl
,
kr
],
'nTestPoints'
:
500
,
'Delta'
:
0
,
'greedyTol'
:
1e-2
,
'S'
:
2
,
'polybasis'
:
polyBasis
,
'robustTol'
:
2e-16
,
'interpRcond'
:
None
,
'errorEstimatorKind'
:
'EXACT'
}
theta
=
20.
*
np
.
pi
/
180.
phi
=
10.
*
np
.
pi
/
180.
mesh
=
fen
.
Mesh
(
"../data/mesh/diapason_1.xml"
)
subdomains
=
fen
.
MeshFunction
(
"size_t"
,
mesh
,
"../data/mesh/diapason_1_physical_region.xml"
)
meshBall
=
fen
.
SubMesh
(
mesh
,
subdomains
,
2
)
meshFork
=
fen
.
SubMesh
(
mesh
,
subdomains
,
1
)
Hball
=
np
.
max
(
meshBall
.
coordinates
()[:,
1
])
#.00257
Ltot
=
np
.
max
(
mesh
.
coordinates
()[:,
1
])
#.1022
Lhandle
=
np
.
max
(
meshFork
.
coordinates
()[:,
1
])
#.026
Lrod
=
Ltot
-
Lhandle
#.0762
L2var
=
(
Lrod
/
4.
)
**
2.
Ehandle_ratio
=
3.
rhohandle_ratio
=
1.5
c
=
3.e2
rho
=
8e3
*
(
2.
*
np
.
pi
)
**
2.
E
=
1.93e11
nu
=
.
3
T
=
1e6
lambda_
=
E
*
nu
/
(
1.
+
nu
)
/
(
1.
-
2.
*
nu
)
mu_
=
E
/
(
1.
+
nu
)
kWave
=
(
np
.
cos
(
theta
)
*
np
.
cos
(
phi
),
np
.
sin
(
phi
),
np
.
sin
(
theta
)
*
np
.
cos
(
phi
))
x
,
y
,
z
=
fen
.
SpatialCoordinate
(
mesh
)[:]
yCorr
=
y
-
Ltot
compPlane
=
kWave
[
0
]
*
x
+
kWave
[
1
]
*
y
+
kWave
[
2
]
*
z
xPlane
,
yPlane
,
zPlane
=
(
xx
-
compPlane
*
xx
for
xx
in
(
x
,
y
,
z
))
xOrtho
,
yOrtho
,
zOrtho
=
(
compPlane
*
xx
for
xx
in
(
x
,
y
,
z
))
forcingBase
=
(
T
/
(
2.
*
np
.
pi
*
L2var
)
**.
5
*
fen
.
exp
(
-
(
xPlane
**
2.
+
yPlane
**
2.
+
zPlane
**
2.
)
/
(
2.
*
L2var
)))
forcingWeight
=
np
.
real
(
k0
)
/
c
*
(
xOrtho
+
yOrtho
+
zOrtho
)
neumannDatum
=
[
ufl
.
as_vector
(
tuple
(
forcingBase
*
fen
.
cos
(
forcingWeight
)
*
kWavedir
for
kWavedir
in
kWave
)),
ufl
.
as_vector
(
tuple
(
forcingBase
*
fen
.
sin
(
forcingWeight
)
*
kWavedir
for
kWavedir
in
kWave
))]
lambda_eff
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
Lhandle
),
fen
.
Constant
(
lambda_
),
fen
.
Constant
(
Ehandle_ratio
*
lambda_
))
mu_eff
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
Lhandle
),
fen
.
Constant
(
mu_
),
fen
.
Constant
(
Ehandle_ratio
*
mu_
))
rho_eff
=
ufl
.
conditional
(
ufl
.
ge
(
y
,
Lhandle
),
fen
.
Constant
(
rho
),
fen
.
Constant
(
rhohandle_ratio
*
rho
))
###
if
dampingEta
>
0
:
solver
=
LEHPED
(
degree_threshold
=
8
,
verbosity
=
0
)
solver
.
eta
=
dampingEta
else
:
solver
=
LEHPE
(
degree_threshold
=
8
,
verbosity
=
0
)
solver
.
omega
=
np
.
real
(
k0
)
solver
.
lambda_
=
lambda_eff
solver
.
mu_
=
mu_eff
solver
.
rho_
=
rho_eff
solver
.
V
=
fen
.
VectorFunctionSpace
(
mesh
,
"P"
,
1
)
solver
.
DirichletBoundary
=
lambda
x
,
on_b
:
on_b
and
x
[
1
]
<
Hball
solver
.
NeumannBoundary
=
"REST"
solver
.
forcingTerm
=
neumannDatum
if
algo
==
"Pade"
:
approx
=
Pade
(
solver
,
mu0
=
k0
,
approxParameters
=
params
,
verbosity
=
verb
)
else
:
params
.
pop
(
"Delta"
)
params
.
pop
(
"polybasis"
)
params
.
pop
(
"robustTol"
)
params
.
pop
(
"interpRcond"
)
params
.
pop
(
"errorEstimatorKind"
)
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
)
polesApp
=
approx
.
getPoles
()
print
(
"Poles:
\n
"
,
polesApp
)
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
])
/
norm
[
j
]
resApp
=
approx
.
errorEstimator
(
k0s
)
res
[
res
<
1e-5
*
approx
.
greedyTol
]
=
np
.
nan
resApp
[
resApp
<
1e-5
*
approx
.
greedyTol
]
=
np
.
nan
err
[
err
<
1e-8
*
approx
.
greedyTol
]
=
np
.
nan
plt
.
figure
()
plt
.
semilogy
(
k0s
,
norm
)
plt
.
semilogy
(
k0s
,
normApp
,
'--'
)
plt
.
semilogy
(
np
.
real
(
approx
.
mus
),
1.05
*
np
.
max
(
norm
)
*
np
.
ones_like
(
approx
.
mus
,
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
),
approx
.
greedyTol
*
np
.
ones_like
(
approx
.
mus
,
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
()
mask
=
(
np
.
real
(
polesApp
)
<
kl
)
|
(
np
.
real
(
polesApp
)
>
kr
)
print
(
"Outliers:"
,
polesApp
[
mask
])
polesAppEff
=
polesApp
[
~
mask
]
plt
.
figure
()
plt
.
plot
(
np
.
real
(
polesAppEff
),
np
.
imag
(
polesAppEff
),
'kx'
)
plt
.
axis
(
'equal'
)
plt
.
grid
()
plt
.
show
()
plt
.
close
()
Event Timeline
Log In to Comment