Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F71858977
pod_scatteringAirfoil.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, Jul 13, 11:13
Size
5 KB
Mime Type
text/x-python
Expires
Mon, Jul 15, 11:13 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
19008704
Attached To
R6746 RationalROMPy
pod_scatteringAirfoil.py
View Options
import
numpy
as
np
import
fenics
as
fen
import
ufl
from
rrompy.hfengines.linear_problem
import
\
HelmholtzBoxScatteringProblemEngine
as
HSP
from
rrompy.reduction_methods.centered
import
RationalPade
as
PC
from
rrompy.reduction_methods.centered
import
RBCentered
as
RBC
from
rrompy.reduction_methods.distributed
import
RationalInterpolant
as
PD
from
rrompy.reduction_methods.distributed
import
RBDistributed
as
RBD
from
rrompy.utilities.parameter_sampling
import
QuadratureSampler
as
QS
from
rrompy.solver.fenics
import
fenONE
from
operator
import
itemgetter
def
subdict
(
d
,
ks
):
return
dict
(
zip
(
ks
,
itemgetter
(
*
ks
)(
d
)))
verb
=
0
####################
homog
=
True
#homog = False
####################
test
=
"solve"
test
=
"Centered"
test
=
"Distributed"
plotSamples
=
True
k0
=
10
kLeft
,
kRight
=
8
+
0.j
,
12
+
0.j
ktar
=
11
ktars
=
np
.
linspace
(
8
,
12
,
21
)
+
0.j
PI
=
np
.
pi
R
=
2
def
Dboundary
(
x
,
on_boundary
):
return
on_boundary
and
(
x
[
0
]
**
2
+
x
[
1
]
**
2
)
**.
5
<
.
95
*
R
kappa
=
10
theta
=
PI
*
-
45
/
180.
mu
=
1.1
epsilon
=
.
1
mesh
=
fen
.
Mesh
(
'../data/mesh/airfoil.xml'
)
c
,
s
=
np
.
cos
(
theta
),
np
.
sin
(
theta
)
x
,
y
=
fen
.
SpatialCoordinate
(
mesh
)[:]
u0R
=
-
fen
.
cos
(
kappa
*
(
c
*
x
+
s
*
y
))
u0I
=
-
fen
.
sin
(
kappa
*
(
c
*
x
+
s
*
y
))
checkReal
=
x
**
2
-
x
+
y
**
2
rhop5
=
((
x
**
2
+
y
**
2
)
/
((
x
-
1
)
**
2
+
y
**
2
))
**.
25
phiroot1
=
fen
.
atan
(
-
y
/
(
x
**
2
-
x
+
y
**
2
))
/
2
phiroot2
=
fen
.
atan
(
-
y
/
(
x
**
2
-
x
+
y
**
2
))
/
2
-
PI
*
ufl
.
sign
(
-
y
/
(
x
**
2
-
x
+
y
**
2
))
/
2
kappam1
=
(((
rhop5
*
fen
.
cos
(
phiroot1
)
+.
5
)
**
2.
+
(
rhop5
*
fen
.
sin
(
phiroot1
))
**
2.
)
/
((
rhop5
*
fen
.
cos
(
phiroot1
)
-
1.
)
**
2.
+
(
rhop5
*
fen
.
sin
(
phiroot1
))
**
2.
)
)
**.
5
-
mu
kappam2
=
(((
rhop5
*
fen
.
cos
(
phiroot2
)
+.
5
)
**
2.
+
(
rhop5
*
fen
.
sin
(
phiroot2
))
**
2.
)
/
((
rhop5
*
fen
.
cos
(
phiroot2
)
-
1.
)
**
2.
+
(
rhop5
*
fen
.
sin
(
phiroot2
))
**
2.
)
)
**.
5
-
mu
Heps1
=
.
9
*
.
5
*
(
1
+
kappam1
/
epsilon
+
fen
.
sin
(
PI
*
kappam1
/
epsilon
)
/
PI
)
+
.
1
Heps2
=
.
9
*
.
5
*
(
1
+
kappam2
/
epsilon
+
fen
.
sin
(
PI
*
kappam2
/
epsilon
)
/
PI
)
+
.
1
cTT
=
ufl
.
conditional
(
ufl
.
le
(
kappam1
,
epsilon
),
Heps1
,
fenONE
)
c_F
=
fen
.
Constant
(
.
1
)
cFT
=
ufl
.
conditional
(
ufl
.
le
(
kappam2
,
epsilon
),
Heps2
,
fenONE
)
c_F
=
fen
.
Constant
(
.
1
)
cT
=
ufl
.
conditional
(
ufl
.
ge
(
kappam1
,
-
epsilon
),
cTT
,
c_F
)
cF
=
ufl
.
conditional
(
ufl
.
ge
(
kappam2
,
-
epsilon
),
cFT
,
c_F
)
a
=
ufl
.
conditional
(
ufl
.
ge
(
checkReal
,
0.
),
cT
,
cF
)
###
solver
=
HSP
(
R
,
np
.
abs
(
k0
),
theta
,
n
=
1
,
verbosity
=
verb
,
degree_threshold
=
8
)
solver
.
V
=
fen
.
FunctionSpace
(
mesh
,
"P"
,
3
)
solver
.
diffusivity
=
a
solver
.
DirichletBoundary
=
Dboundary
solver
.
RobinBoundary
=
"REST"
solver
.
DirichletDatum
=
[
u0R
,
u0I
]
###
if
test
==
"solve"
:
uinc
=
solver
.
liftDirichletData
(
k0
)
if
homog
:
uhtot
=
solver
.
solve
(
k0
,
homogeneized
=
homog
)
uh
=
uhtot
+
uinc
else
:
uh
=
solver
.
solve
(
k0
,
homogeneized
=
homog
)
uhtot
=
uh
-
uinc
print
(
solver
.
norm
(
uh
))
print
(
solver
.
norm
(
uhtot
))
solver
.
plot
(
fen
.
project
(
a
,
solver
.
V
)
.
vector
(),
what
=
'Real'
,
name
=
'a'
)
solver
.
plot
(
uinc
,
what
=
'Real'
,
name
=
'u_inc'
)
solver
.
plot
(
uh
,
what
=
'ABS'
)
solver
.
plot
(
uhtot
,
what
=
'ABS'
,
name
=
'u + u_inc'
)
elif
test
in
[
"Centered"
,
"Distributed"
]:
if
test
==
"Centered"
:
params
=
{
'N'
:
8
,
'M'
:
8
,
'R'
:
8
,
'E'
:
8
,
'sampleType'
:
'Arnoldi'
,
'POD'
:
True
}
parPade
=
subdict
(
params
,
[
'N'
,
'M'
,
'E'
,
'sampleType'
,
'POD'
])
parRB
=
subdict
(
params
,
[
'R'
,
'E'
,
'sampleType'
,
'POD'
])
approxPade
=
PC
(
solver
,
mu0
=
k0
,
approxParameters
=
parPade
,
verbosity
=
verb
,
homogeneized
=
homog
)
approxRB
=
RBC
(
solver
,
mu0
=
k0
,
approxParameters
=
parRB
,
verbosity
=
verb
,
homogeneized
=
homog
)
else
:
params
=
{
'N'
:
8
,
'M'
:
8
,
'R'
:
9
,
'S'
:
9
,
'POD'
:
True
,
'basis'
:
"CHEBYSHEV"
,
'sampler'
:
QS
([
kLeft
,
kRight
],
"CHEBYSHEV"
)}
parPade
=
subdict
(
params
,
[
'N'
,
'M'
,
'S'
,
'POD'
,
'basis'
,
'sampler'
])
parRB
=
subdict
(
params
,
[
'R'
,
'S'
,
'POD'
,
'sampler'
])
approxPade
=
PD
(
solver
,
mu0
=
np
.
mean
([
kLeft
,
kRight
]),
approxParameters
=
parPade
,
verbosity
=
verb
,
homogeneized
=
homog
)
approxRB
=
RBD
(
solver
,
mu0
=
np
.
mean
([
kLeft
,
kRight
]),
approxParameters
=
parRB
,
verbosity
=
verb
,
homogeneized
=
homog
)
approxPade
.
setupApprox
()
approxRB
.
setupApprox
()
if
plotSamples
:
approxPade
.
plotSamples
()
approxPade
.
plotHF
(
ktar
,
name
=
'u_HF'
)
approxPade
.
plotApprox
(
ktar
,
name
=
'u_Pade'''
)
approxPade
.
plotErr
(
ktar
,
name
=
'err_Pade'''
)
approxPade
.
plotRes
(
ktar
,
name
=
'res_Pade'''
)
approxRB
.
plotApprox
(
ktar
,
name
=
'u_RB'
)
approxRB
.
plotErr
(
ktar
,
name
=
'err_RB'
)
approxRB
.
plotRes
(
ktar
,
name
=
'res_RB'
)
HFNorm
,
RHSNorm
=
approxPade
.
normHF
(
ktar
),
approxPade
.
normRHS
(
ktar
)
PadeRes
,
PadeErr
=
approxPade
.
normRes
(
ktar
),
approxPade
.
normErr
(
ktar
)
RBRes
,
RBErr
=
approxRB
.
normRes
(
ktar
),
approxRB
.
normErr
(
ktar
)
print
(
'HFNorm:
\t
{}
\n
RHSNorm:
\t
{}'
.
format
(
HFNorm
,
RHSNorm
))
print
(
'PadeRes:
\t
{}
\n
PadeErr:
\t
{}'
.
format
(
PadeRes
,
PadeErr
))
print
(
'RBRes:
\t
{}
\n
RBErr:
\t
{}'
.
format
(
RBRes
,
RBErr
))
print
(
'
\n
Poles Pade'':'
)
print
(
approxPade
.
getPoles
())
Event Timeline
Log In to Comment