Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78369337
HelmholtzAirfoilTaylor.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, Aug 20, 08:20
Size
9 KB
Mime Type
text/x-python
Expires
Thu, Aug 22, 08:20 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20034132
Attached To
R6746 RationalROMPy
HelmholtzAirfoilTaylor.py
View Options
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import
numpy
as
np
from
context
import
FreeFemHelmholtzScatteringEngine
as
HFS
from
context
import
FreeFemHelmholtzScatteringAugmentedEngine
as
HFSA
from
context
import
FreeFemHSEngine
as
HS
from
context
import
FreeFemHSAugmentedEngine
as
HSA
from
context
import
ROMApproximantTaylorPade
as
TP
from
context
import
ROMApproximantTaylorRB
as
TRB
from
context
import
ROMApproximantLagrangePade
as
LP
from
context
import
ROMApproximantLagrangeRB
as
LRB
from
context
import
ROMApproximantSweeper
as
Sweeper
from
matplotlib
import
pyplot
as
plt
from
operator
import
itemgetter
def
subdict
(
d
,
ks
):
return
dict
(
zip
(
ks
,
itemgetter
(
*
ks
)(
d
)))
####################
test
=
"solve"
test
=
"Taylor"
#test = "Lagrange"
test
=
"TaylorSweep"
#test = "LagrangeSweep"
ttype
=
"simple"
#ttype = "augmentedI"
#ttype = "augmentedM"
plotSamples
=
'ALL'
#plotSamples = []
k0
=
10
+
1.j
kLeft
,
kRight
=
8
+
1.j
,
12
+
1.j
ktar
=
11
ktars
=
np
.
linspace
(
8
,
12
,
33
)
-
.
5j
ktars
=
np
.
linspace
(
12.125
,
12.5
,
4
)
-
.
5j
####################
PI
=
np
.
pi
kappa
=
10
theta
=
PI
*
-
30
/
180.
mu
=
1.1
epsilon
=
.
1
V
=
(
"real epsilon = {}, mu = {};
\n
"
"macro sign(f) (f >= 0 ? 1. : -1.)//
\n
"
"macro mBase1() (x^2 - x + y^2)//
\n
"
"macro mBase2() (((x^2 + y^2) / ((x - 1)^2 + y^2))^.25)//
\n
"
"macro mBase3() (- y / mBase1)//
\n
"
"macro mPhi1() (atan(mBase3) / 2)//
\n
"
"macro mPhi2() (mPhi1 - pi/2 * sign(mBase3))//
\n
"
"macro mK1() ((((1 + cos(mPhi1) / mBase2 + (.5 / mBase2)^2.) / (1 - 2. * cos(mPhi1) / mBase2 + mBase2^-2.))^.5 - mu)/ epsilon)//
\n
"
"macro mK2() ((((1 + cos(mPhi2) / mBase2 + (.5 / mBase2)^2.) / (1 - 2. * cos(mPhi2) / mBase2 + mBase2^-2.))^.5 - mu)/ epsilon)//
\n
"
"macro mHep1() (.9 * .5 * (1. + mK1 + sin(mK1) / pi) + .1)//
\n
"
"macro mHep2() (.9 * .5 * (1. + mK2 + sin(mK2) / pi) + .1)//
\n
"
"int n = 100, R = 2;
\n
"
"border wing(t = 0, 2 * pi){{x = cos(t)/3 + 5./12 + (3*cos(t) - 3/4) / (17 - 8*cos(t)); y = sin(t)/3 - 3*sin(t) / (17 - 8*cos(t)); label = 1;}}
\n
"
"border Inf(t = 0, 2 * pi){{x = R * cos(t); y = R * sin(t); label = 2;}}
\n
"
"mesh Th = buildmesh(wing(-n) + Inf(n));
\n
"
"load
\"
Element_P3
\"
;
\n
"
"fespace V(Th, P3);"
)
.
format
(
epsilon
,
mu
)
u0
=
"- exp(1.i * {} * ({} * x + {} * y))"
.
format
(
kappa
,
np
.
cos
(
theta
),
np
.
sin
(
theta
))
a
=
(
"(mBase1 >= 0 ? (mK1>= -1. ? (mK1 <= 1. ? mHep1 : 1) : .1) : "
"(mK2 >= -1. ? (mK2 <= 1. ? mHep2 : 1) : .1))"
)
if
ttype
==
"simple"
:
solver
=
HFS
(
V
,
k0
,
diffusivity
=
a
,
DirichletBoundary
=
"1"
,
RobinBoundary
=
"2"
,
DirichletDatum
=
u0
)
plotter
=
HS
(
solver
.
V
)
else
:
if
ttype
[
-
1
]
==
"I"
:
constraintType
=
"IDENTITY"
else
:
constraintType
=
"MASS"
solver
=
HFSA
(
V
,
k0
,
diffusivity
=
a
,
DirichletBoundary
=
"1"
,
RobinBoundary
=
"2"
,
DirichletDatum
=
u0
,
constraintType
=
constraintType
)
plotter
=
HSA
(
solver
.
V
,
d
=
2
)
uinc
=
solver
.
liftDirichletData
()
if
ttype
[:
-
1
]
==
"augmented"
:
uinc
=
[
uinc
[:
int
(
len
(
uinc
)
/
2
)],
kappa
*
uinc
[:
int
(
len
(
uinc
)
/
2
)]]
if
test
==
"solve"
:
av
=
plotter
.
convExpression
(
a
)
uh
=
solver
.
solve
()
print
(
plotter
.
norm
(
uh
,
kappa
))
if
ttype
==
"simple"
:
uhtot
=
uh
-
uinc
else
:
uhtot
=
[
x
-
y
for
x
,
y
in
zip
(
uh
,
uinc
)]
print
(
plotter
.
norm
(
uhtot
,
kappa
))
plotter
.
plot
(
av
,
split
=
False
,
what
=
'Real'
,
name
=
'a'
)
plotter
.
plot
(
uhtot
-
uh
,
what
=
'Real'
,
name
=
'u_inc'
)
plotter
.
plot
(
uh
,
what
=
'ABS'
)
plotter
.
plot
(
uhtot
,
what
=
'ABS'
,
name
=
'u + u_inc'
)
elif
test
in
[
"Taylor"
,
"Lagrange"
]:
if
test
==
"Taylor"
:
params
=
{
'N'
:
8
,
'M'
:
7
,
'R'
:
8
,
'E'
:
8
,
'sampleType'
:
'Krylov'
,
'POD'
:
False
}
parPade
=
subdict
(
params
,
[
'N'
,
'M'
,
'E'
,
'sampleType'
,
'POD'
])
parRB
=
subdict
(
params
,
[
'R'
,
'E'
,
'sampleType'
,
'POD'
])
approxPade
=
TP
(
solver
,
plotter
,
k0
=
k0
,
plotDer
=
plotSamples
,
approxParameters
=
parPade
)
approxRB
=
TRB
(
solver
,
plotter
,
k0
=
k0
,
approxParameters
=
parRB
)
else
:
params
=
{
'N'
:
8
,
'M'
:
7
,
'R'
:
9
,
'S'
:
9
,
'polyBasis'
:
'CHEBYSHEV'
,
'nodesType'
:
'CHEBYSHEV'
,
'POD'
:
True
}
parPade
=
subdict
(
params
,
[
'N'
,
'M'
,
'S'
,
'polyBasis'
,
'POD'
])
parRB
=
subdict
(
params
,
[
'R'
,
'S'
,
'nodesType'
,
'POD'
])
approxPade
=
LP
(
solver
,
plotter
,
ks
=
[
kLeft
,
kRight
],
w
=
kappa
,
plotSnap
=
plotSamples
,
approxParameters
=
parPade
)
approxRB
=
LRB
(
solver
,
plotter
,
ks
=
[
kLeft
,
kRight
],
w
=
kappa
,
approxParameters
=
parRB
)
PadeErr
,
solNorm
=
approxPade
.
approxError
(
ktar
),
approxPade
.
HFNorm
(
ktar
)
RBErr
=
approxRB
.
approxError
(
ktar
)
print
((
'SolNorm:
\t
{}
\n
ErrPade:
\t
{}
\n
ErrRelPade:
\t
{}
\n
ErrRB:
\t\t
{}'
'
\n
ErrRelRB:
\t
{}'
)
.
format
(
solNorm
,
PadeErr
,
np
.
divide
(
PadeErr
,
solNorm
),
RBErr
,
np
.
divide
(
RBErr
,
solNorm
)))
print
(
'
\n
Poles Pade'':'
)
print
(
approxPade
.
getPoles
(
True
))
print
(
'
\n
Poles RB:'
)
print
(
approxRB
.
getPoles
(
True
))
uHF
=
approxPade
.
getHF
(
ktar
)
plotter
.
plot
(
uHF
,
name
=
'u_ex'
)
approxPade
.
plotApp
(
ktar
,
name
=
'u_Pade'''
)
approxPade
.
plotErr
(
ktar
,
name
=
'errPade'''
)
approxRB
.
plotApp
(
ktar
,
name
=
'u_RB'
)
approxRB
.
plotErr
(
ktar
,
name
=
'errRB'
)
elif
test
in
[
"TaylorSweep"
,
"LagrangeSweep"
]:
if
test
==
"TaylorSweep"
:
shift
=
2
nsets
=
4
stride
=
3
Emax
=
stride
*
(
nsets
-
1
)
+
shift
+
1
params
=
{
'Emax'
:
Emax
,
'sampleType'
:
'Krylov'
,
'POD'
:
False
}
paramsSetsPade
=
[
None
]
*
nsets
paramsSetsRB
=
[
None
]
*
nsets
for
i
in
range
(
nsets
):
paramsSetsPade
[
i
]
=
{
'N'
:
stride
*
i
+
shift
+
1
,
'M'
:
stride
*
i
+
shift
,
'E'
:
stride
*
i
+
shift
+
1
}
paramsSetsRB
[
i
]
=
{
'E'
:
stride
*
i
+
shift
,
'R'
:
stride
*
i
+
shift
+
1
}
approxPade
=
TP
(
solver
,
plotter
,
k0
=
kappa
,
approxParameters
=
params
)
approxRB
=
TRB
(
solver
,
plotter
,
k0
=
kappa
,
approxParameters
=
params
)
else
:
kLeft
,
kRight
=
8
,
12
shift
=
3
nsets
=
3
stride
=
3
Smax
=
stride
*
(
nsets
-
1
)
+
shift
+
2
paramsPade
=
{
'S'
:
Smax
,
'polyBasis'
:
'CHEBYSHEV'
,
'POD'
:
True
}
paramsRB
=
{
'S'
:
Smax
,
'nodesType'
:
'CHEBYSHEV'
,
'POD'
:
True
}
paramsSetsPade
=
[
None
]
*
nsets
paramsSetsRB
=
[
None
]
*
nsets
for
i
in
range
(
nsets
):
paramsSetsPade
[
i
]
=
{
'N'
:
stride
*
i
+
shift
+
1
,
'M'
:
stride
*
i
+
shift
,
'S'
:
stride
*
i
+
shift
+
2
}
paramsSetsRB
[
i
]
=
{
'R'
:
stride
*
i
+
shift
+
1
,
'S'
:
stride
*
i
+
shift
+
2
}
approxPade
=
LP
(
solver
,
plotter
,
ks
=
[
kLeft
,
kRight
],
w
=
kappa
,
approxParameters
=
paramsPade
)
approxRB
=
LRB
(
solver
,
plotter
,
ks
=
[
kLeft
,
kRight
],
w
=
kappa
,
approxParameters
=
paramsRB
)
sweeper
=
Sweeper
.
ROMApproximantSweeper
(
ktars
=
ktars
,
mostExpensive
=
'Approx'
)
sweeper
.
ROMEngine
=
approxPade
sweeper
.
params
=
paramsSetsPade
filenamePade
=
sweeper
.
sweep
(
'../Data/HelmholtzAirfoil'
+
test
[:
-
5
]
+
'PadeFF.dat'
)
sweeper
.
ROMEngine
=
approxRB
sweeper
.
params
=
paramsSetsRB
filenameRB
=
sweeper
.
sweep
(
'../Data/HelmholtzAirfoil'
+
test
[:
-
5
]
+
'RBFF.dat'
)
for
i
in
range
(
nsets
):
if
test
==
"TaylorSweep"
:
val
=
stride
*
i
+
shift
PadeOutput
=
sweeper
.
read
(
filenamePade
,
{
'M'
:
val
},
[
'kRe'
,
'HFNorm'
,
'AppNorm'
,
'ErrNorm'
])
RBOutput
=
sweeper
.
read
(
filenameRB
,
{
'E'
:
val
},
[
'kRe'
,
'AppNorm'
,
'ErrNorm'
])
let
=
'E'
else
:
val
=
stride
*
i
+
shift
+
2
PadeOutput
=
sweeper
.
read
(
filenamePade
,
{
'S'
:
val
},
[
'kRe'
,
'HFNorm'
,
'AppNorm'
,
'ErrNorm'
])
RBOutput
=
sweeper
.
read
(
filenameRB
,
{
'S'
:
val
},
[
'kRe'
,
'AppNorm'
,
'ErrNorm'
])
let
=
'S'
ktarsF
=
PadeOutput
[
'kRe'
]
solNormF
=
PadeOutput
[
'HFNorm'
]
PadektarsF
=
PadeOutput
[
'kRe'
]
PadeNormF
=
PadeOutput
[
'AppNorm'
]
PadeErrorF
=
PadeOutput
[
'ErrNorm'
]
RBktarsF
=
RBOutput
[
'kRe'
]
RBNormF
=
RBOutput
[
'AppNorm'
]
RBErrorF
=
RBOutput
[
'ErrNorm'
]
plt
.
figure
(
figsize
=
(
10
,
5
))
plt
.
plot
(
ktarsF
,
solNormF
,
'k-'
,
label
=
'Sol norm'
)
plt
.
plot
(
PadektarsF
,
PadeNormF
,
'b--'
,
label
=
'Pade'' norm, {1} = {0}'
.
format
(
val
,
let
))
plt
.
plot
(
RBktarsF
,
RBNormF
,
'g--'
,
label
=
'RB norm, {1} = {0}'
.
format
(
val
,
let
))
plt
.
legend
()
plt
.
grid
()
p
=
plt
.
legend
(
loc
=
'upper left'
)
plt
.
savefig
(
'./normA'
+
str
(
i
)
+
'.eps'
,
format
=
'eps'
,
dpi
=
1000
)
plt
.
figure
(
figsize
=
(
10
,
5
))
plt
.
semilogy
(
PadektarsF
,
PadeErrorF
,
'b'
,
label
=
'Pade'' error, {1} = {0}'
.
format
(
val
,
let
))
plt
.
semilogy
(
RBktarsF
,
RBErrorF
,
'g'
,
label
=
'RB error, {1} = {0}'
.
format
(
val
,
let
))
plt
.
legend
()
plt
.
grid
()
p
=
plt
.
legend
(
loc
=
'lower right'
)
plt
.
savefig
(
'./errorA'
+
str
(
i
)
+
'.eps'
,
format
=
'eps'
,
dpi
=
1000
)
plt
.
figure
(
figsize
=
(
10
,
5
))
plt
.
semilogy
(
ktarsF
,
np
.
divide
(
PadeErrorF
,
solNormF
),
'b'
,
label
=
'Pade'' relative error, {1} = {0}'
.
format
(
val
,
let
))
plt
.
semilogy
(
RBktarsF
,
np
.
divide
(
RBErrorF
,
solNormF
),
'g'
,
label
=
'RB relative error, {1} = {0}'
.
format
(
val
,
let
))
plt
.
legend
()
plt
.
grid
()
p
=
plt
.
legend
(
loc
=
'lower right'
)
plt
.
savefig
(
'./errorAR'
+
str
(
i
)
+
'.eps'
,
format
=
'eps'
,
dpi
=
1000
)
Event Timeline
Log In to Comment