Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/Python_FreeFem++/HelmholtzMembraneTaylor.py
6e6c171cf128master

Last Change
Raw File

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()
c4science · Help