Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F90034395
HelmholtzMembraneTaylor.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
Mon, Oct 28, 16:47
Size
5 KB
Mime Type
text/x-python
Expires
Wed, Oct 30, 16:47 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21997059
Attached To
R6746 RationalROMPy
HelmholtzMembraneTaylor.py
View Options
#!/usr/bin/env python3
import
numpy
as
np
from
context
import
FreeFemHelmholtzEngine
as
HF
from
context
import
FreeFemHSEngine
as
HS
from
context
import
ROMApproximantTaylorPade
as
TP
from
context
import
ROMApproximantSweeper
as
Sweeper
from
matplotlib
import
pyplot
as
plt
####################
test
=
"solve"
#test = "Taylor"
#test = "TaylorSweep"
ttype
=
"simple"
plotSamples
=
'ABS'
#plotSamples = []
z0
=
100
+
1.j
ztar
=
95
ztars
=
np
.
linspace
(
78
,
122
,
5
)
####################
def
boundaryNeumann
(
x
,
on_boundary
):
return
on_boundary
and
x
[
1
]
>
.
25
and
x
[
0
]
>
0.995
and
x
[
0
]
<
1.005
PI
=
np
.
pi
V
=
(
"int n = 10;
\n
"
"real x0 = .5, y0 = .5, A = 1., B = 0., Rr = .1, Ri = .1;
\n
"
"border B1(t = 0, 1){x = 2. * t; y = 0.; label = 1;}
\n
"
"border B2(t = 0, 1){x = 2.; y = t; label = 1;}
\n
"
"border B3(t = 0, 1){x = 2. - .995 * t; y = 1.; label = 1;}
\n
"
"border B4(t = 0, 1){x = 1.005 - .0045 * t; y = 1. - .5 * t; label = 2;}
\n
"
"border B45(t = 0, 1){x = 1.0005 - .001 * t; y = .5; label = 2;}
\n
"
"border B5(t = 0, 1){x = .9995 - .0045 * t; y = .5 + .5 * t; label = 2;}
\n
"
"border B6(t = 0, 1){x = .995 - .995 * t; y = 1.; label = 1;}
\n
"
"border B7(t = 0, 1){x = 0.; y = 1. - t; label = 1;}
\n
"
"mesh Th = buildmesh(B1(4*n) + B2(2*n) + B3(3*n) + B4(3*n) + B45(1) + B5(3*n) + B6(3*n) + B7(2*n));
\n
"
"load
\"
Element_P3
\"
;
\n
"
"fespace V(Th, P3);"
)
f
=
(
"(A * exp(- ((x - x0)^2 + (y - y0)^2) / (2 * Rr^2)) + "
"1.i * B * exp(- ((x - x0)^2 + (y - y0)^2) / (2 * Ri^2)))"
)
solver
=
HF
(
V
,
ztar
**.
5
,
forcingTerm
=
f
,
DirichletBoundary
=
"1"
,
NeumannBoundary
=
"2"
)
plotter
=
HS
(
solver
.
V
)
plotter
.
plotmesh
()
def
crack_disp
(
u
):
import
os
import
subprocess
from
context
import
utilities
filename
=
utilities
.
getNewFilename
(
"temp"
,
"edp"
)
uFFfile
=
FFCT
.
npvec2ffvec
(
u
)
with
open
(
filename
,
'w'
)
as
f
:
f
.
write
(
solver
.
V
+
"
\n
"
)
f
.
write
(
"ifstream fin(
\"
{}
\"
);
\n
"
.
format
(
uFFfile
))
f
.
write
(
"V<complex> u;
\n
"
)
f
.
write
(
"fin >> u[];
\n
"
)
f
.
write
(
"cout.precision(15);"
)
f
.
write
(
"cout.scientific << sqrt(int1d(Th, 2)(abs(u)^2));
\n
"
)
bashCommand
=
"FreeFem++ -v 0 {}"
.
format
(
filename
)
process
=
subprocess
.
Popen
(
bashCommand
.
split
(),
stdout
=
subprocess
.
PIPE
)
output
,
error
=
process
.
communicate
()
try
:
func
=
np
.
float
(
output
)
except
:
raise
Exception
((
"Error was encountered in the execution of "
"FreeFem++ code:
\n
{}
\n
See files {} and {}."
)
\
.
format
(
output
.
decode
(
'utf-8'
),
filename
,
uFFfile
))
os
.
remove
(
filename
)
os
.
remove
(
uFFfile
)
return
func
solver
.
functional
=
crack_disp
if
test
==
"solve"
:
uh
=
solver
.
solve
()
plotter
.
plot
(
uh
,
what
=
'REAL'
)
print
(
solver
.
functional
(
uh
))
elif
test
==
"Taylor"
:
params
=
{
'N'
:
6
,
'M'
:
5
,
'E'
:
6
,
'sampleType'
:
'Arnoldi'
,
'POD'
:
True
}
approxPade
=
TP
(
solver
,
plotter
,
k0
=
z0
,
w
=
np
.
real
(
z0
**.
5
),
plotDer
=
plotSamples
,
approxParameters
=
params
)
PadeErr
,
solNorm
=
approxPade
.
approxError
(
ztar
),
approxPade
.
HFNorm
(
ztar
)
print
(
'SolNorm:
\t
{}
\n
ErrPade:
\t
{}
\n
ErrRelPade:
\t
{}'
.
format
(
solNorm
,
PadeErr
,
np
.
divide
(
PadeErr
,
solNorm
)))
print
(
'
\n
Poles Pade'':'
)
print
(
approxPade
.
getPoles
(
True
))
uHF
=
approxPade
.
getHF
(
ztar
)
plotter
.
plot
(
uHF
,
name
=
'u_ex'
)
approxPade
.
plotApp
(
ztar
,
name
=
'u_Pade'''
)
approxPade
.
plotErr
(
ztar
,
name
=
'errPade'''
)
elif
test
==
"TaylorSweep"
:
shift
=
2
nsets
=
3
stride
=
3
Emax
=
stride
*
(
nsets
-
1
)
+
shift
+
1
params
=
{
'Emax'
:
Emax
,
'sampleType'
:
'Arnoldi'
,
'POD'
:
True
}
paramsSetsPade
=
[
None
]
*
nsets
for
i
in
range
(
nsets
):
paramsSetsPade
[
i
]
=
{
'N'
:
stride
*
i
+
shift
+
1
,
'M'
:
stride
*
i
+
shift
,
'E'
:
stride
*
i
+
shift
+
1
}
approxPade
=
TP
(
solver
,
plotter
,
k0
=
z0
,
approxParameters
=
params
)
sweeper
=
Sweeper
.
ROMApproximantSweeper
(
ktars
=
ztars
,
mostExpensive
=
'Approx'
)
sweeper
.
ROMEngine
=
approxPade
sweeper
.
params
=
paramsSetsPade
filenamePade
=
sweeper
.
sweep
(
'../Data/HelmholtzMembraneTPFE.dat'
,
outputs
=
'ALL'
)
for
i
in
range
(
nsets
):
val
=
stride
*
i
+
shift
+
1
PadeOutput
=
sweeper
.
read
(
filenamePade
,
{
'E'
:
val
},
[
'kRe'
,
'HFNorm'
,
'AppNorm'
,
'ErrNorm'
])
ktarsF
=
PadeOutput
[
'kRe'
]
solNormF
=
PadeOutput
[
'HFNorm'
]
PadektarsF
=
PadeOutput
[
'kRe'
]
PadeNormF
=
PadeOutput
[
'AppNorm'
]
PadeErrorF
=
PadeOutput
[
'ErrNorm'
]
plt
.
figure
(
figsize
=
(
10
,
5
))
plt
.
semilogy
(
ktarsF
,
solNormF
,
'k-'
,
label
=
'Sol norm'
)
plt
.
semilogy
(
PadektarsF
,
PadeNormF
,
'b--'
,
label
=
'Pade'' norm, E = {0}'
.
format
(
val
))
plt
.
legend
()
plt
.
grid
()
p
=
plt
.
legend
(
loc
=
'upper left'
)
#plt.savefig('./Membrane/normA' + str(i) + '.eps', format='eps', dpi=1000)
plt
.
show
()
plt
.
close
()
plt
.
figure
(
figsize
=
(
10
,
5
))
plt
.
semilogy
(
PadektarsF
,
PadeErrorF
,
'b'
,
label
=
'Pade'' error, E = {0}'
.
format
(
val
))
plt
.
legend
()
plt
.
grid
()
p
=
plt
.
legend
(
loc
=
'upper left'
)
#plt.savefig('./Membrane/errorA' + str(i) + '.eps', format='eps', dpi=1000)
plt
.
show
()
plt
.
close
()
plt
.
figure
(
figsize
=
(
10
,
5
))
plt
.
semilogy
(
ktarsF
,
np
.
divide
(
PadeErrorF
,
solNormF
),
'b'
,
label
=
'Pade'' relative error, E = {0}'
.
format
(
val
))
plt
.
legend
()
plt
.
grid
()
p
=
plt
.
legend
(
loc
=
'upper left'
)
#plt.savefig('./Membrane/errorAR' + str(i) + '.eps', format='eps', dpi=1000)
plt
.
show
()
plt
.
close
()
Event Timeline
Log In to Comment