Page MenuHomec4science
Files F63791793

HelmholtzMembraneTaylor.py
No OneTemporary
Actions

File Metadata

Created
Wed, May 22, 12:29

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{}\nSee 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{}\nErrPade:\t{}\nErrRelPade:\t{}'.format(
solNorm, PadeErr, np.divide(PadeErr, solNorm)))
print('\nPoles 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

c4science · Help