Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F63588930
greedy_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
Tue, May 21, 04:52
Size
4 KB
Mime Type
text/x-python
Expires
Thu, May 23, 04:52 (2 d)
Engine
blob
Format
Raw Data
Handle
17792031
Attached To
R6746 RationalROMPy
greedy_scatteringAirfoil.py
View Options
import
numpy
as
np
import
fenics
as
fen
import
ufl
from
helmholtz_box_scattering_problem_engine
import
\
HelmholtzBoxScatteringProblemEngine
as
HSP
from
rrompy.reduction_methods.greedy
import
RationalInterpolantGreedy
as
RI
from
rrompy.reduction_methods.greedy
import
ReducedBasisGreedy
as
RB
from
rrompy.solver.fenics
import
fenONE
,
L2NormMatrix
from
rrompy.parameter.parameter_sampling
import
QuadratureSampler
as
QS
verb
=
2
timed
=
False
algo
=
"RI"
#algo = "RB"
polyBasis
=
"LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
k0s
=
np
.
linspace
(
5
,
20
,
25
)
k0
=
np
.
mean
(
k0s
)
kl
,
kr
=
min
(
k0s
),
max
(
k0s
)
params
=
{
'sampler'
:
QS
([
kl
,
kr
],
"UNIFORM"
),
'nTestPoints'
:
500
,
'greedyTol'
:
1e-2
,
'S'
:
10
,
'polybasis'
:
polyBasis
,
'errorEstimatorKind'
:
'DIAGONAL'
}
# 'errorEstimatorKind':'INTERPOLATORY'}
# 'errorEstimatorKind':'AFFINE'}
#########
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
.
real
(
k0
),
theta
,
n
=
1
,
verbosity
=
verb
,
degree_threshold
=
8
)
solver
.
omega
=
np
.
real
(
k0
)
solver
.
V
=
fen
.
FunctionSpace
(
mesh
,
"P"
,
2
)
solver
.
diffusivity
=
a
solver
.
DirichletBoundary
=
Dboundary
solver
.
RobinBoundary
=
"REST"
solver
.
DirichletDatum
=
[
u0R
,
u0I
]
#########
if
algo
==
"RI"
:
approx
=
RI
(
solver
,
mu0
=
k0
,
approxParameters
=
params
,
verbosity
=
verb
)
else
:
params
.
pop
(
"polybasis"
)
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
)
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