Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F61523932
pod.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, 05:11
Size
6 KB
Mime Type
text/x-python
Expires
Thu, May 9, 05:11 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17523228
Attached To
R6746 RationalROMPy
pod.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.lagrange
import
ApproximantLagrangePade
as
Pade
from
rrompy.reduction_methods.lagrange
import
ApproximantLagrangeRB
as
RB
from
rrompy.utilities.parameter_sampling
import
QuadratureSampler
as
QS
verb
=
100
sol
=
"single"
#sol = "sweep"
algo
=
"Pade"
#algo = "RB"
polyBasis
=
"LEGENDRE"
polyBasis
=
"CHEBYSHEV"
#polyBasis = "MONOMIAL"
dampingEta
=
0
*
1e4
/
2.
/
np
.
pi
ktar
=
1.e4
# [Hz]
k0s
=
np
.
array
([
2.5e2
,
1.0e4
])
k0s
=
np
.
array
([
2.5e3
,
1.5e4
])
#k0s = np.array([5.0e4, 1.0e5])
#k0s = np.array([2.0e5, 3.0e5])
k0
=
np
.
mean
(
np
.
power
(
k0s
,
2.
))
**
.
5
###
if
dampingEta
>
0
:
rescaling
=
lambda
x
:
x
rescalingInv
=
lambda
x
:
x
else
:
rescaling
=
lambda
x
:
np
.
power
(
x
,
2.
)
rescalingInv
=
lambda
x
:
np
.
power
(
x
,
.
5
)
params
=
{
'N'
:
15
,
'M'
:
14
,
'S'
:
25
,
'POD'
:
True
,
'polybasis'
:
polyBasis
}
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
:
approx
=
RB
(
solver
,
mu0
=
k0
,
approxParameters
=
params
,
verbosity
=
verb
)
approx
.
sampler
=
QS
(
k0s
,
"CHEBYSHEV"
,
rescaling
,
rescalingInv
)
approx
.
setupApprox
()
if
sol
==
"single"
:
approx
.
outParaviewTimeDomainSamples
(
filename
=
"out/outSamples{}"
.
format
(
dampingEta
),
forceNewFile
=
False
,
folders
=
True
)
nameBase
=
"{}_{}"
.
format
(
ktar
,
dampingEta
)
approx
.
outParaviewTimeDomainApprox
(
ktar
,
omega
=
2.
*
np
.
pi
*
ktar
,
filename
=
"out/outTApp{}"
.
format
(
nameBase
),
forceNewFile
=
False
,
folder
=
True
)
approx
.
outParaviewTimeDomainHF
(
ktar
,
omega
=
2.
*
np
.
pi
*
ktar
,
filename
=
"out/outTHF{}"
.
format
(
nameBase
),
forceNewFile
=
False
,
folder
=
True
)
approx
.
outParaviewTimeDomainErr
(
ktar
,
omega
=
2.
*
np
.
pi
*
ktar
,
filename
=
"out/outTErr{}"
.
format
(
nameBase
),
forceNewFile
=
False
,
folder
=
True
)
approx
.
outParaviewTimeDomainRes
(
ktar
,
omega
=
2.
*
np
.
pi
*
ktar
,
filename
=
"out/outTRes{}"
.
format
(
nameBase
),
forceNewFile
=
False
,
folder
=
True
)
appErr
,
solNorm
=
approx
.
normErr
(
ktar
),
approx
.
normHF
(
ktar
)
resNorm
,
RHSNorm
=
approx
.
normRes
(
ktar
),
approx
.
normRHS
(
ktar
)
print
((
'SolNorm:
\t
{}
\n
Err:
\t
{}
\n
ErrRel:
\t
{}'
)
.
format
(
solNorm
,
appErr
,
np
.
divide
(
appErr
,
solNorm
)))
print
((
'RHSNorm:
\t
{}
\n
Res:
\t
{}
\n
ResRel:
\t
{}'
)
.
format
(
RHSNorm
,
resNorm
,
np
.
divide
(
resNorm
,
RHSNorm
)))
print
(
'
\n
Poles:'
)
print
(
approx
.
getPoles
())
if
sol
==
"sweep"
:
k0s
=
np
.
linspace
(
k0s
[
0
],
k0s
[
1
],
100
)
kl
,
kr
=
min
(
k0s
),
max
(
k0s
)
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
)
err
=
np
.
zeros_like
(
normApp
)
for
j
in
range
(
len
(
k0s
)):
normApp
[
j
]
=
approx
.
normApprox
(
k0s
[
j
])
norm
[
j
]
=
approx
.
normHF
(
k0s
[
j
])
err
[
j
]
=
approx
.
normErr
(
k0s
[
j
])
/
norm
[
j
]
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
,
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