Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/from_papers/pod_membrane_centered.py
b9ced405bab8master

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
from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper
verb = 0
test = "poles"
test = "error"
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'
if test == "poles":
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()
elif test == "error":
M0 = 5
Emax = 8
params = {'E':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, mu0 = k0, approxParameters = params,
verbosity = verb)
sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
sweeper.ROMEngine = approxPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep('membrane_error.dat',
outputs = 'ALL')
sweeper.plot(filenamePade, ['muRe'], ['normHF', 'normApprox'], ['M'],
onePlot = True)
sweeper.plot(filenamePade, ['muRe'], ['normErr'], ['M'])
c4science · Help