Homec4science
Diffusion RationalROMPy (feature_greedy_coarsening) Browse

RROMPy/examples/pod/parametricDomain.py
51cd9a0a331afeature_greedy_coarseni…

Last Change
Raw File

parametricDomain.py
View Options

import numpy as np
from rrompy.hfengines.scipy import HelmholtzSquareBubbleDomainProblemEngine as HSBDPE
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade
from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB
from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper
from operator import itemgetter
def subdict(d, ks):
return dict(zip(ks, itemgetter(*ks)(d)))
testNo = 3
verb = 0
if testNo == 1:
mu = 7 ** .5
solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, mu0 = mu,
degree_threshold = 15, verbosity = verb)
uh = solver.solve(mu)
solver.plotmesh()
print(solver.norm(uh))
solver.plot(uh)
############
if testNo == 2:
params = {'N':8, 'M':8, 'E':8, 'sampleType':'Arnoldi', 'POD':True}
# params = {'N':7, 'M':8, 'E':8, 'sampleType':'Krylov', 'POD':True}
mu0 = 7 ** .5
mutar = (7. + .1j) ** .5
solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, mu0 = mu0,
degree_threshold = 15, verbosity = verb)
approxP = Pade(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
paramsRB = subdict(params, ['E', 'sampleType', 'POD'])
approxR = RB(solver, mu0 = mu0, approxParameters = paramsRB,
verbosity = verb)
approxP.setupApprox()
approxR.setupApprox()
# approxP.plotSamples()
approxP.plotHF(mutar, name = 'u_HF')
approxP.plotApp(mutar, name = 'u_Pade''')
approxR.plotApp(mutar, name = 'u_RB')
approxP.plotErr(mutar, name = 'err_Pade''')
approxR.plotErr(mutar, name = 'err_RB')
solNorm = approxP.normHF(mutar)
appPErr = approxP.normErr(mutar)
appRErr = approxR.normErr(mutar)
print(('SolNorm:\t{}\nErrRelP:\t{}\nErrRelR:\t{}').format(solNorm,
appPErr / solNorm, appRErr / solNorm))
print('\nPoles Pade'':')
print(approxP.getPoles())
############
elif testNo == 3:
mu0 = 14 ** .5
mutars = np.linspace(9**.5, 19**.5, 100)
solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20, mu0 = mu0,
degree_threshold = 15, verbosity = verb)
shift = 10
nsets = 1
stride = 2
Emax = stride * (nsets - 1) + shift
params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True}
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
paramsSetsPoly = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N':stride*i+shift, 'M':stride*i+shift,
'E':stride*i+shift}
paramsSetsRB[i] = {'E':stride*i+shift}
paramsSetsPoly[i] = {'N':0, 'M':stride*i+shift,
'E':stride*i+shift}
approxPade = Pade(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approxRB = RB(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approxPoly = Pade(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
filenamebase = '../data/output/domainTaylor'
sweeper = Sweeper(mutars = mutars, 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')
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normHF', 'normApp'], ['E'], onePlot = True,
save = filenamebase + 'Norm',
saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normResRel'], ['E'], save = filenamebase + 'Res',
saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normErrRel'], ['E'], save = filenamebase + 'Err',
saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
c4science · Help