Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/examples/from_papers/pod_membrane_centered.py
f972e1cdd3b2feature_greedy_coarseni…

Last Change
Raw File

pod_membrane_centered.py
View Options

import fenics as fen
import numpy as np
from rrompy.hfengines.linear_problem import HelmholtzProblemEngine as HPE
from rrompy.reduction_methods.centered import RationalPade as TP
verb = 0
k0 = 10
ktars = np.linspace(78**.5, 122**.5, 50)
def boundaryNeumann(x, on_boundary):
return on_boundary and x[1] > .25 and x[0] > 0.995 and x[0] < 1.005
meshname = '../data/mesh/crack_coarse.xml'
#meshname = '../data/mesh/crack_fine.xml'
mesh = fen.Mesh(meshname)
x, y = fen.SpatialCoordinate(mesh)[:]
x0, y0 = .5, .5
Rr, Ri = .1, .1
forcingTerm = fen.exp(- ((x - x0)**2 + (y - y0)**2) / 2 / Rr**2)
solver = HPE(verbosity = verb)
solver.omega = np.real(k0)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.forcingTerm = forcingTerm
solver.NeumannBoundary = boundaryNeumann
solver.DirichletBoundary = 'rest'
appPoles = {}
Emax = 13
params = {'N':6, 'M':0, 'E':6, 'sampleType':'Arnoldi',
'POD':True}
approxPade = TP(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
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 - k0))
poles = poles[II]
print('Exact', end = ': ')
[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()
c4science · Help