Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/fenics/HelmholtzMembraneTaylor.py
6513a8e1212amaster

Last Change
Raw File

HelmholtzMembraneTaylor.py
View Options

import fenics as fen
import numpy as np
import sympy as sp
from rrompy.hfengines.fenics import HelmholtzProblemEngine as HPE
from rrompy.hsengines.fenics import HSEngine as HS
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as TP
from rrompy.reduction_methods.base import ParameterSweeper as Sweeper
test = "poles"
#test = "error"
#test = "norm"
z0 = 100 + 0.j
ztars = np.array([80, 95, 98])
ztarsNorm = np.linspace(78, 122, 250)
def boundaryNeumann(x, on_boundary):
return on_boundary and x[1] > .25 and x[0] > 0.995 and x[0] < 1.005
x0, y0 = .5, .5
Rr, Ri = .1, .1
x, y = sp.symbols('x[0] x[1]', real=True)
fex = sp.exp(- ((x - x0)**2 + (y - y0)**2) / 2 / Rr**2)
meshname = '../data/mesh/crack_coarse.xml'
#meshname = '../data/mesh/crack_fine.xml'
mesh = fen.Mesh(meshname)
forcingTerm = fen.Expression(sp.printing.ccode(sp.simplify(fex)), degree = 3)
solver = HPE()
solver.omega = np.real(z0 ** .5)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.forcingTerm = forcingTerm
solver.NeumannBoundary = boundaryNeumann
solver.DirichletBoundary = 'rest'
plotter = HS(solver.V)
if test == "poles":
appPoles = {}
Emax = 13
params = {'N':6, 'M':0, 'E':6, 'Emax':Emax, 'sampleType':'Arnoldi',
'POD':True}
approxPade = TP(solver, plotter, mu0 = z0, approxParameters = params)
for E in range(6, Emax + 1):
approxPade.E = E
appPoles[E] = np.sort(approxPade.getPoles())
a = fen.dot(fen.grad(solver.u), fen.grad(solver.v)) * fen.dx
A = fen.assemble(a)
fen.DirichletBC(solver.V, fen.Constant(0.),
solver.DirichletBoundary).apply(A)
AMat = fen.as_backend_type(A).mat()
Ar, Ac, Av = AMat.getValuesCSR()
import scipy.sparse as scsp
A = scsp.csr_matrix((Av, Ac, Ar), shape = AMat.size)
m = fen.dot(solver.u, solver.v) * fen.dx
M = fen.assemble(m)
fen.DirichletBC(solver.V, fen.Constant(0.),
solver.DirichletBoundary).apply(M)
MMat = fen.as_backend_type(M).mat()
Mr, Mc, Mv = MMat.getValuesCSR()
import scipy.sparse as scsp
M = scsp.csr_matrix((Mv, Mc, Mr), shape = MMat.size)
poles = scsp.linalg.eigs(A, k = 7, M = M, sigma = 100.,
return_eigenvectors = False)
II = np.argsort(np.abs(poles - z0))
poles = poles[II]
[print('{},{}'.format(np.real(x), np.imag(x)), end = ',') for x in poles]
print()
for E in range(6, Emax + 1):
print(E, end = ',')
[print('{},{}'.format(np.real(x), np.imag(x)), end = ',')\
for x in np.sort(appPoles[E])]
print()
elif test == "error":
M0 = 5
Emax = 13
params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True}
paramsSetsPade = [None] * (Emax - M0 + 1)
for M in range(M0, Emax + 1):
paramsSetsPade[M - M0] = {'N':M0 + 1, 'M':M, 'E':max(M, M0 + 1)}
approxPade = TP(solver, plotter, mu0 = z0, approxParameters = params)
sweeper = Sweeper(mutars = ztars, mostExpensive = 'Approx')
sweeper.ROMEngine = approxPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep('./error.dat', outputs = 'ALL')
elif test == "norm":
params = [{'N':6, 'M':10, 'E':10, 'sampleType':'Arnoldi', 'POD':True}]
approxPade = TP(solver, plotter, mu0 = z0, approxParameters = params[0])
sweeper = Sweeper(mutars = ztarsNorm, mostExpensive = 'Approx')
sweeper.ROMEngine = approxPade
sweeper.params = params
filenamePade = sweeper.sweep('./norm.dat', outputs = ["HFNorm", "AppNorm",
"ErrNorm"])
c4science · Help