Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/examples/pod/scatteringSquare.py
88e97ff5f796feature_greedy_coarseni…

Last Change
Raw File

scatteringSquare.py
View Options

from copy import copy
import numpy as np
from rrompy.hfengines.scipy import (
HelmholtzCavityScatteringProblemEngine as CSPE)
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as TP
from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as LP
from rrompy.reduction_methods.taylor import ApproximantTaylorRB as TRB
from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as LRB
from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper
from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
from operator import itemgetter
def subdict(d, ks):
return dict(zip(ks, itemgetter(*ks)(d)))
verb = 0
####################
test = "solve"
test = "Taylor"
test = "Lagrange"
test = "TaylorSweep"
#test = "LagrangeSweep"
plotSamples = True
k0 = 10
kLeft, kRight = 9, 11
ktar = 9.5
ktars = np.linspace(8.5, 11.5, 125)
#ktars = np.array([k0])
kappa = 5
n = 50
solver = CSPE(kappa = kappa, n = n, verbosity = verb)
solver.omega = k0
if test == "solve":
uh = solver.solve(k0)
print(solver.norm(uh))
solver.plot(uh, what = ['ABS', 'REAL'])
elif test in ["Taylor", "Lagrange"]:
if test == "Taylor":
params = {'N':8, 'M':7, 'R':8, 'E':8, 'sampleType':'Krylov',
'POD':True}
params = {'N':8, 'M':7, 'R':8, 'E':8, 'sampleType':'Arnoldi',
'POD':True}
parPade = subdict(params, ['N', 'M', 'E', 'sampleType', 'POD'])
parRB = subdict(params, ['R', 'E', 'sampleType', 'POD'])
approxPade = TP(solver, mu0 = k0, approxParameters = parPade,
verbosity = verb)
approxRB = TRB(solver, mu0 = k0, approxParameters = parRB,
verbosity = verb)
else:
params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True, 'basis':"MONOMIAL",
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True, 'basis':"CHEBYSHEV",
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
parPade = subdict(params, ['N', 'M', 'S', 'POD', 'sampler', 'basis'])
parRB = subdict(params, ['R', 'S', 'POD', 'sampler'])
approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parPade,
verbosity = verb)
approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parRB,
verbosity = verb)
approxPade.setupApprox()
approxRB.setupApprox()
if plotSamples:
approxPade.plotSamples()
PadeErr, solNorm = approxPade.normErr(ktar), approxPade.normHF(ktar)
RBErr = approxRB.normErr(ktar)
print(('SolNorm:\t{}\nErrPade:\t{}\nErrRelPade:\t{}\nErrRB:\t\t{}'
'\nErrRelRB:\t{}').format(solNorm, PadeErr,
np.divide(PadeErr, solNorm), RBErr,
np.divide(RBErr, solNorm)))
print('\nPoles Pade'':')
print(approxPade.getPoles())
print('\nPoles RB:')
print(approxRB.getPoles())
approxPade.plotHF(ktar, name = 'u_ex')
approxPade.plotApp(ktar, name = 'u_Pade''')
approxRB.plotApp(ktar, name = 'u_RB')
approxPade.plotErr(ktar, name = 'errPade''')
approxRB.plotErr(ktar, name = 'errRB')
elif test in ["TaylorSweep", "LagrangeSweep"]:
if test == "TaylorSweep":
shift = 5
nsets = 4
stride = 3
Emax = stride * (nsets - 1) + shift + 1
params = {'Emax':Emax, 'sampleType':'Krylov', 'POD':True}
params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True}
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N': stride*i+shift+1,
'M': stride*i+shift,#+1,
'E': stride*i+shift+1}
paramsSetsRB[i] = {'R': stride*i+shift+1,#+1,
'E': stride*i+shift+1}
approxPade = TP(solver, mu0 = k0,approxParameters = params,
verbosity = verb)
approxRB = TRB(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
else:
shift = 7
nsets = 4
stride = 3
Smax = stride * (nsets - 1) + shift + 2
paramsPade = {'S':Smax, 'POD':True, 'basis':"MONOMIAL",
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
paramsPade = {'S':Smax, 'POD':True, 'basis':"CHEBYSHEV",
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
paramsRB = copy(paramsPade)
paramsPoly = copy(paramsPade)
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
paramsSetsPoly = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N': stride*i+shift+1,
'M': stride*i+shift+1,
'S': stride*i+shift+2}
paramsSetsRB[i] = {'R': stride*i+shift+2,
'S': stride*i+shift+2}
paramsSetsPoly[i] = {'N': 0,
'M': stride*i+shift+1,
'S': stride*i+shift+2}
approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = paramsPade,
verbosity = verb)
approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = paramsRB,
verbosity = verb)
approxPoly = LP(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = paramsPoly,
verbosity = verb)
filenamebase = '../data/output/scatSquare' + test[:-5]
sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
sweeper.ROMEngine = approxPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep(filenamebase + 'Pade.dat')
sweeper.ROMEngine = approxRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep(filenamebase + 'RB.dat')
sweeper.ROMEngine = approxPoly
sweeper.params = paramsSetsPoly
filenamePoly = sweeper.sweep(filenamebase + 'Poly.dat')
if test == "TaylorSweep":
constr = ['E']
else:
constr = ['S']
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normHF', 'normApp'], constr, onePlot = True,
save = filenamebase + 'Norm',
saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normResRel'], constr, save = filenamebase + 'Res',
saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normErrRel'], constr, save = filenamebase + 'Err',
saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
c4science · Help