diff --git a/examples/fenics/LagrangeSweep.py b/examples/fenics/LagrangeSweep.py
deleted file mode 100644
index 2b6cf79..0000000
--- a/examples/fenics/LagrangeSweep.py
+++ /dev/null
@@ -1,59 +0,0 @@
-from copy import copy
-import numpy as np
-from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.fenics import HelmholtzBoxScatteringProblemEngine as HBSPE
-from rrompy.hsengines.fenics import HSEngine as HS
-from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade
-from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as RB
-from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper
-from rrompy.sampling import QuadratureSampler as QS
-
-testNo = 2
-npoints = 31
-
-if testNo == 1:
- ks = [4 + .5j, 14 + .5j]
- solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 10)
- solver.omega = np.real(np.mean(ks))
- mutars = np.linspace(0, 21, npoints)
- filenamebase = '../data/output/HelmholtzBubbleLagrange'
-elif testNo == 2:
- ks = [10 + .25j, 14 + .25j]
- solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 10)
- solver.omega = np.real(np.mean(ks))
- mutars = np.linspace(9, 14, npoints)
- filenamebase = '../data/output/HelmholtzBoxLagrange'
-plotter = HS(solver.V)
-
-k0 = np.mean(ks)
-shift = 12
-nsets = 5
-stride = 3
-Smax = stride * (nsets - 1) + shift + 2
-paramsPade = {'S':Smax, 'POD':True, 'sampler':QS(ks, "CHEBYSHEV")}
-paramsRB = copy(paramsPade)
-paramsSetsPade = [None] * nsets
-paramsSetsRB = [None] * nsets
-for i in range(nsets):
- paramsSetsPade[i] = {'N': stride * i + shift + 1, 'M': stride * i + shift,
- 'S': stride * i + shift + 2}
- paramsSetsRB[i] = {'R': stride * i + shift + 1, 'S': stride * i + shift + 2}
-
-appPade = Pade(solver, plotter, mu0 = k0, approxParameters = paramsPade)
-appRB = RB(solver, plotter, mu0 = k0, approxParameters = paramsRB)
-
-sweeper = Sweeper(mutars = mutars, mostExpensive = 'Approx')
-sweeper.ROMEngine = appPade
-sweeper.params = paramsSetsPade
-filenamePade = sweeper.sweep(filenamebase + 'PadeFE.dat', outputs = 'ALL')
-
-sweeper.ROMEngine = appRB
-sweeper.params = paramsSetsRB
-filenameRB = sweeper.sweep(filenamebase + 'RBFE.dat', outputs = 'ALL')
-
-sweeper.plot(filenamePade, ['muRe'], ['HFNorm', 'AppNorm'], ['S'],
- onePlot = True)
-sweeper.plot(filenameRB, ['muRe'], ['HFNorm', 'AppNorm'], ['S'],
- onePlot = True)
-sweeper.plot(filenamePade, ['muRe'], ['ErrNorm'], ['S'])
-sweeper.plot(filenameRB, ['muRe'], ['ErrNorm'], ['S'])
diff --git a/examples/fenics/PadeLagrange.py b/examples/fenics/PadeLagrange.py
deleted file mode 100644
index b9f3e45..0000000
--- a/examples/fenics/PadeLagrange.py
+++ /dev/null
@@ -1,78 +0,0 @@
-import numpy as np
-from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.fenics import HelmholtzSquareTransmissionProblemEngine as HSTPE
-from rrompy.hfengines.fenics import HelmholtzBoxScatteringProblemEngine as HBSPE
-from rrompy.hsengines.fenics import HSEngine as HS
-from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade
-from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
-
-testNo = 1
-
-if testNo == 1:
- z0s = [10 + 0.j, 14 + 0.j]
- z0 = np.mean(z0s)
- ztar = 11 + .5j
- params = {'N':4, 'M':3, 'S':5, 'POD':True, 'sampler':QS(z0s, "CHEBYSHEV")}
-
- solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
- solver.omega = np.real(z0**.5)
- plotter = HS(solver.V)
- approx = Pade(solver, plotter, mu0 = z0, approxParameters = params)#,
-# plotSnap = 'ALL')
-
- approx.plotApp(ztar, name = 'u_Pade''')
- approx.plotHF(ztar, name = 'u_HF')
- approx.plotErr(ztar, name = 'err')
-
- appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles Pade'':')
- print(approx.getPoles())
-
-############
-elif testNo == 2:
- z0s = np.power([3.85 + 0.j, 4.15 + 0.j], 2.)
- z0 = np.mean(z0s)
- ztar = (4 + .15j) ** 2.
- params = {'N':9, 'M':8, 'S':10, 'POD':True, 'sampler':QS(z0s, "CHEBYSHEV")}
-
- solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50)
- solver.omega = np.real(z0**.5)
- plotter = HS(solver.V)
- approx = Pade(solver, plotter, mu0 = z0, approxParameters = params,
- plotSnap = 'ALL')
-
- approx.plotApp(ztar, name = 'u_Pade''')
- approx.plotHF(ztar, name = 'u_HF')
- approx.plotErr(ztar, name = 'err')
-
- appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles Pade'':')
- print(approx.getPoles())
-
-############
-elif testNo == 3:
- k0s = [2, 5]
- k0 = np.mean(k0s)
- ktar = 4.5 - .2j
- params = {'N':10, 'M':9, 'S':15, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV")}
-
- solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40)
- solver.omega = np.real(k0)
- plotter = HS(solver.V)
- approx = Pade(solver, plotter, mu0 = k0, approxParameters = params)#,
-# plotSnap = 'ALL')
-
- approx.plotApp(ktar, name = 'u_Pade''')
- approx.plotHF(ktar, name = 'u_HF')
- approx.plotErr(ktar, name = 'err')
-
- appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles Pade'':')
- print(approx.getPoles())
-
diff --git a/examples/fenics/PadeTaylor.py b/examples/fenics/PadeTaylor.py
deleted file mode 100644
index 9b7c813..0000000
--- a/examples/fenics/PadeTaylor.py
+++ /dev/null
@@ -1,75 +0,0 @@
-import numpy as np
-from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.fenics import HelmholtzSquareTransmissionProblemEngine as HSTPE
-from rrompy.hfengines.fenics import HelmholtzBoxScatteringProblemEngine as HBSPE
-from rrompy.hsengines.fenics import HSEngine as HS
-from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade
-
-testNo = 4
-
-if testNo == 1:
- params = {'N':4, 'M':3, 'E':4, 'sampleType':'Arnoldi', 'POD':True}
- z0 = 12
- ztar = 10.5
-
- solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
- solver.omega = np.real(z0 ** .5)
- plotter = HS(solver.V)
- approx = Pade(solver, plotter, mu0 = z0, approxParameters = params)
-
- approx.plotApp(ztar, name = 'u_Pade''')
- approx.plotHF(ztar, name = 'u_HF')
- approx.plotErr(ztar, name = 'err')
-
- appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles Pade'':')
- print(approx.getPoles())
-
-############
-elif testNo == 2:
- params = {'N':7, 'M':6, 'E':7, 'sampleType':'Arnoldi', 'POD':True}
- z0 = 16
- ztar = 14
-
- solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50)
- solver.omega = np.real(z0 ** .5)
- plotter = HS(solver.V)
- approx = Pade(solver, plotter, mu0 = z0, approxParameters = params,
- plotDer = 'ALL')
-
- approx.plotApp(ztar, name = 'u_Pade''')
- approx.plotHF(ztar, name = 'u_HF')
- approx.plotErr(ztar, name = 'err')
-
- appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles Pade'':')
- print(approx.getPoles())
-
-############
-elif testNo in [3, 4]:
- if testNo == 3:
- params = {'N':8, 'M':7, 'E':8, 'sampleType':'Krylov', 'POD':True}
- else:
- params = {'N':8, 'M':7, 'E':8, 'sampleType':'Arnoldi', 'POD':True}
- k0 = 3
- ktar = 4.25+.5j
-
- solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 30)
- solver.omega = np.real(k0)
- plotter = HS(solver.V)
- approx = Pade(solver, plotter, mu0 = k0, approxParameters = params, plotDer = 'ALL')
-
- approx.plotApp(ktar, name = 'u_Pade''')
- approx.plotHF(ktar, name = 'u_HF')
- approx.plotErr(ktar, name = 'err')
-
- appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles Pade'':')
- print(approx.getPoles())
-
diff --git a/examples/fenics/RBLagrange.py b/examples/fenics/RBLagrange.py
deleted file mode 100644
index 2dc844b..0000000
--- a/examples/fenics/RBLagrange.py
+++ /dev/null
@@ -1,77 +0,0 @@
-import numpy as np
-from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.fenics import HelmholtzSquareTransmissionProblemEngine as HSTPE
-from rrompy.hfengines.fenics import HelmholtzBoxScatteringProblemEngine as HBSPE
-from rrompy.hsengines.fenics import HSEngine as HS
-from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as RB
-from rrompy.sampling import QuadratureSampler as QS
-
-testNo = 3
-
-if testNo == 1:
- z0s = [10 + 0.j, 14 + 0.j]
- z0 = np.mean(z0s)
- ztar = 11 + .5j
- params = {'S':5, 'R':4, 'POD':True, 'sampler':QS(z0s, "CHEBYSHEV")}
-
- solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
- solver.omega = np.real(z0**.5)
- plotter = HS(solver.V)
- approx = RB(solver, plotter, mu0 = z0, approxParameters = params)#,
-# plotSnap = 'ALL')
-
- approx.plotApp(ztar, name = 'u_RB')
- approx.plotHF(ztar, name = 'u_HF')
- approx.plotErr(ztar, name = 'err')
-
- appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles RB:')
- print(approx.getPoles())
-
-############
-elif testNo == 2:
- z0s = np.power([3.85 + 0.j, 4.15 + 0.j], 2.)
- z0 = np.mean(z0s)
- ztar = (4 + .15j) ** 2.
- params = {'S':10, 'R':9, 'POD':True, 'sampler':QS(z0s, "CHEBYSHEV")}
-
- solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50)
- solver.omega = np.real(z0**.5)
- plotter = HS(solver.V)
- approx = RB(solver, plotter, mu0 = z0, approxParameters = params,
- plotSnap = 'ALL')
-
- approx.plotApp(ztar, name = 'u_RB')
- approx.plotHF(ztar, name = 'u_HF')
- approx.plotErr(ztar, name = 'err')
-
- appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles RB:')
- print(approx.getPoles())
-
-############
-elif testNo == 3:
- k0s = [2, 5]
- k0 = np.mean(k0s)
- ktar = 4.5 - .2j
- params = {'S':15, 'R':10, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV")}
-
- solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40)
- solver.omega = np.real(k0)
- plotter = HS(solver.V)
- approx = RB(solver, plotter, mu0 = k0, approxParameters = params)
-
- approx.plotApp(ktar, name = 'u_RB')
- approx.plotHF(ktar, name = 'u_HF')
- approx.plotErr(ktar, name = 'err')
-
- appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles RB:')
- print(approx.getPoles())
-
diff --git a/examples/fenics/RBTaylor.py b/examples/fenics/RBTaylor.py
deleted file mode 100644
index e246c61..0000000
--- a/examples/fenics/RBTaylor.py
+++ /dev/null
@@ -1,72 +0,0 @@
-import numpy as np
-from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.fenics import HelmholtzSquareTransmissionProblemEngine as HSTPE
-from rrompy.hfengines.fenics import HelmholtzBoxScatteringProblemEngine as HBSPE
-from rrompy.hsengines.fenics import HSEngine as HS
-from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB
-
-testNo = 2
-
-if testNo == 1:
- params = {'E':4, 'R':4, 'sampleType':'Arnoldi', 'POD':True}
- z0 = 12+.5j
- ztar = 11
-
- solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
- solver.omega = np.real(z0**.5)
- plotter = HS(solver.V)
- approx = RB(solver, plotter, mu0 = z0, approxParameters = params)
-
- approx.plotApp(ztar, name = 'u_RB')
- approx.plotHF(ztar, name = 'u_HF')
- approx.plotErr(ztar, name = 'err')
-
- appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
- print('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}'.format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles RB:')
- print(approx.getPoles())
-
-############
-elif testNo == 2:
- params = {'E':20, 'R':21, 'sampleType':'Arnoldi', 'POD':True}
- z0 = 16
- ztar = 14
-
- solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50)
- solver.omega = np.real(z0**.5)
- plotter = HS(solver.V)
- approx = RB(solver, plotter, mu0 = z0, approxParameters = params,
- plotDer = 'ALL')
-
- approx.plotApp(ztar, name = 'u_RB')
- approx.plotHF(ztar, name = 'u_HF')
- approx.plotErr(ztar, name = 'err')
-
- appErr, solNorm = approx.approxError(ztar), approx.HFNorm(ztar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles RB:')
- print(approx.getPoles())
-
-############
-elif testNo == 3:
- params = {'E':8, 'R':8, 'sampleType':'Krylov', 'POD':True}
- k0 = 3
- ktar = 4+.5j
-
- solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 30)
- solver.omega = np.real(k0)
- plotter = HS(solver.V)
- approx = RB(solver, plotter, mu0 = k0, approxParameters = params)
-
- approx.plotApp(ktar, name = 'u_RB')
- approx.plotHF(ktar, name = 'u_HF')
- approx.plotErr(ktar, name = 'err')
-
- appErr, solNorm = approx.approxError(ktar), approx.HFNorm(ktar)
- print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
- np.divide(appErr, solNorm)))
- print('\nPoles RB:')
- print(approx.getPoles())
-
diff --git a/examples/fenics/TaylorPoles.py b/examples/fenics/TaylorPoles.py
deleted file mode 100644
index 538900f..0000000
--- a/examples/fenics/TaylorPoles.py
+++ /dev/null
@@ -1,63 +0,0 @@
-import numpy as np
-from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hsengines.fenics import HSEngine as HS
-from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade
-from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB
-from rrompy.utilities import squareResonances
-
-z0 = 12+1.j
-Nmin, Nmax = 2, 10
-Nvals = np.arange(Nmin, Nmax + 1, 2)
-
-params = {'N':Nmin, 'M':0, 'Emax':Nmax, 'POD':True, 'sampleType':'Arnoldi'}
-#, 'robustTol':1e-14}
-
-#boolCon = lambda x : np.abs(np.imag(x)) < 1e-1 * np.abs(np.real(x) - np.real(z0))
-#cleanupParameters = {'boolCondition':boolCon, 'residueCheck':True}
-
-solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 25)
-solver.omega = np.real(z0**.5)
-plotter = HS(solver.V)
-approxP = Pade(solver, plotter, mu0 = z0, approxParameters = params)#,
-# equilibration = True, cleanupParameters = cleanupParameters)
-approxR = RB(solver, plotter, mu0 = z0, approxParameters = params)
-
-rP, rE = [None] * len(Nvals), [None] * len(Nvals)
-
-verbose = 1
-for j, N in enumerate(Nvals):
- if verbose > 0:
- print('N = E = {}'.format(N))
- approxP.approxParameters = {'N':N, 'E':N}
- approxR.approxParameters = {'R':N, 'E':N}
- if verbose > 1:
- print(approxP.approxParameters)
- print(approxR.approxParameters)
-
- rP[j] = approxP.getPoles()
- rE[j] = approxR.getPoles()
- if verbose > 2:
- print(rP)
- print(rE)
-
-from matplotlib import pyplot as plt
-plotRows = int(np.ceil(len(Nvals) / 3))
-fig, axes = plt.subplots(plotRows, 3, figsize = (15, 3.5 * plotRows))
-for j, N in enumerate(Nvals):
- i1, i2 = int(np.floor(j / 3)), j % 3
- axes[i1, i2].set_title('N = E = {}'.format(N))
- axes[i1, i2].plot(np.real(rP[j]), np.imag(rP[j]), 'Xb',
- label="Pade'", markersize = 8)
- axes[i1, i2].plot(np.real(rE[j]), np.imag(rE[j]), '*r',
- label="RB", markersize = 10)
- axes[i1, i2].axhline(linewidth=1, color='k')
- xmin, xmax = axes[i1, i2].get_xlim()
- res = squareResonances(xmin, xmax, False)
- axes[i1, i2].plot(res, np.zeros_like(res), 'ok', markersize = 4)
- axes[i1, i2].grid()
- axes[i1, i2].set_xlim(xmin, xmax)
- axes[i1, i2].axis('equal')
- p = axes[i1, i2].legend()
-plt.tight_layout()
-for j in range((len(Nvals) - 1) % 3 + 1, 3):
- axes[plotRows - 1, j].axis('off')
diff --git a/examples/fenics/TaylorSweep.py b/examples/fenics/TaylorSweep.py
deleted file mode 100644
index 58bc0aa..0000000
--- a/examples/fenics/TaylorSweep.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import numpy as np
-from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.fenics import HelmholtzBoxScatteringProblemEngine as HBSPE
-from rrompy.hsengines.fenics import HSEngine as HS
-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
-
-testNo = 1
-
-z0 = 12 + .25j
-npoints = 31
-
-shift = 5
-nsets = 3
-stride = 2
-Emax = stride * (nsets - 1) + shift + 2
-
-if testNo == 1:
- solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 10)
- solver.omega = np.real(z0**.5)
- params = {'Emax':Emax, 'sampleType':'ARNOLDI', 'POD':True}
- ktars = np.linspace(7, 16, npoints)
- filenamebase = '../data/output/HelmholtzBubbleTaylor'
-elif testNo == 2:
- solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 10)
- solver.omega = np.real(z0**.5)
- params = {'Emax':Emax, 'sampleType':'KRYLOV', 'POD':True}
- ktars = np.linspace(11, 13, npoints)
- filenamebase = '../data/output/HelmholtzBoxTaylor'
-plotter = HS(solver.V)
-
-paramsSetsPade = [None] * nsets
-paramsSetsRB = [None] * nsets
-for i in range(nsets):
- paramsSetsPade[i] = {'N':stride*i+shift+1, 'M':stride*i+shift,
- 'E':stride*i+shift+1}
- paramsSetsRB[i] = {'E':stride*i+shift+1,'R':stride*i+shift+1}
-
-appPade = Pade(solver, plotter, mu0 = z0, approxParameters = params)
-appRB = RB(solver, plotter, mu0 = z0, approxParameters = params)
-
-sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
-sweeper.ROMEngine = appPade
-sweeper.params = paramsSetsPade
-filenamePade = sweeper.sweep(filenamebase + 'PadeFE.dat', outputs = 'ALL')
-
-sweeper.ROMEngine = appRB
-sweeper.params = paramsSetsRB
-filenameRB = sweeper.sweep(filenamebase + 'RBFE.dat', outputs = 'ALL')
-
-sweeper.plot(filenamePade, ['muRe'], ['HFNorm', 'AppNorm'], ['E'],
- onePlot = True)
-sweeper.plot(filenameRB, ['muRe'], ['HFNorm', 'AppNorm'], ['E'],
- onePlot = True)
-sweeper.plot(filenamePade, ['muRe'], ['ErrNorm'], ['E'])
-sweeper.plot(filenameRB, ['muRe'], ['ErrNorm'], ['E'])
-
diff --git a/examples/fenics/airfoil.py b/examples/fenics/airfoil.py
deleted file mode 100644
index 1b78b42..0000000
--- a/examples/fenics/airfoil.py
+++ /dev/null
@@ -1,212 +0,0 @@
-from copy import copy
-import fenics as fen
-import numpy as np
-import sympy as sp
-from rrompy.hfengines.fenics import ScatteringProblemEngine as SPE
-from rrompy.hsengines.fenics import HSEngine as HS
-from rrompy.reduction_methods.taylor import ApproximantTaylorPade as TP
-from rrompy.reduction_methods.taylor import ApproximantTaylorRB as TRB
-from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as LP
-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 matplotlib import pyplot as plt
-from operator import itemgetter
-def subdict(d, ks):
- return dict(zip(ks, itemgetter(*ks)(d)))
-
-####################
-
-test = "solve"
-#test = "Taylor"
-#test = "Lagrange"
-#test = "TaylorSweep"
-#test = "LagrangeSweep"
-
-ttype = "simple"
-###ttype = "augmentedI"
-###ttype = "augmentedM"
-
-plotSamples = 'ALL'
-#plotSamples = []
-
-k0 = 10 + 1.j
-kLeft, kRight = 8 + 1.j, 12 + 1.j
-ktar = 11
-ktars = np.linspace(8, 12, 33) - .5j
-
-PI = np.pi
-R = 2
-def Dboundary(x, on_boundary):
- return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R
-
-kappa = 10
-theta = PI * - 45 / 180.
-mu = 1.1
-epsilon = .1
-
-x, y = sp.symbols('x[0] x[1]', real=True)
-phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
-u0ex = - sp.exp(1.j * phiex)
-
-checkReal = x**2-x+y**2
-rhoroot = ((x**2+y**2)/((x-1)**2+y**2))**.25
-phiroot1 = sp.atan(-y/(x**2-x+y**2)) / 2
-phiroot2 = sp.atan(-y/(x**2-x+y**2)) / 2 - PI * sp.sign(-y/(x**2-x+y**2)) / 2
-kappam1 = (((rhoroot*sp.cos(phiroot1)+.5)**2.+(rhoroot*sp.sin(phiroot1))**2.)/
- ((rhoroot*sp.cos(phiroot1)-1.)**2.+(rhoroot*sp.sin(phiroot1))**2.)
- )**.5 - mu
-kappam2 = (((rhoroot*sp.cos(phiroot2)+.5)**2.+(rhoroot*sp.sin(phiroot2))**2.)/
- ((rhoroot*sp.cos(phiroot2)-1.)**2.+(rhoroot*sp.sin(phiroot2))**2.)
- )**.5 - mu
-Heps1 = .9 * .5 * (1 + kappam1/epsilon + sp.sin(PI*kappam1/epsilon) / PI) + .1
-Heps2 = .9 * .5 * (1 + kappam2/epsilon + sp.sin(PI*kappam2/epsilon) / PI) + .1
-
-mesh = fen.Mesh('../data/mesh/airfoil.xml')
-
-a = fen.Expression(('{0}>=0 ? ({2}>=-{1} ? ({2}<={1} ? {4} : 1) : .1) : '
- '({3}>=-{1} ? ({3}<={1} ? {5} : 1) : .1)')\
- .format(sp.printing.ccode(checkReal), str(epsilon),
- sp.printing.ccode(kappam1),
- sp.printing.ccode(kappam2),
- sp.printing.ccode(Heps1),
- sp.printing.ccode(Heps2)), degree = 3)
-
-DirichletTerm = [sp.printing.ccode(sp.simplify(sp.re(u0ex))),
- sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
-
-if ttype == "simple":
- solver = SPE()
- solver.omega = k0
- solver.V = fen.FunctionSpace(mesh, "P", 3)
- solver.diffusivity = a
- solver.DirichletBoundary = Dboundary
- solver.RobinBoundary = 'rest'
- solver.DirichletDatum = [fen.Expression(x, degree = 3)\
- for x in DirichletTerm]
-
- plotter = HS(solver.V)
-else:
- print('NOPE')
-# if ttype[-1] == "I": constraintType = "IDENTITY"
-# else: constraintType = "MASS"
-# solver = HFSA(mesh = mesh, wavenumber = k0, diffusivity = a,
-# forcingTerm = 0, FEDegree = 3, DirichletBoundary = Dboundary,
-# RobinBoundary = 'rest', DirichletDatum = DirichletTerm,
-# constraintType = constraintType)
-# plotter = HSA(solver.V, d = 2)
-
-baseRe, baseIm = solver.DirichletDatum
-baseRe = fen.project(baseRe, solver.V)
-baseIm = fen.project(baseIm, solver.V)
-uinc = np.array(baseRe.vector()) + 1.j * np.array(baseIm.vector())
-if ttype[: -1] == "augmented":
- print('NOPE')
-# uinc = [uinc, kappa * uinc]
-
-if test == "solve":
- aF = fen.interpolate(a, solver.V)
- av = aF.vector()
- uh = solver.solve(k0)
- print(plotter.norm(uh, kappa))
- if ttype == "simple":
- uhtot = uh - uinc
- else:
- uhtot = [x - y for x, y in zip(uh, uinc)]
- print(plotter.norm(uhtot, kappa))
- plotter.plot(av, what = 'Real', name = 'a')
- plotter.plot(uhtot - uh, what = 'Real', name = 'u_inc')
- plotter.plot(uh, what = 'ABS')
- plotter.plot(uhtot, what = 'ABS', name = 'u + u_inc')
-
-elif test in ["Taylor", "Lagrange"]:
- if test == "Taylor":
- params = {'N':8, 'M':7, 'R':8, 'E':8, 'sampleType':'Krylov', 'POD':False}
- parPade = subdict(params, ['N', 'M', 'E', 'sampleType', 'POD'])
- parRB = subdict(params, ['R', 'E', 'sampleType', 'POD'])
- approxPade = TP(solver, plotter, mu0 = k0, plotDer = plotSamples,
- approxParameters = parPade)
- approxRB = TRB(solver, plotter, mu0 = k0, approxParameters = parRB)
- else:
- params = {'N':8, 'M':7, 'R':9, 'S':9, 'POD':True,
- 'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
- parPade = subdict(params, ['N', 'M', 'S', 'POD', 'sampler'])
- parRB = subdict(params, ['R', 'S', 'POD', 'sampler'])
- approxPade = LP(solver, plotter, mu0 = np.mean([kLeft, kRight]),
- plotSnap = plotSamples, approxParameters = parPade)
- approxRB = LRB(solver, plotter, mu0 = np.mean([kLeft, kRight]),
- approxParameters = parRB)
-
- PadeErr, solNorm = approxPade.approxError(ktar), approxPade.HFNorm(ktar)
- RBErr = approxRB.approxError(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())
-
- uHF = approxPade.getHF(ktar)
- plotter.plot(uHF, 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 = 10
- nsets = 3
- stride = 3
- Emax = stride * (nsets - 1) + shift + 1
- params = {'Emax':Emax, 'sampleType':'Krylov', 'POD':False}
- paramsSetsPade = [None] * nsets
- paramsSetsRB = [None] * nsets
- for i in range(nsets):
- paramsSetsPade[i] = {'N':stride*i+shift + 1, 'M':stride*i+shift,
- 'E':stride*i+shift + 1}
- paramsSetsRB[i] = {'E':stride*i+shift + 1,'R':stride*i+shift + 1}
- approxPade = TP(solver, plotter, mu0 = kappa,approxParameters = params)
- approxRB = TRB(solver, plotter, mu0 = kappa, approxParameters = params)
- else:
- shift = 10
- nsets = 3
- stride = 3
- Smax = stride * (nsets - 1) + shift + 2
- paramsPade = {'S':Smax, 'POD':True,
- 'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
- paramsRB = copy(paramsPade)
- paramsSetsPade = [None] * nsets
- paramsSetsRB = [None] * nsets
- for i in range(nsets):
- paramsSetsPade[i] = {'N': stride*i+shift+1, 'M': stride*i+shift,
- 'S': stride*i+shift+2}
- paramsSetsRB[i] = {'R': stride*i+shift+1, 'S': stride*i+shift+2}
- approxPade = LP(solver, plotter, mu0 = np.mean([kLeft, kRight]),
- approxParameters = paramsPade)
- approxRB = LRB(solver, plotter, mu0 = np.mean([kLeft, kRight]),
- approxParameters = paramsRB)
-
- sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
-
- sweeper.ROMEngine = approxPade
- sweeper.params = paramsSetsPade
- filenamePade = sweeper.sweep('../data/output/airfoil'+test[:-5]+'Pade.dat')
-
- sweeper.ROMEngine = approxRB
- sweeper.params = paramsSetsRB
- filenameRB = sweeper.sweep('../data/output/airfoil'+test[:-5]+'RB.dat')
-
- if test == "TaylorSweep":
- constr = ['E']
- else:
- constr = ['S']
- sweeper.plot(filenamePade, ['muRe'], ['HFNorm', 'AppNorm'], constr,
- onePlot = True)
- sweeper.plot(filenameRB, ['muRe'], ['HFNorm', 'AppNorm'], constr,
- onePlot = True)
- sweeper.plot(filenamePade, ['muRe'], ['ErrNorm'], constr)
- sweeper.plot(filenameRB, ['muRe'], ['ErrNorm'], constr)
diff --git a/examples/fenics/membraneTaylor.py b/examples/fenics/membraneTaylor.py
deleted file mode 100644
index 3c29410..0000000
--- a/examples/fenics/membraneTaylor.py
+++ /dev/null
@@ -1,112 +0,0 @@
-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.utilities.parameter_sweeper 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('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 = 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('../data/output/membrane_error.dat',
- outputs = 'ALL')
- sweeper.plot(filenamePade, ['muRe'], ['HFNorm', 'AppNorm'], ['E'],
- onePlot = True)
- sweeper.plot(filenamePade, ['muRe'], ['ErrNorm'], ['E'])
-
-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('../data/output/membrane_norm.dat',
- outputs = ["HFNorm", "AppNorm", "ErrNorm"])
- sweeper.plot(filenamePade, ['muRe'], ['HFNorm', 'AppNorm'], ['E'],
- onePlot = True)
- sweeper.plot(filenamePade, ['muRe'], ['ErrNorm'], ['E'])
-
diff --git a/examples/fenics/parametricDomain.py b/examples/fenics/parametricDomain.py
deleted file mode 100644
index 09cc680..0000000
--- a/examples/fenics/parametricDomain.py
+++ /dev/null
@@ -1,106 +0,0 @@
-import numpy as np
-from rrompy.hfengines.fenics import HelmholtzSquareBubbleDomainProblemEngine as HSBDPE
-from rrompy.hfengines.fenics import HelmholtzSquareBubbleDomainInverseProblemEngine as HSBDIPE
-from rrompy.hsengines.fenics import HSEngine as HS
-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
-
-# testNo > 0 -> mu = 1/a^2
-# testNo < 0 -> mu = a^2
-
-testNo = -3
-
-if testNo in [1, -1]:
- a2 = 7
- if testNo > 0:
- solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
- mu = 1 / a2
- else:
- solver = HSBDIPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
- mu = a2
-
- uh = solver.solve(mu)
- plotter = HS(solver.V)
- plotter.plotmesh()
- print(plotter.norm(uh))
- plotter.plot(uh)
-
-############
-if testNo in [2, -2]:
- params = {'N':7, 'M':8, 'E':8, 'sampleType':'Arnoldi', 'POD':True}
-# params = {'N':7, 'M':8, 'E':8, 'sampleType':'Krylov', 'POD':True}
- a20 = 7
- a2tar = 8.5
- if testNo > 0:
- solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
- mu0 = 1 / a20
- mutar = 1 / a2tar
- else:
- solver = HSBDIPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
- mu0 = a20
- mutar = a2tar
- plotter = HS(solver.V)
- approxP = Pade(solver, plotter, mu0 = mu0, approxParameters = params,
- plotDer = 'ALL')
- approxR = RB(solver, plotter, mu0 = mu0, approxParameters = params,
- plotDer = 'ALL')
-
- 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.HFNorm(mutar)
- appPErr = approxP.approxError(mutar)
- appRErr = approxR.approxError(mutar)
- print(('SolNorm:\t{}\nErrRelP:\t{}\nErrRelR:\t{}').format(solNorm,
- appPErr / solNorm, appRErr / solNorm))
- print('\nPoles Pade'':')
- print(approxP.getPoles())
- print('\nPoles RB:')
- print(approxR.getPoles())
-
-############
-elif testNo in [3, -3]:
- a20 = 14
- a2tars = np.linspace(9, 19, 250)
- if testNo > 0:
- solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20)
- mu0 = 1 / a20
- mutars = np.power(a2tars, -1.)
- else:
- solver = HSBDIPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20)
- mu0 = a20
- mutars = a2tars
- plotter = HS(solver.V)
- shift = 10
- nsets = 3
- stride = 2
- Emax = stride * (nsets - 1) + shift
- params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True}
- paramsSetsPade = [None] * nsets
- paramsSetsRB = [None] * nsets
- for i in range(nsets):
- paramsSetsPade[i] = {'N':10, 'M':stride*i+shift, 'E':stride*i+shift}
- paramsSetsRB[i] = {'R':11, 'E':stride*i+shift}
- approxPade = Pade(solver, plotter, mu0 = mu0, approxParameters = params)
- approxRB = RB(solver, plotter, mu0 = mu0, approxParameters = params)
-
- sweeper = Sweeper(mutars = mutars, mostExpensive = 'Approx')
- sweeper.ROMEngine = approxPade
- sweeper.params = paramsSetsPade
- filenamePade = sweeper.sweep('../data/output/domain_errorPade.dat',
- outputs = 'ALL')
- sweeper.ROMEngine = approxRB
- sweeper.params = paramsSetsRB
- filenameRB = sweeper.sweep('../data/output/domain_errorRB.dat',
- outputs = 'ALL')
-
- sweeper.plot(filenamePade, ['muRe'], ['HFNorm', 'AppNorm'], ['E'],
- onePlot = True)
- sweeper.plot(filenameRB, ['muRe'], ['HFNorm', 'AppNorm'], ['E'],
- onePlot = True)
- sweeper.plot(filenamePade, ['muRe'], ['ErrNorm'], ['E'])
- sweeper.plot(filenameRB, ['muRe'], ['ErrNorm'], ['E'])
diff --git a/examples/fenics/scatteringSquare.py b/examples/fenics/scatteringSquare.py
deleted file mode 100644
index a1a75cb..0000000
--- a/examples/fenics/scatteringSquare.py
+++ /dev/null
@@ -1,154 +0,0 @@
-from copy import copy
-import numpy as np
-from rrompy.hfengines.fenics import HelmholtzCavityScatteringProblemEngine as CSPE
-from rrompy.hsengines.fenics import HSEngine as HS
-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.sampling import QuadratureSampler as QS
-from operator import itemgetter
-def subdict(d, ks):
- return dict(zip(ks, itemgetter(*ks)(d)))
-
-####################
-
-test = "solve"
-#test = "Taylor"
-#test = "Lagrange"
-test = "TaylorSweep"
-#test = "LagrangeSweep"
-
-ttype = "simple"
-###ttype = "augmentedI"
-###ttype = "augmentedM"
-
-plotSamples = 'ALL'
-#plotSamples = []
-
-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
-
-if ttype == "simple":
- solver = CSPE(kappa = kappa, n = n)
- solver.omega = k0
- plotter = HS(solver.V)
-else:
- print('NOPE')
-# if ttype[-1] == "I": constraintType = "IDENTITY"
-# else: constraintType = "MASS"
-# solver = HFSA(mesh = mesh, wavenumber = k0, diffusivity = a,
-# forcingTerm = 0, FEDegree = 3, DirichletBoundary = Dboundary,
-# RobinBoundary = 'rest', DirichletDatum = DirichletTerm,
-# constraintType = constraintType)
-# plotter = HSA(solver.V, d = 2)
-
-if test == "solve":
- uh = solver.solve(k0)
- print(plotter.norm(uh, kappa))
- plotter.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}
- parPade = subdict(params, ['N', 'M', 'E', 'sampleType', 'POD'])
- parRB = subdict(params, ['R', 'E', 'sampleType', 'POD'])
- approxPade = TP(solver, plotter, mu0 = k0, plotDer = plotSamples,
- approxParameters = parPade)
- approxRB = TRB(solver, plotter, mu0 = k0, approxParameters = parRB)
- else:
- params = {'N':8, 'M':7, 'R':8, 'S':9, 'POD':True,
- 'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
- parPade = subdict(params, ['N', 'M', 'S', 'POD', 'sampler'])
- parRB = subdict(params, ['R', 'S', 'POD', 'sampler'])
- approxPade = LP(solver, plotter, mu0 = np.mean([kLeft, kRight]),
- plotSnap = plotSamples, approxParameters = parPade)
- approxRB = LRB(solver, plotter, mu0 = np.mean([kLeft, kRight]),
- approxParameters = parRB)
-
- PadeErr, solNorm = approxPade.approxError(ktar), approxPade.HFNorm(ktar)
- RBErr = approxRB.approxError(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())
-
- uHF = approxPade.getHF(ktar)
- plotter.plot(uHF, 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}
- 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, plotter, mu0 = k0,approxParameters = params)
- approxRB = TRB(solver, plotter, mu0 = k0, approxParameters = params)
- else:
- shift = 7
- nsets = 4
- stride = 3
- Smax = stride * (nsets - 1) + shift + 2
- paramsPade = {'S':Smax, 'POD':True,
- 'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
- paramsRB = copy(paramsPade)
- paramsSetsPade = [None] * nsets
- paramsSetsRB = [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+1,#+1,
- 'S': stride*i+shift+2}
- approxPade = LP(solver, plotter, mu0 = np.mean([kLeft, kRight]),
- approxParameters = paramsPade)
- approxRB = LRB(solver, plotter, mu0 = np.mean([kLeft, kRight]),
- approxParameters = paramsRB)
-
- sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
-
- sweeper.ROMEngine = approxPade
- sweeper.params = paramsSetsPade
- filenamePade = sweeper.sweep('../data/output/ssquare'+test[:-5]+'Pade.dat',
- outputs = ["Poles"])
-
- sweeper.ROMEngine = approxRB
- sweeper.params = paramsSetsRB
- filenameRB = sweeper.sweep('../data/output/ssquare'+test[:-5]+'RB.dat',
- outputs = ["Poles"])
-
- if test == "TaylorSweep":
- constr = ['E']
- else:
- constr = ['S']
- sweeper.plot(filenamePade, ['muRe'], ['HFNorm', 'AppNorm'], constr,
- onePlot = True)
- sweeper.plot(filenameRB, ['muRe'], ['HFNorm', 'AppNorm'], constr,
- onePlot = True)
- sweeper.plot(filenamePade, ['muRe'], ['ErrNorm'], constr)
- sweeper.plot(filenameRB, ['muRe'], ['ErrNorm'], constr)
diff --git a/examples/fenics/solver.py b/examples/fenics/solver.py
deleted file mode 100644
index 4c71388..0000000
--- a/examples/fenics/solver.py
+++ /dev/null
@@ -1,38 +0,0 @@
-import numpy as np
-from rrompy.hfengines.fenics import HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.fenics import HelmholtzSquareTransmissionProblemEngine as HSTPE
-from rrompy.hfengines.fenics import HelmholtzBoxScatteringProblemEngine as HBSPE
-from rrompy.hsengines.fenics import HSEngine as HS
-
-testNo = 3
-
-if testNo == 1:
- solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40)
-
- uh = solver.solve(12. + 0.j)
- plotter = HS(solver.V)
- plotter.plotmesh(save = True)
- print(plotter.norm(uh))
- plotter.plot(uh, save = True)
-
-###########
-elif testNo == 2:
- solver = HSTPE(nT = 1, nB = 2, theta = np.pi * 20 / 180,
- kappa = 4., n = 50)
-
- uh = solver.solve(16.)
- plotter = HS(solver.V)
- print(plotter.norm(uh))
- plotter.plot(uh)
-
-###########
-elif testNo == 3:
- solver = HBSPE(R = 5, kappa = 12**.5, theta = - np.pi * 60 / 180, n = 30)
- plotter = HS(solver.V)
- uinc = - solver.liftDirichletData()
- uh = solver.solve(12**.5)
- plotter.plotmesh()
- print(plotter.norm(uh, 12**.5))
- print(plotter.norm(uh + uinc, 12**.5))
- plotter.plot(uh)
- plotter.plot(uh + uinc, name = 'u_tot')
diff --git a/rrompy/hfengines/fenics/__init__.py b/rrompy/hfengines/fenics/__init__.py
deleted file mode 100644
index 79c6492..0000000
--- a/rrompy/hfengines/fenics/__init__.py
+++ /dev/null
@@ -1,44 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-from rrompy.hfengines.fenics.generic_problem_engine import GenericProblemEngine
-from rrompy.hfengines.fenics.helmholtz_base_problem_engine import HelmholtzBaseProblemEngine
-from rrompy.hfengines.fenics.helmholtz_box_scattering_problem_engine import HelmholtzBoxScatteringProblemEngine
-from rrompy.hfengines.fenics.helmholtz_cavity_scattering_problem_engine import HelmholtzCavityScatteringProblemEngine
-from rrompy.hfengines.fenics.helmholtz_problem_engine import HelmholtzProblemEngine
-from rrompy.hfengines.fenics.helmholtz_square_bubble_problem_engine import HelmholtzSquareBubbleProblemEngine
-from rrompy.hfengines.fenics.helmholtz_square_bubble_domain_problem_engine import HelmholtzSquareBubbleDomainProblemEngine
-from rrompy.hfengines.fenics.helmholtz_square_bubble_domain_inverse_problem_engine import HelmholtzSquareBubbleDomainInverseProblemEngine
-from rrompy.hfengines.fenics.helmholtz_square_transmission_problem_engine import HelmholtzSquareTransmissionProblemEngine
-from rrompy.hfengines.fenics.scattering_problem_engine import ScatteringProblemEngine
-
-__all__ = [
- 'GenericProblemEngine',
- 'HelmholtzBaseProblemEngine',
- 'HelmholtzBoxScatteringProblemEngine',
- 'HelmholtzCavityScatteringProblemEngine',
- 'HelmholtzProblemEngine',
- 'HelmholtzSquareBubbleProblemEngine',
- 'HelmholtzSquareBubbleDomainProblemEngine',
- 'HelmholtzSquareBubbleDomainInverseProblemEngine',
- 'HelmholtzSquareTransmissionProblemEngine',
- 'ScatteringProblemEngine'
- ]
-
-
-
diff --git a/rrompy/hfengines/fenics/generic_problem_engine.py b/rrompy/hfengines/fenics/generic_problem_engine.py
deleted file mode 100644
index aad494c..0000000
--- a/rrompy/hfengines/fenics/generic_problem_engine.py
+++ /dev/null
@@ -1,87 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-from abc import abstractmethod
-import numpy as np
-import scipy.sparse.linalg as scspla
-from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List
-
-__all__ = ['GenericProblemEngine']
-
-class GenericProblemEngine:
- """
- Generic solver for parametric problems.
-
- Attributes:
- As: Scipy sparse array representation (in CSC format) of As.
- bs: Numpy array representation of bs.
- """
- Aterms, bterms = 1, 1
- mu0Blocks = 0.
-
- functional = lambda self, u: 0.
-
- def name(self) -> str:
- return self.__class__.__name__
-
- def __str__(self) -> str:
- return self.name()
-
- @abstractmethod
- def energyNormMatrix(self) -> Np2D:
- """Sparse matrix (in CSR format) representative of scalar product."""
- pass
-
- @abstractmethod
- def assembleA(self):
- """Assemble matrix blocks of linear system."""
- self.As = [None] * 1
- pass
-
- @abstractmethod
- def assembleb(self):
- """Assemble matrix blocks of linear system."""
- self.bs = [None] * 1
- pass
-
- def A(self, mu:complex) -> Np2D:
- """Assemble matrix of linear system."""
- self.assembleA()
- A = self.As[0]
- for j in range(1, len(self.As)):
- A = A + np.power(mu, j) * self.As[j]
- return A
-
- def b(self, mu:complex) -> Np1D:
- """Assemble RHS of linear system."""
- self.assembleb()
- b = self.bs[0]
- for j in range(1, len(self.bs)):
- b = b + np.power(mu, j) * self.bs[j]
- return b
-
- def solve(self, mu:complex) -> Np1D:
- """Find solution of linear system."""
- return scspla.spsolve(self.A(mu), self.b(mu))
-
- def getLSBlocks(self) -> Tuple[List[Np2D], List[Np1D]]:
- """Get blocks of linear system."""
- self.assembleA()
- self.assembleb()
- return self.As, self.bs
-
diff --git a/rrompy/hfengines/fenics/helmholtz_base_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_base_problem_engine.py
deleted file mode 100644
index b7d0eb9..0000000
--- a/rrompy/hfengines/fenics/helmholtz_base_problem_engine.py
+++ /dev/null
@@ -1,307 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-import scipy.sparse as scsp
-import fenics as fen
-from rrompy.hfengines.fenics.generic_problem_engine import GenericProblemEngine
-from rrompy.utilities.base.types import Np1D, Np2D
-from rrompy.utilities.base.fenics import fenZERO, fenONE, bdrTrue, bdrFalse
-
-__all__ = ['HelmholtzBaseProblemEngine']
-
-class HelmholtzBaseProblemEngine(GenericProblemEngine):
- """
- ABSTRACT
- Solver for generic Helmholtz problems with parametric wavenumber.
- - \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega
- u = u0 on \Gamma_D
- \partial_nu = g1 on \Gamma_N
- \partial_nu + h u = g2 on \Gamma_R
-
- Attributes:
- omega: Value of omega.
- diffusivity: Value of a.
- refractionIndex: Value of n.
- forcingTerm: Value of f.
- DirichletDatum: Value of u0.
- NeumannDatum: Value of g1.
- RobinDatumG: Value of g2.
- DirichletBoundary: Function handle to \Gamma_D.
- NeumannBoundary: Function handle to \Gamma_N.
- RobinBoundary: Function handle to \Gamma_R.
- V: Real FE space.
- u: Generic trial functions for variational form evaluation.
- v: Generic test functions for variational form evaluation.
- ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
- A0: Scipy sparse array representation (in CSC format) of A0.
- A1: Scipy sparse array representation (in CSC format) of A1.
- b0: Numpy array representation of b0.
- As: Scipy sparse array representation (in CSC format) of As.
- bs: Numpy array representation of bs.
- """
- Aterms = 2
- omega = 1.
- _diffusivity = [fenONE, fenZERO]
- _refractionIndex = [fenONE, fenZERO]
- _forcingTerm = [fenZERO, fenZERO]
- _DirichletDatum = [fenZERO, fenZERO]
- _NeumannDatum = [fenZERO, fenZERO]
- _RobinDatumG = [fenZERO, fenZERO]
-
- def __init__(self):
- self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
- self.DirichletBoundary = "ALL"
-
- @property
- def V(self):
- """Value of V."""
- return self._V
- @V.setter
- def V(self, V):
- if hasattr(self, "A0"): del self.A0
- if hasattr(self, "A1"): del self.A1
- if hasattr(self, "b0"): del self.b0
- if not type(V).__name__ == 'FunctionSpace':
- raise Exception("V type not recognized.")
- self._V = V
- self.u = fen.TrialFunction(V)
- self.v = fen.TestFunction(V)
- self.dsToBeSet = True
-
- @property
- def diffusivity(self):
- """Value of a."""
- return self._diffusivity
- @diffusivity.setter
- def diffusivity(self, diffusivity):
- if hasattr(self, "A0"): del self.A0
- if not isinstance(diffusivity, (list, tuple,)):
- diffusivity = [diffusivity, fenZERO]
- self._diffusivity = diffusivity
-
- @property
- def refractionIndex(self):
- """Value of n."""
- return self._refractionIndex
- @refractionIndex.setter
- def refractionIndex(self, refractionIndex):
- if hasattr(self, "A1"): del self.A1
- if not isinstance(refractionIndex, (list, tuple,)):
- refractionIndex = [refractionIndex, fenZERO]
- self._refractionIndex = refractionIndex
-
- @property
- def forcingTerm(self):
- """Value of f."""
- return self._forcingTerm
- @forcingTerm.setter
- def forcingTerm(self, forcingTerm):
- if hasattr(self, "b0"): del self.b0
- if not isinstance(forcingTerm, (list, tuple,)):
- forcingTerm = [forcingTerm, fenZERO]
- self._forcingTerm = forcingTerm
-
- @property
- def DirichletDatum(self):
- """
- Value of u0. Its assignment changes u0Re, u0Im, DirichletBCRe and
- DirichletBCIm.
- """
- return self._DirichletDatum
- @DirichletDatum.setter
- def DirichletDatum(self, DirichletDatum):
- if hasattr(self, "b0"): del self.b0
- if not isinstance(DirichletDatum, (list, tuple,)):
- DirichletDatum = [DirichletDatum, fenZERO]
- self._DirichletDatum = DirichletDatum
-
- @property
- def NeumannDatum(self):
- """Value of g1."""
- return self._NeumannDatum
- @NeumannDatum.setter
- def NeumannDatum(self, NeumannDatum):
- if hasattr(self, "b0"): del self.b0
- if not isinstance(NeumannDatum, (list, tuple,)):
- NeumannDatum = [NeumannDatum, fenZERO]
- self._NeumannDatum = NeumannDatum
-
- @property
- def RobinDatumG(self):
- """Value of g2."""
- return self._RobinDatumG
- @RobinDatumG.setter
- def RobinDatumG(self, RobinDatumG):
- if hasattr(self, "b0"): del self.b0
- if not isinstance(RobinDatumG, (list, tuple,)):
- RobinDatumG = [RobinDatumG, fenZERO]
- self._RobinDatumG = RobinDatumG
-
- @property
- def DirichletBoundary(self):
- """Function handle to DirichletBoundary."""
- return self._DirichletBoundary
- @DirichletBoundary.setter
- def DirichletBoundary(self, DirichletBoundary):
- if hasattr(self, "A0"): del self.A0
- if hasattr(self, "A1"): del self.A1
- if isinstance(DirichletBoundary, (str,)):
- if DirichletBoundary.upper() == "NONE":
- self._DirichletBoundary = bdrFalse
- self.DREST = 0
- elif DirichletBoundary.upper() == "ALL":
- self._DirichletBoundary = bdrTrue
- self.NeumannBoundary = "NONE"
- self.RobinBoundary = "NONE"
- self.DREST = 0
- elif DirichletBoundary.upper() == "REST":
- if self.NREST + self.RREST > 0:
- raise Exception("Only 1 'REST' wildcard can be specified.")
- self._DirichletBoundary = lambda x, on_boundary : (on_boundary
- and not self.NeumannBoundary(x, on_boundary)
- and not self.RobinBoundary(x, on_boundary))
- self.DREST = 1
- else:
- raise Exception("DirichletBoundary label not recognized.")
- elif callable(DirichletBoundary):
- self._DirichletBoundary = DirichletBoundary
- self.DREST = 0
- else:
- raise Exception("DirichletBoundary type not recognized.")
-
- @property
- def NeumannBoundary(self):
- """Function handle to NeumannBoundary."""
- return self._NeumannBoundary
- @NeumannBoundary.setter
- def NeumannBoundary(self, NeumannBoundary):
- if hasattr(self, "A1"): del self.A1
- self.dsToBeSet = True
- if isinstance(NeumannBoundary, (str,)):
- if NeumannBoundary.upper() == "NONE":
- self._NeumannBoundary = bdrFalse
- self.NREST = 0
- elif NeumannBoundary.upper() == "ALL":
- self._NeumannBoundary = bdrTrue
- self.DirichletBoundary = "NONE"
- self.RobinBoundary = "NONE"
- self.NREST = 0
- elif NeumannBoundary.upper() == "REST":
- if self.DREST + self.RREST > 0:
- raise Exception("Only 1 'REST' wildcard can be specified.")
- self._NeumannBoundary = lambda x, on_boundary : (on_boundary
- and not self.DirichletBoundary(x, on_boundary)
- and not self.RobinBoundary(x, on_boundary))
- self.NREST = 1
- else:
- raise Exception("NeumannBoundary label not recognized.")
- elif callable(NeumannBoundary):
- self._NeumannBoundary = NeumannBoundary
- self.NREST = 0
- else:
- raise Exception("DirichletBoundary type not recognized.")
-
- @property
- def RobinBoundary(self):
- """Function handle to RobinBoundary."""
- return self._RobinBoundary
- @RobinBoundary.setter
- def RobinBoundary(self, RobinBoundary):
- if hasattr(self, "A0"): del self.A0
- if hasattr(self, "A1"): del self.A1
- self.dsToBeSet = True
- if isinstance(RobinBoundary, (str,)):
- if RobinBoundary.upper() == "NONE":
- self._RobinBoundary = bdrFalse
- self.RREST = 0
- elif RobinBoundary.upper() == "ALL":
- self._RobinBoundary = bdrTrue
- self.DirichletBoundary = "NONE"
- self.NeumannBoundary = "NONE"
- self.RREST = 0
- elif RobinBoundary.upper() == "REST":
- if self.DREST + self.NREST > 0:
- raise Exception("Only 1 'REST' wildcard can be specified.")
- self._RobinBoundary = lambda x, on_boundary : (on_boundary
- and not self.DirichletBoundary(x, on_boundary)
- and not self.NeumannBoundary(x, on_boundary))
- self.RREST = 1
- else:
- raise Exception("RobinBoundary label not recognized.")
- return
- elif callable(RobinBoundary):
- self._RobinBoundary = RobinBoundary
- self.RREST = 0
- else:
- raise Exception("RobinBoundary type not recognized.")
-
- def autoSetDS(self):
- """Set FEniCS boundary measure based on boundary function handles."""
- if self.dsToBeSet:
- NB = self.NeumannBoundary
- RB = self.RobinBoundary
- class NBoundary(fen.SubDomain):
- def inside(self, x, on_boundary):
- return NB(x, on_boundary)
- class RBoundary(fen.SubDomain):
- def inside(self, x, on_boundary):
- return RB(x, on_boundary)
- boundary_markers = fen.MeshFunction("size_t", self.V.mesh(),
- self.V.mesh().topology().dim() - 1)
- NBoundary().mark(boundary_markers, 0)
- RBoundary().mark(boundary_markers, 1)
- self.ds = fen.Measure("ds", domain = self.V.mesh(),
- subdomain_data = boundary_markers)
- self.dsToBeSet = False
-
- def energyNormMatrix(self) -> Np2D:
- """Sparse matrix (in CSR format) representative of scalar product."""
- normMatFen = fen.assemble((fen.dot(fen.grad(self.u), fen.grad(self.v))
- + np.abs(self.omega)**2 * fen.dot(self.u, self.v)) * fen.dx)
- normMat = fen.as_backend_type(normMatFen).mat()
- return scsp.csr_matrix(normMat.getValuesCSR()[::-1],
- shape = normMat.size)
-
- def liftDirichletData(self) -> Np1D:
- """Lift Dirichlet datum."""
- solLRe = fen.interpolate(self.DirichletDatum[0], self.V)
- solLIm = fen.interpolate(self.DirichletDatum[1], self.V)
- return np.array(solLRe.vector()) + 1.j * np.array(solLIm.vector())
-
- def assembleb(self):
- """Assemble RHS blocks of linear system."""
- if not hasattr(self, "b0"):
- fRe, fIm = self.forcingTerm
- g1Re, g1Im = self.NeumannDatum
- g2Re, g2Im = self.RobinDatumG
- L0Re = (fen.dot(fRe, self.v) * fen.dx
- + fen.dot(g1Re, self.v) * self.ds(0)
- + fen.dot(g2Re, self.v) * self.ds(1))
- L0Im = (fen.dot(fIm, self.v) * fen.dx
- + fen.dot(g1Im, self.v) * self.ds(0)
- + fen.dot(g2Im, self.v) * self.ds(1))
- b0Re = fen.assemble(L0Re)
- b0Im = fen.assemble(L0Im)
- fen.DirichletBC(self.V, self.DirichletDatum[0],
- self.DirichletBoundary).apply(b0Re)
- fen.DirichletBC(self.V, self.DirichletDatum[1],
- self.DirichletBoundary).apply(b0Im)
- self.b0 = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
- self.bs = [self.b0]
-
diff --git a/rrompy/hfengines/fenics/helmholtz_box_scattering_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_box_scattering_problem_engine.py
deleted file mode 100644
index 408ebdc..0000000
--- a/rrompy/hfengines/fenics/helmholtz_box_scattering_problem_engine.py
+++ /dev/null
@@ -1,59 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-import fenics as fen
-from rrompy.hfengines.fenics.scattering_problem_engine import ScatteringProblemEngine
-
-__all__ = ['HelmholtzBoxScatteringProblemEngine']
-
-class HelmholtzBoxScatteringProblemEngine(ScatteringProblemEngine):
- """
- Solver for scattering problem outside a box with parametric wavenumber.
- - \Delta u - omega^2 * n^2 * u = 0 in \Omega
- u = 0 on \Gamma_D
- \partial_nu - i k u = 0 on \Gamma_R
- with exact solution a transmitted plane wave.
- """
- def __init__(self, R:float, kappa:float, theta:float, n:int):
- super().__init__()
-
- self.omega = kappa
-
- import mshr
- scatterer = mshr.Polygon([fen.Point(-1, -.5), fen.Point(1, -.5),
- fen.Point(1, .5), fen.Point(.8, .5),
- fen.Point(.8, -.3), fen.Point(-.8, -.3),
- fen.Point(-.8, .5), fen.Point(-1, .5),])
- mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0, 0), R)-scatterer, n)
- self.V = fen.FunctionSpace(mesh, "P", 3)
-
- self.DirichletBoundary = (lambda x, on_boundary:
- on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R)
- self.RobinBoundary = (lambda x, on_boundary:
- on_boundary and (x[0]**2+x[1]**2)**.5 > .95 * R)
-
- import sympy as sp
- x, y = sp.symbols('x[0] x[1]', real=True)
- phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
- u0ex = - sp.exp(1.j * phiex)
- u0 = [sp.printing.ccode(sp.simplify(sp.re(u0ex))),
- sp.printing.ccode(sp.simplify(sp.im(u0ex)))]
- self.DirichletDatum = [fen.Expression(x, degree = 3) for x in u0]
-
-
diff --git a/rrompy/hfengines/fenics/helmholtz_cavity_scattering_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_cavity_scattering_problem_engine.py
deleted file mode 100644
index 22262ef..0000000
--- a/rrompy/hfengines/fenics/helmholtz_cavity_scattering_problem_engine.py
+++ /dev/null
@@ -1,54 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-import fenics as fen
-from rrompy.hfengines.fenics.scattering_problem_engine import ScatteringProblemEngine
-
-__all__ = ['HelmholtzCavityScatteringProblemEngine']
-
-class HelmholtzCavityScatteringProblemEngine(ScatteringProblemEngine):
- """
- Solver for scattering problem inside a cavity with parametric wavenumber.
- - \Delta u - omega^2 * n^2 * u = 0 in \Omega
- u = 0 on \Gamma_D
- \partial_nu - i k u = 0 on \Gamma_R
- with exact solution a transmitted plane wave.
- """
- def __init__(self, kappa:float, n:int, gamma : float = 0.,
- signR : int = -1):
- super().__init__()
-
- self.omega = kappa
-
- mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(np.pi,np.pi), n, n)
- self.V = fen.FunctionSpace(mesh, "P", 3)
-
- self.RobinBoundary = (lambda x, on_boundary: on_boundary
- and np.isclose(x[1], np.pi))
- self.DirichletBoundary = "REST"
-
- import sympy as sp
- x, y = sp.symbols('x[0] x[1]', real=True)
- wex = 4/np.pi**4 * x * y * (x - np.pi) * (y - 2 * np.pi)
- phiex = signR * kappa * (x * gamma + y)
- uex = wex * sp.exp(-1.j * phiex)
- fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
- forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))),
- sp.printing.ccode(sp.simplify(sp.im(fex)))]
- self.forcingTerm = [fen.Expression(x, degree = 3) for x in forcingTerm]
diff --git a/rrompy/hfengines/fenics/helmholtz_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_problem_engine.py
deleted file mode 100644
index e6aebc5..0000000
--- a/rrompy/hfengines/fenics/helmholtz_problem_engine.py
+++ /dev/null
@@ -1,118 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import scipy.sparse as scsp
-import fenics as fen
-from rrompy.hfengines.fenics.helmholtz_base_problem_engine import HelmholtzBaseProblemEngine, fenZERO
-
-__all__ = ['HelmholtzProblemEngine']
-
-class HelmholtzProblemEngine(HelmholtzBaseProblemEngine):
- """
- Solver for Helmholtz problems with parametric wavenumber.
- - \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega
- u = u0 on \Gamma_D
- \partial_nu = g1 on \Gamma_N
- \partial_nu + h u = g2 on \Gamma_R
-
- Attributes:
- omega: Value of omega.
- diffusivity: Value of a.
- refractionIndex: Value of n.
- forcingTerm: Value of f.
- DirichletDatum: Value of u0.
- NeumannDatum: Value of g1.
- RobinDatumH: Value of h.
- RobinDatumG: Value of g2.
- DirichletBoundary: Function handle to \Gamma_D.
- NeumannBoundary: Function handle to \Gamma_N.
- RobinBoundary: Function handle to \Gamma_R.
- V: Real FE space.
- u: Generic trial functions for variational form evaluation.
- v: Generic test functions for variational form evaluation.
- ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
- A0: Scipy sparse array representation (in CSC format) of A0.
- A1: Scipy sparse array representation (in CSC format) of A1.
- b0: Numpy array representation of b0.
- As: Scipy sparse array representation (in CSC format) of As.
- bs: Numpy array representation of bs.
- """
- _RobinDatumH = [fenZERO, fenZERO]
-
- def __init__(self):
- super().__init__()
- self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
- self.DirichletBoundary = "ALL"
-
- @property
- def RobinDatumH(self):
- """Value of h."""
- return self._RobinDatumH
- @RobinDatumH.setter
- def RobinDatumH(self, RobinDatumH):
- if hasattr(self, "A0"): del self.A0
- if not isinstance(RobinDatumH, (list, tuple,)):
- RobinDatumH = [RobinDatumH, fenZERO]
- self._RobinDatumH = RobinDatumH
-
- def assembleA(self):
- """Assemble matrix blocks of linear system."""
- self.autoSetDS()
- if not hasattr(self, "A0"):
- DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- aRe, aIm = self.diffusivity
- hRe, hIm = self.RobinDatumH
- a0Re = (aRe * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
- + hRe * fen.dot(self.u, self.v) * self.ds(1))
- a0Im = (aIm * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
- + hIm * fen.dot(self.u, self.v) * self.ds(1))
- A0Re = fen.assemble(a0Re)
- A0Im = fen.assemble(a0Im)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
- if not hasattr(self, "A1"):
- DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- nRe, nIm = self.refractionIndex
- n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
- a1Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
- a1Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
- A1Re = fen.assemble(a1Re)
- A1Im = fen.assemble(a1Im)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.A1 = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
- self.As = [self.A0, self.A1]
-
-
diff --git a/rrompy/hfengines/fenics/helmholtz_square_bubble_domain_inverse_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_square_bubble_domain_inverse_problem_engine.py
deleted file mode 100644
index 83d693b..0000000
--- a/rrompy/hfengines/fenics/helmholtz_square_bubble_domain_inverse_problem_engine.py
+++ /dev/null
@@ -1,127 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-import scipy.sparse as scsp
-import fenics as fen
-from rrompy.hfengines.fenics.helmholtz_problem_engine import (
- HelmholtzProblemEngine, fenZERO)
-
-__all__ = ['HelmholtzSquareBubbleDomainInverseProblemEngine']
-
-class HelmholtzSquareBubbleDomainInverseProblemEngine(HelmholtzProblemEngine):
- """
- Solver for square bubble Helmholtz problems with parametric domain heigth.
- - \Delta u - kappa^2 * u = f in \Omega_mu = [0,\pi] x [0,\mu\pi]
- u = 0 on \Gamma_mu = \partial\Omega_mu
- with exact solution square bubble times plane wave.
- """
- bterms = 2
-
- def __init__(self, kappa:float, theta:float, n:int):
- super().__init__()
-
- self.omega = kappa
-
- mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(np.pi,np.pi), n, n)
- self.V = fen.FunctionSpace(mesh, "P", 3)
-
- import sympy as sp
- x, y = sp.symbols('x[0] x[1]', real=True)
- wex = 16/np.pi**4 * x * y * (x - np.pi) * (y - np.pi)
- phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
- uex = wex * sp.exp(-1.j * phiex)
- fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
- forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))),
- sp.printing.ccode(sp.simplify(sp.im(fex)))]
- self.forcingTerm = [fen.Expression(x, degree = 3) for x in forcingTerm]
-
- def assembleA(self):
- """Assemble matrix blocks of linear system."""
- self.autoSetDS()
- if not hasattr(self, "A0"):
- DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- aRe, aIm = self.diffusivity
- a0Re = aRe * fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
- a0Im = aIm * fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
- A0Re = fen.assemble(a0Re)
- A0Im = fen.assemble(a0Im)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
- if not hasattr(self, "A1"):
- DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- aRe, aIm = self.diffusivity
- hRe, hIm = self.RobinDatumH
- nRe, nIm = self.refractionIndex
- n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
- k2Re, k2Im = np.real(self.omega ** 2), np.imag(self.omega ** 2)
- k2n2Re = k2Re * n2Re - k2Im * n2Im
- k2n2Im = k2Re * n2Im + k2Im * n2Re
- a1Re = ((aRe * fen.dot(self.u.dx(0), self.v.dx(0))
- - k2n2Re * fen.dot(self.u, self.v)) * fen.dx
- + hRe * fen.dot(self.u, self.v) * self.ds(1))
- a1Im = ((aIm * fen.dot(self.u.dx(0), self.v.dx(0))
- - k2n2Im * fen.dot(self.u, self.v)) * fen.dx
- + hIm * fen.dot(self.u, self.v) * self.ds(1))
- A1Re = fen.assemble(a1Re)
- A1Im = fen.assemble(a1Im)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.A1 = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
- self.As = [self.A0, self.A1]
-
- def assembleb(self):
- """Assemble RHS blocks of linear system."""
- if not hasattr(self, "b1"):
- fRe, fIm = self.forcingTerm
- g1Re, g1Im = self.NeumannDatum
- g2Re, g2Im = self.RobinDatumG
- L1Re = (fen.dot(fRe, self.v) * fen.dx
- + fen.dot(g1Re, self.v) * self.ds(0)
- + fen.dot(g2Re, self.v) * self.ds(1))
- L1Im = (fen.dot(fIm, self.v) * fen.dx
- + fen.dot(g1Im, self.v) * self.ds(0)
- + fen.dot(g2Im, self.v) * self.ds(1))
- b1Re = fen.assemble(L1Re)
- b1Im = fen.assemble(L1Im)
- fen.DirichletBC(self.V, self.DirichletDatum[0],
- self.DirichletBoundary).apply(b1Re)
- fen.DirichletBC(self.V, self.DirichletDatum[1],
- self.DirichletBoundary).apply(b1Im)
- self.b1 = np.array(b1Re[:] + 1.j * b1Im[:], dtype = np.complex)
- if not hasattr(self, "b0"):
- self.b0 = np.zeros_like(self.b1, dtype = np.complex)
- self.bs = [self.b0, self.b1]
-
diff --git a/rrompy/hfengines/fenics/helmholtz_square_bubble_domain_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_square_bubble_domain_problem_engine.py
deleted file mode 100644
index 5a61d62..0000000
--- a/rrompy/hfengines/fenics/helmholtz_square_bubble_domain_problem_engine.py
+++ /dev/null
@@ -1,102 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-import scipy.sparse as scsp
-import fenics as fen
-from rrompy.hfengines.fenics.helmholtz_problem_engine import (
- HelmholtzProblemEngine, fenZERO)
-
-__all__ = ['HelmholtzSquareBubbleDomainProblemEngine']
-
-class HelmholtzSquareBubbleDomainProblemEngine(HelmholtzProblemEngine):
- """
- Solver for square bubble Helmholtz problems with parametric domain heigth.
- - \Delta u - kappa^2 * u = f in \Omega_mu = [0,\pi] x [0,\mu\pi]
- u = 0 on \Gamma_mu = \partial\Omega_mu
- with exact solution square bubble times plane wave.
- """
- def __init__(self, kappa:float, theta:float, n:int):
- super().__init__()
-
- self.omega = kappa
-
- mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(np.pi,np.pi), n, n)
- self.V = fen.FunctionSpace(mesh, "P", 3)
-
- import sympy as sp
- x, y = sp.symbols('x[0] x[1]', real=True)
- wex = 16/np.pi**4 * x * y * (x - np.pi) * (y - np.pi)
- phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
- uex = wex * sp.exp(-1.j * phiex)
- fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
- forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))),
- sp.printing.ccode(sp.simplify(sp.im(fex)))]
- self.forcingTerm = [fen.Expression(x, degree = 3) for x in forcingTerm]
-
- def assembleA(self):
- """Assemble matrix blocks of linear system."""
- self.autoSetDS()
- if not hasattr(self, "A0"):
- DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- aRe, aIm = self.diffusivity
- hRe, hIm = self.RobinDatumH
- nRe, nIm = self.refractionIndex
- n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
- k2Re, k2Im = np.real(self.omega ** 2), np.imag(self.omega ** 2)
- k2n2Re = k2Re * n2Re - k2Im * n2Im
- k2n2Im = k2Re * n2Im + k2Im * n2Re
- a0Re = ((aRe * fen.dot(self.u.dx(0), self.v.dx(0))
- - k2n2Re * fen.dot(self.u, self.v)) * fen.dx
- + hRe * fen.dot(self.u, self.v) * self.ds(1))
- a0Im = ((aIm * fen.dot(self.u.dx(0), self.v.dx(0))
- - k2n2Im * fen.dot(self.u, self.v)) * fen.dx
- + hIm * fen.dot(self.u, self.v) * self.ds(1))
- A0Re = fen.assemble(a0Re)
- A0Im = fen.assemble(a0Im)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
- if not hasattr(self, "A1"):
- DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- aRe, aIm = self.diffusivity
- a1Re = aRe * fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
- a1Im = aIm * fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
- A1Re = fen.assemble(a1Re)
- A1Im = fen.assemble(a1Im)
- DirichletBC0.zero(A1Re)
- DirichletBC0.zero(A1Im)
- A1ReMat = fen.as_backend_type(A1Re).mat()
- A1ImMat = fen.as_backend_type(A1Im).mat()
- A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
- A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
- self.A1 = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
- shape = A1ReMat.size)
- + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
- shape = A1ImMat.size))
- self.As = [self.A0, self.A1]
-
diff --git a/rrompy/hfengines/fenics/helmholtz_square_bubble_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_square_bubble_problem_engine.py
deleted file mode 100644
index 3626ecd..0000000
--- a/rrompy/hfengines/fenics/helmholtz_square_bubble_problem_engine.py
+++ /dev/null
@@ -1,49 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-import fenics as fen
-from rrompy.hfengines.fenics.helmholtz_problem_engine import HelmholtzProblemEngine
-
-__all__ = ['HelmholtzSquareBubbleProblemEngine']
-
-class HelmholtzSquareBubbleProblemEngine(HelmholtzProblemEngine):
- """
- Solver for square bubble Helmholtz problems with parametric wavenumber.
- - \Delta u - omega^2 * u = f in \Omega
- u = 0 on \Gamma_D
- with exact solution square bubble times plane wave.
- """
- def __init__(self, kappa:float, theta:float, n:int):
- super().__init__()
-
- self.omega = kappa
-
- mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(np.pi,np.pi), n, n)
- self.V = fen.FunctionSpace(mesh, "P", 3)
-
- import sympy as sp
- x, y = sp.symbols('x[0] x[1]', real=True)
- wex = 16/np.pi**4 * x * y * (x - np.pi) * (y - np.pi)
- phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
- uex = wex * sp.exp(-1.j * phiex)
- fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
- forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))),
- sp.printing.ccode(sp.simplify(sp.im(fex)))]
- self.forcingTerm = [fen.Expression(x, degree = 3) for x in forcingTerm]
-
diff --git a/rrompy/hfengines/fenics/helmholtz_square_transmission_problem_engine.py b/rrompy/hfengines/fenics/helmholtz_square_transmission_problem_engine.py
deleted file mode 100644
index 0ec2ada..0000000
--- a/rrompy/hfengines/fenics/helmholtz_square_transmission_problem_engine.py
+++ /dev/null
@@ -1,62 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-import fenics as fen
-from rrompy.hfengines.fenics.helmholtz_problem_engine import HelmholtzProblemEngine
-
-__all__ = ['HelmholtzSquareTransmissionProblemEngine']
-
-class HelmholtzSquareTransmissionProblemEngine(HelmholtzProblemEngine):
- """
- Solver for square transmission Helmholtz problems with parametric
- wavenumber.
- - \Delta u - omega^2 * n^2 * u = 0 in \Omega
- u = 0 on \Gamma_D
- with exact solution a transmitted plane wave.
- """
- def __init__(self, nT:float, nB:float, kappa:float, theta:float, n:int):
- super().__init__()
-
- self.omega = kappa
-
- mesh = fen.RectangleMesh(fen.Point(-np.pi/2, -np.pi/2),
- fen.Point(np.pi/2, np.pi/2), n, n)
- self.V = fen.FunctionSpace(mesh, "P", 3)
-
- import sympy as sp
- dx, dy = np.cos(theta), np.sin(theta)
- Kx = kappa * nB * dx
- Ky = kappa * (nT**2. - (nB * dx)**2. + 0.j)**.5
- T = 2 * kappa * nB * dy / (Ky + kappa * nB * dy)
- x, y = sp.symbols('x[0] x[1]', real=True)
- uT = T * sp.exp(1.j * (Kx*x + Ky*y))
- uB = ( sp.exp(1.j * kappa * nB * (dx*x + dy*y))
- + (T - 1) * sp.exp(1.j * kappa * nB * (dx*x - dy*y)))
- uRe = fen.Expression('x[1]>=0 ? {} : {}'.format(
- sp.printing.ccode(sp.re(uT)),
- sp.printing.ccode(sp.re(uB))),
- degree = 3)
- uIm = fen.Expression('x[1]>=0 ? {} : {}'.format(
- sp.printing.ccode(sp.im(uT)),
- sp.printing.ccode(sp.im(uB))),
- degree = 3)
-
- self.refractionIndex = fen.Expression('x[1] >= 0 ? nT : nB',
- nT = nT, nB = nB, degree = 0)
- self.DirichletDatum = [uRe, uIm]
diff --git a/rrompy/hfengines/fenics/scattering_problem_engine.py b/rrompy/hfengines/fenics/scattering_problem_engine.py
deleted file mode 100644
index a49a573..0000000
--- a/rrompy/hfengines/fenics/scattering_problem_engine.py
+++ /dev/null
@@ -1,113 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import scipy.sparse as scsp
-import fenics as fen
-from rrompy.hfengines.fenics.helmholtz_base_problem_engine import HelmholtzBaseProblemEngine, fenZERO
-
-__all__ = ['ScatteringProblemEngine']
-
-class ScatteringProblemEngine(HelmholtzBaseProblemEngine):
- """
- Solver for scattering problems with parametric wavenumber.
- - \nabla \cdot (a \nabla u) - omega^2 * n**2 * u = f in \Omega
- u = u0 on \Gamma_D
- \partial_nu = g1 on \Gamma_N
- \partial_nu +- i k u = g2 on \Gamma_R
-
- Attributes:
- signR: Sign in ABC.
- omega: Value of omega.
- diffusivity: Value of a.
- refractionIndex: Value of n.
- forcingTerm: Value of f.
- DirichletDatum: Value of u0.
- NeumannDatum: Value of g1.
- RobinDatumG: Value of g2.
- DirichletBoundary: Function handle to \Gamma_D.
- NeumannBoundary: Function handle to \Gamma_N.
- RobinBoundary: Function handle to \Gamma_R.
- V: Real FE space.
- u: Generic trial functions for variational form evaluation.
- v: Generic test functions for variational form evaluation.
- ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
- A0: Scipy sparse array representation (in CSC format) of A0.
- A1: Scipy sparse array representation (in CSC format) of A1.
- A2: Scipy sparse array representation (in CSC format) of A1.
- b0: Numpy array representation of b0.
- As: Scipy sparse array representation (in CSC format) of As.
- bs: Numpy array representation of bs.
- """
- Aterms = 3
- signR = - 1.
-
- def __init__(self):
- self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
- self.DirichletBoundary = "ALL"
-
- def assembleA(self):
- """Assemble matrix blocks of linear system."""
- self.autoSetDS()
- if not hasattr(self, "A0"):
- DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- aRe, aIm = self.diffusivity
- a0Re = aRe * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
- a0Im = aIm * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
- A0Re = fen.assemble(a0Re)
- A0Im = fen.assemble(a0Im)
- DirichletBC0.apply(A0Re)
- DirichletBC0.zero(A0Im)
- A0ReMat = fen.as_backend_type(A0Re).mat()
- A0ImMat = fen.as_backend_type(A0Im).mat()
- A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
- A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
- self.A0 = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
- shape = A0ReMat.size)
- + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
- shape = A0ImMat.size))
- if not hasattr(self, "A1"):
- DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- a1 = fen.dot(self.u, self.v) * self.ds(1)
- A1 = fen.assemble(a1)
- DirichletBC0.zero(A1)
- A1Mat = fen.as_backend_type(A1).mat()
- A1r, A1c, A1v = A1Mat.getValuesCSR()
- self.A1 = self.signR * 1.j * scsp.csr_matrix((A1v, A1c, A1r),
- shape = A1Mat.size)
- if not hasattr(self, "A2"):
- DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
- self.DirichletBoundary)
- nRe, nIm = self.refractionIndex
- n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
- a2Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
- a2Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
- A2Re = fen.assemble(a2Re)
- A2Im = fen.assemble(a2Im)
- DirichletBC0.zero(A2Re)
- DirichletBC0.zero(A2Im)
- A2ReMat = fen.as_backend_type(A2Re).mat()
- A2ImMat = fen.as_backend_type(A2Im).mat()
- A2Rer, A2Rec, A2Rev = A2ReMat.getValuesCSR()
- A2Imr, A2Imc, A2Imv = A2ImMat.getValuesCSR()
- self.A2 = (scsp.csr_matrix((A2Rev, A2Rec, A2Rer),
- shape = A2ReMat.size)
- + 1.j * scsp.csr_matrix((A2Imv, A2Imc, A2Imr),
- shape = A2ImMat.size))
- self.As = [self.A0, self.A1, self.A2]
diff --git a/rrompy/hsengines/__init__.py b/rrompy/hsengines/__init__.py
deleted file mode 100644
index ed60590..0000000
--- a/rrompy/hsengines/__init__.py
+++ /dev/null
@@ -1,18 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
diff --git a/rrompy/hsengines/fenics/__init__.py b/rrompy/hsengines/fenics/__init__.py
deleted file mode 100644
index d43cf54..0000000
--- a/rrompy/hsengines/fenics/__init__.py
+++ /dev/null
@@ -1,25 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-from rrompy.hsengines.fenics.hsengine import HSEngine
-
-__all__ = [
- 'HSEngine'
- ]
-
-
diff --git a/rrompy/hsengines/fenics/hsengine.py b/rrompy/hsengines/fenics/hsengine.py
deleted file mode 100644
index 15fb72e..0000000
--- a/rrompy/hsengines/fenics/hsengine.py
+++ /dev/null
@@ -1,153 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-from matplotlib import pyplot as plt
-import fenics as fen
-from rrompy.utilities.base.types import Np1D, strLst, FenSpace, N2FSExpr
-from rrompy.utilities.base import purgeList, getNewFilename
-
-__all__ = ['HSEngine']
-
-class HSEngine:
- """
- Fenics-based Hilbert space engine.
- """
-
- def __init__(self, V:FenSpace):
- self.V = V
-
- def name(self) -> str:
- return self.__class__.__name__
-
- def __str__(self) -> str:
- return self.name()
-
- def norm(self, u:Np1D, normType : N2FSExpr = "H10") -> float:
- """
- Compute general norm of complex-valued function with given dofs.
-
- Args:
- u: numpy complex array with function dofs.
- normType(optional): Target norm identifier. If matrix, target norm
- is one induced by normType. If number, target norm is weighted
- H^1 norm with given weight. If string, must be recognizable by
- Fenics norm command. Defaults to 'H10'.
-
- Returns:
- Norm of the function (non-negative).
- """
- if type(normType).__name__[-6:] == "matrix":
- return np.abs(u.dot(normType.dot(u).conj())) ** .5
- if isinstance(normType, (int, float)):
- return ( HSEngine.norm(self, u, "H10")**2
- + (normType * HSEngine.norm(self, u, "L2"))**2)**.5
- uAbs = fen.Function(self.V)
- uAbs.vector()[:] = np.array(np.abs(u), dtype = float)
- return fen.norm(uAbs, normType)
-
- def plot(self, u:Np1D, name : str = "u", save : bool = False,
- what : strLst = 'all', **figspecs):
- """
- Do some nice plots of the complex-valued function with given dofs.
-
- Args:
- u: numpy complex array with function dofs.
- name(optional): Name to be shown as title of the plots. Defaults to
- 'u'.
- what(optional): Which plots to do. If list, can contain 'ABS',
- 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
- Defaults to 'ALL'.
- save(optional): Whether to save plot(s). Defaults to False.
- figspecs(optional key args): Optional arguments for matplotlib
- figure creation.
- """
- if isinstance(what, (str,)):
- if what.upper() == 'ALL':
- what = ['ABS', 'PHASE', 'REAL', 'IMAG']
- else:
- what = [what]
- what = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'],
- listname = self.name() + ".what", baselevel = 1)
- if len(what) == 0: return
- if 'figsize' not in figspecs.keys():
- figspecs['figsize'] = (13. * len(what) / 4, 3)
- subplotcode = 100 + len(what) * 10
- if 'REAL' in what and 'IMAG' in what:
- combined_data_ReIm = np.concatenate([np.real(u), np.imag(u)])
- _min = np.amin(combined_data_ReIm)
- _max = np.amax(combined_data_ReIm)
- pars = {'vmin' : _min, 'vmax' : _max}
- else:
- pars = {}
-
- plt.figure(**figspecs)
- plt.jet()
- if 'ABS' in what:
- uAb = fen.Function(self.V)
- uAb.vector()[:] = np.array(np.abs(u), dtype = float)
- subplotcode = subplotcode + 1
- plt.subplot(subplotcode)
- p = fen.plot(uAb, title = "|{0}|".format(name))
- plt.colorbar(p)
- if 'PHASE' in what:
- uPh = fen.Function(self.V)
- uPh.vector()[:] = np.array(np.angle(u), dtype = float)
- subplotcode = subplotcode + 1
- plt.subplot(subplotcode)
- p = fen.plot(uPh, title = "phase({0})".format(name))
- plt.colorbar(p)
- if 'REAL' in what:
- uRe = fen.Function(self.V)
- uRe.vector()[:] = np.array(np.real(u), dtype = float)
- subplotcode = subplotcode + 1
- plt.subplot(subplotcode)
- p = fen.plot(uRe, title = "Re({0})".format(name), **pars)
- plt.colorbar(p)
- if 'IMAG' in what:
- uIm = fen.Function(self.V)
- uIm.vector()[:] = np.array(np.imag(u), dtype = float)
- subplotcode = subplotcode + 1
- plt.subplot(subplotcode)
- p = fen.plot(uIm, title = "Im({0})".format(name), **pars)
- plt.colorbar(p)
- if save:
- plt.savefig(getNewFilename("fig", "eps"), format='eps', dpi=1000)
- plt.show()
- plt.close()
-
- def plotmesh(self, name : str = "Mesh", save : bool = False, **figspecs):
- """
- Do a nice plot of the mesh.
-
- Args:
- u: numpy complex array with function dofs.
- name(optional): Name to be shown as title of the plots. Defaults to
- 'u'.
- save(optional): Whether to save plot(s). Defaults to False.
- figspecs(optional key args): Optional arguments for matplotlib
- figure creation.
- """
- plt.figure(**figspecs)
- fen.plot(self.V.mesh())
- if save:
- plt.savefig(getNewFilename("msh", "eps"), format='eps', dpi=1000)
- plt.show()
- plt.close()
-
-
diff --git a/rrompy/linalg/__init__.py b/rrompy/linalg/__init__.py
deleted file mode 100644
index 941e0f2..0000000
--- a/rrompy/linalg/__init__.py
+++ /dev/null
@@ -1,31 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-from rrompy.linalg.linalg_base_engine import LinAlgBaseEngine
-from rrompy.linalg.linalg_krylov_engine import LinAlgKrylovEngine
-from rrompy.linalg.linalg_arnoldi_engine import LinAlgArnoldiEngine
-from rrompy.linalg.pod_engine import PODEngine
-
-__all__ = [
- 'LinAlgBaseEngine',
- 'LinAlgKrylovEngine',
- 'LinAlgArnoldiEngine',
- 'PODEngine'
- ]
-
-
diff --git a/rrompy/linalg/linalg_arnoldi_engine.py b/rrompy/linalg/linalg_arnoldi_engine.py
deleted file mode 100644
index acda57f..0000000
--- a/rrompy/linalg/linalg_arnoldi_engine.py
+++ /dev/null
@@ -1,124 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-from copy import copy
-import numpy as np
-from rrompy.linalg.pod_engine import PODEngine
-from rrompy.linalg.linalg_krylov_engine import LinAlgKrylovEngine
-from rrompy.utilities.base.types import Np1D, Np2D, List
-
-__all__ = ['LinAlgArnoldiEngine']
-
-class LinAlgArnoldiEngine(LinAlgKrylovEngine):
- """HERE"""
-
- def __init__(self, As:List[Np2D], bs:List[Np1D], energyNormMatrix:Np2D,
- mu0 : complex = 0.):
- super().__init__(As, bs, mu0)
- self.energyNormMatrix = energyNormMatrix
-
- def resetHistory(self):
- super().resetHistory()
- self.HArnoldi = None
- self.RArnoldi = None
- self.RHSs = None
- self.samplesAug = None
-
- @property
- def energyNormMatrix(self):
- """Value of energyNormMatrix. Its assignment resets history."""
- return self._energyNormMatrix
- @energyNormMatrix.setter
- def energyNormMatrix(self, energyNormMatrix):
- self._energyNormMatrix = energyNormMatrix
- self.PODEngine = PODEngine(self.energyNormMatrix)
-
- def preprocesssamples(self):
- ns = self.nsamples
- if ns <= 0 or self.nAs() <= 1: return
- return self.samplesAug[:, ns - 1].reshape((self.nAs() - 1, - 1)).T
-
- def preprocessb(self, bs : Np1D = None, overwrite : bool = False):
- if bs is None: bs = self.bs
- ns = self.nsamples
- if ns < len(bs):
- r = copy(bs[ns])
- else:
- r = np.zeros_like(bs[0])
- if min(len(bs), ns + 1) > 1:
- if ns == 1:
- r = r / self.RArnoldi[0, 0]
- else:
- r = ((r - self.RHSs[:, :ns-1].dot(self.RArnoldi[:ns-1, ns-1]))
- / self.RArnoldi[ns-1, ns-1])
- if overwrite:
- self.RHSs[:, ns - 1] = r
- else:
- if ns == 1:
- self.RHSs = r.reshape((- 1, 1))
- else:
- self.RHSs = np.hstack((self.RHSs, r[:, None]))
- return r
-
- def postprocessu(self, u:Np1D, overwrite : bool = False):
- ns = self.nsamples
- if ns == 0:
- u, h, uAug = self.PODEngine.GS(u, np.empty((0, 0)))
- r = h[0]
- uAug = np.pad(u, ((self.nAs() - 2) * u.size, 0), "constant")
- else:
- uAug = np.concatenate((self.samplesAug[u.size :, ns - 1], u),
- axis = None)
- u, h, uAug = self.PODEngine.GS(u, self.samples[:, : ns], ns, uAug,
- self.samplesAug[:, : ns])
- if overwrite:
- self.HArnoldi[: ns + 1, ns] = h
- if ns > 0:
- r = self.HArnoldi[: ns + 1, 1 : ns + 1].dot(
- self.RArnoldi[: ns, ns - 1])
- self.RArnoldi[: ns + 1, ns] = r
- self.samplesAug[:, ns] = uAug
- else:
- if self.nsamples == 0:
- self.HArnoldi = h.reshape((1, 1))
- self.RArnoldi = r.reshape((1, 1))
- self.samplesAug = uAug.reshape((-1, 1))
- else:
- self.HArnoldi=np.block([[self.HArnoldi, h[:-1, None]],
- [np.zeros((1, self.nsamples)), h[-1]]])
- if ns > 0:
- r = self.HArnoldi[: ns + 1, 1 : ns + 1].dot(
- self.RArnoldi[: ns, ns - 1])
- self.RArnoldi=np.block([[self.RArnoldi, r[:-1, None]],
- [np.zeros((1, self.nsamples)), r[-1]]])
- self.samplesAug = np.hstack((self.samplesAug, uAug[:, None]))
- return u
-
- def preallocateSamples(self, u:Np1D, n:int):
- super().preallocateSamples(u, n)
- h = self.HArnoldi
- r = self.RArnoldi
- saug = self.samplesAug
- self.HArnoldi = np.zeros((n, n), dtype = u.dtype)
- self.HArnoldi[0, 0] = h[0, 0]
- self.RArnoldi = np.zeros((n, n), dtype = u.dtype)
- self.RArnoldi[0, 0] = r[0, 0]
- self.RHSs = np.empty((u.size, n - 1), dtype = u.dtype)
- self.samplesAug = np.empty((saug.shape[0], n), dtype = u.dtype)
- self.samplesAug[:, 0] = saug[:, 0]
-
\ No newline at end of file
diff --git a/rrompy/linalg/linalg_base_engine.py b/rrompy/linalg/linalg_base_engine.py
deleted file mode 100644
index 37b4775..0000000
--- a/rrompy/linalg/linalg_base_engine.py
+++ /dev/null
@@ -1,147 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-from copy import copy
-import scipy.sparse.linalg as spla
-from rrompy.utilities.base.types import Np1D, Np2D, List, Tuple
-
-__all__ = ['LinAlgBaseEngine']
-
-class LinAlgBaseEngine:
- """HERE"""
-
- def __init__(self, As:List[Np2D], bs:List[Np1D], mu0 : complex = 0.):
- self.As = As
- self.bs = bs
- self.mu0 = mu0
-
- def name(self) -> str:
- return self.__class__.__name__
-
- def __str__(self) -> str:
- return self.name()
-
- def resetHistory(self):
- pass
-
- @property
- def As(self):
- """Value of As. Its assignment resets history."""
- return self._As
- @As.setter
- def As(self, As):
- self._As = As
- self.resetHistory()
-
- @property
- def bs(self):
- """Value of bs. Its assignment resets history."""
- return self._bs
- @bs.setter
- def bs(self, bs):
- self._bs = bs
- self.resetHistory()
-
- def nAs(self):
- """Number of terms in As."""
- return len(self._As)
-
- def nbs(self):
- """Number of terms in bs."""
- return len(self._bs)
-
- def shiftCenter(self, mu:complex, As : List[Np2D] = None,
- bs : List[Np1D] = None, update : bool = False,
- which : Np1D = None)-> Tuple[List[Np2D], List[Np1D]]:
- """
- Build linear system from polynomial blocks.
-
- Args:
- mu: Parameter value.
- As: Matrix blocks. If None, set to self.As. Defaults to None.
- bs: RHS blocks. If None, set to self.bs. Defaults to None.
- which: What terms to compute. If None, all are computed. Defaults
- to None.
-
- Returns:
- Shifted matrix and RHS coefficients of system.
- """
- if update and not (As is None and bs is None and which is None):
- raise Exception(("Default values for As, bs, and which are "
- "required for update."))
- if As is None: As = self.As
- if bs is None: bs = self.bs
- lAs, lbs = len(As), len(bs)
- maxL = max(lAs, lbs)
- if which is None: which = np.arange(maxL)
- mu0 = self.mu0
- if np.isclose(mu, mu0):
- return ([copy(As[term]) for term in which],
- [copy(bs[term]) for term in which])
- AsNew = []
- bsNew = []
- for term in which:
- rho = 1.
- if term < lAs:
- AsNew += [copy(As[term])]
- if term < lbs:
- bsNew += [copy(bs[term])]
- for j in range(1, maxL - term):
- rho *= (term + j) / j * (mu - mu0)
- if j < lAs:
- AsNew[-1] += rho * As[term + j]
- if j < lbs:
- bsNew[-1] += rho * bs[term + j]
- if update:
- self.As = AsNew
- self.bs = bsNew
- self.mu0 = mu
- return AsNew, bsNew
-
- def buildLS(self, mu:complex, As : List[Np2D] = None,
- bs : List[Np1D] = None) -> Tuple[Np2D, Np1D]:
- """
- Build linear system from polynomial blocks.
-
- Args:
- mu: Parameter value.
- As: Matrix blocks. If None, set to self.As. Defaults to None.
- bs: RHS blocks. If None, set to self.bs. Defaults to None.
-
- Returns:
- Matrix and RHS of system.
- """
- As, bs = self.shiftCenter(mu, As, bs, which = [0])
- return As[0], bs[0]
-
- def solveLS(self, mu:complex, As : List[Np2D] = None,
- bs : List[Np1D] = None) -> Np1D:
- """
- Solve linear system.
-
- Args:
- mu: Parameter value.
-
- Returns:
- Solution of system.
- """
- A, b = self.buildLS(mu, As, bs)
- u = spla.spsolve(A, b)
- return u
-
diff --git a/rrompy/linalg/linalg_krylov_engine.py b/rrompy/linalg/linalg_krylov_engine.py
deleted file mode 100644
index f2628fb..0000000
--- a/rrompy/linalg/linalg_krylov_engine.py
+++ /dev/null
@@ -1,81 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-from copy import copy
-import numpy as np
-from rrompy.linalg.linalg_base_engine import LinAlgBaseEngine
-from rrompy.utilities.base.types import Np1D, Np2D
-
-__all__ = ['LinAlgKrylovEngine']
-
-class LinAlgKrylovEngine(LinAlgBaseEngine):
- """HERE"""
-
- def resetHistory(self):
- self.samples = None
- self.nsamples = 0
-
- def preprocesssamples(self):
- if self.samples is None or self.nAs() <= 1: return
- if self.nAs() > self.nsamples:
- return self.samples[:, : self.nsamples]
- return self.samples[:, self.nsamples + 1 - self.nAs() : self.nsamples]
-
- def preprocessb(self, bs : Np1D = None, overwrite : bool = False):
- if bs is None: bs = self.bs
- if self.nsamples < len(bs):
- return copy(bs[self.nsamples])
- return np.zeros_like(bs[0])
-
- def postprocessu(self, u:Np1D, overwrite : bool = False):
- return u
-
- def preallocateSamples(self, u:Np1D, n:int):
- self.samples = np.empty((u.size, n), dtype = u.dtype)
- self.samples[:, 0] = u
-
- def nextSample(self, mu:complex, overwrite : bool = False) -> Np1D:
- if not np.isclose(mu, self.mu0):
- self.shiftCenter(mu, update = True)
- samplesOld = self.preprocesssamples()
- RHS = self.preprocessb(overwrite = overwrite)
- for i in range(1, min(self.nAs(), self.nsamples + 1)):
- RHS -= self.As[i].dot(samplesOld[:, - i])
- u = self.postprocessu(self.solveLS(mu, bs = [RHS]),
- overwrite = overwrite)
- if overwrite:
- self.samples[:, self.nsamples] = u
- else:
- if self.nsamples == 0:
- self.samples = u[:, None]
- else:
- self.samples = np.hstack((self.samples, u[:, None]))
- self.nsamples += 1
- return u
-
- def iterSample(self, mu:complex, n:int) -> Np2D:
- if n <= 0:
- raise Exception(("Number of Krylov iterations must be positive."))
- self.resetHistory()
- u = self.nextSample(mu)
- if n > 1:
- self.preallocateSamples(u, n)
- for _ in range(1, n):
- self.nextSample(mu, overwrite = True)
- return self.samples
-
\ No newline at end of file
diff --git a/rrompy/linalg/pod_engine.py b/rrompy/linalg/pod_engine.py
deleted file mode 100644
index 0b014b1..0000000
--- a/rrompy/linalg/pod_engine.py
+++ /dev/null
@@ -1,149 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-import numpy as np
-from copy import copy
-from rrompy.utilities.base.types import Np1D, Np2D, Tuple
-
-__all__ = ['PODEngine']
-
-class PODEngine:
- """
- POD engine for general matrix orthogonalization.
-
- Args/Attributes:
- M: Square matrix associated to scalar product.
- """
-
- def __init__(self, M:Np2D):
- self.M = M
-
- def name(self) -> str:
- return self.__class__.__name__
-
- def __str__(self) -> str:
- return self.name()
-
- def norm(self, a:Np1D) -> float:
- """Compute norm of a vector."""
- return np.sqrt(np.abs(a.conj().T.dot(self.M.dot(a))))
-
- def GS(self, a:Np1D, Q:Np2D, n : int = None,
- aA:Np1D = None, QA:Np2D = None) -> Tuple[Np1D, Np1D, Np1D]:
- """
- Compute 1 Gram-Schmidt step with given projector.
-
- Args:
- a: vector to be projected;
- Q: orthogonal projection matrix;
- n: number of columns of Q to be considered;
- aA: augmented components of vector to be projected;
- QA: augmented components of projection matrix.
-
- Returns:
- Resulting normalized vector, coefficients of a wrt the updated
- basis.
- """
- if n is None:
- n = Q.shape[1]
- if aA is None != QA is None:
- raise Exception(("Either both or none of augmented components "
- "must be provided."))
- r = np.zeros((n + 1,), dtype = a.dtype)
- if n > 0:
- Q = Q[:, : n]
- for j in range(2): # twice is enough!
- nu = Q.conj().T.dot(self.M.dot(a))
- a = a - Q.dot(nu)
- if aA is not None:
- aA = aA - QA.dot(nu)
- r[:-1] = r[:-1] + nu
- r[-1] = self.norm(a)
- if np.isclose(np.abs(r[-1]), 0.):
- r[-1] = 1.
- a = a / r[-1]
- if aA is not None:
- aA = aA / r[-1]
- return a, r, aA
-
- def QRGramSchmidt(self, A:Np2D,
- only_R : bool = False) -> Tuple[Np1D, Np1D]:
- """
- Compute QR decomposition of a matrix through Gram-Schmidt method.
-
- Args:
- A: matrix to be decomposed;
- only_R(optional): whether to skip reconstruction of Q; defaults to
- False.
-
- Returns:
- Resulting orthogonal and upper-triangular factors.
- """
- N = A.shape[1]
- Q = np.zeros_like(A, dtype = A.dtype)
- R = np.zeros((N, N), dtype = A.dtype)
- for k in range(N):
- Q[:, k], R[: k + 1, k], _ = self.GS(A[:, k], Q, k)
- if only_R:
- return R
- return Q, R
-
- def QRHouseholder(self, A:Np2D, Q0 : Np2D = None,
- only_R : bool = False) -> Tuple[Np1D, Np1D]:
- """
- Compute QR decomposition of a matrix through Householder method.
-
- Args:
- A: matrix to be decomposed;
- Q0(optional): initial orthogonal guess for Q; defaults to random;
- only_R(optional): whether to skip reconstruction of Q; defaults to
- False.
-
- Returns:
- Resulting (orthogonal and )upper-triangular factor(s).
- """
- B = copy(A)
- N = B.shape[1]
- V = np.zeros_like(B, dtype = B.dtype)
- R = np.zeros((N, N), dtype = B.dtype)
- if Q0 is None:
- Q = np.zeros_like(B, dtype = B.dtype) + np.random.randn(*(B.shape))
- else:
- Q = copy(Q0)
- for k in range(N):
- if Q0 is None:
- Q[:, k], _, _ = self.GS(Q[:, k], Q, k)
- a = B[:, k]
- R[k, k] = self.norm(a)
- alpha = Q[:, k].conj().T.dot(self.M.dot(a))
- if np.isclose(np.abs(alpha), 0.): s = 1.
- else: s = - alpha / np.abs(alpha)
- Q[:, k] = s * Q[:, k]
- V[:, k], _, _ = self.GS(R[k, k] * Q[:, k] - a, Q, k)
- J = np.arange(k + 1, N)
- vtB = V[:, k].conj().dot(self.M.dot(B[:, J]))
- B[:, J] = B[:, J] - 2 * np.outer(V[:, k], vtB)
- R[k, J] = Q[:, k].conj().T.dot(self.M.dot(B[:, J]))
- B[:, J] = B[:, J] - np.outer(Q[:, k], R[k, J])
- if only_R:
- return R
- for k in range(N - 1, -1, -1):
- J = np.arange(k, N)
- vtQ = V[:, k].conj().T.dot(self.M.dot(Q[:, J]))
- Q[:, J] = Q[:, J] - 2 * np.outer(V[:, k], vtQ)
- return Q, R
diff --git a/rrompy/reduction_methods/base/generic_approximant.py b/rrompy/reduction_methods/base/generic_approximant.py
index 6949dd0..12f8d26 100644
--- a/rrompy/reduction_methods/base/generic_approximant.py
+++ b/rrompy/reduction_methods/base/generic_approximant.py
@@ -1,338 +1,290 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
import numpy as np
from copy import copy
-from rrompy.linalg.linalg_base_engine import LinAlgBaseEngine
-from rrompy.utilities.base.types import (Np1D, Np2D, DictAny, N2FSExpr, HFEng,
- HSEng, Tuple, List, LAEng)
+from rrompy.sampling.base.sampling_engine_base import SamplingEngineBase
+from rrompy.utilities.base.types import Np1D, HS1D, DictAny, HFEng
from rrompy.utilities.base import purgeDict
__all__ = ['GenericApproximant']
class GenericApproximant:
"""
ABSTRACT
ROM approximant computation for parametric problems.
Args:
- HFEngine: HF problem solver. Should have members:
- - energyNormMatrix: assemble sparse matrix (in CSC format)
- associated to weighted H10 norm;
- - problemData: list of HF problem data (excluding k);
- - setProblemData: set HF problem data and k to prescribed values;
- - getLSBlocks: get algebraic system blocks.
- HSEngine: Hilbert space general purpose engine. Should have members:
- - norm: compute Hilbert norm of expression represented as complex
- numpy vector;
- - plot: plot Hilbert expression represented as complex numpy vector.
+ HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
- HSEngine: Hilbert space general purpose engine.
mu0: Default parameter.
w: Weight for norm computation.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots.
POD: Whether to compute POD of snapshots.
- energyNormMatrix: Sparse matrix (in CSC format) assoociated to
- w-weighted H10 norm.
- LAEngine: Linear algebra engine.
+ samplingEngine: Sampling engine.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
- As: List of sparse matrices (in CSC format) representing coefficients
- of linear system matrix wrt mu.
- bs: List of numpy vectors representing coefficients of linear system
- RHS wrt mu.
"""
- def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0,
+ def __init__(self, HFEngine:HFEng, mu0 : complex = 0,
approxParameters : DictAny = {}):
self.HFEngine = HFEngine
- self.HSEngine = HSEngine
self._HFEngine0 = copy(HFEngine)
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["POD"]
self.solLifting = self.HFEngine.liftDirichletData()
self.mu0 = mu0
self.approxParameters = approxParameters
- self.energyNormMatrix = self.HFEngine.energyNormMatrix()
- self.As, self.bs = self.HFEngine.getLSBlocks()
- self.setupLAEngine()
+ self.setupSampling()
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
- def setupLAEngine(self) -> LAEng:
- """Setup linear algebra engine."""
- self.LAEngine = LinAlgBaseEngine(self.As, self.bs,
- mu0 = self.HFEngine.mu0Blocks)
+ def setupSampling(self):
+ """Setup sampling engine."""
+ self.samplingEngine = SamplingEngineBase(self.HFEngine)
@property
def approxParameters(self):
"""Value of approximant parameters."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
if not hasattr(self, "approxParameters"):
self._approxParameters = {}
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
keyList = list(approxParameters.keys())
if "POD" in keyList:
self.POD = approxParameters["POD"]
elif hasattr(self, "POD"):
self.POD = self.POD
else:
self.POD = True
@property
def POD(self):
"""Value of POD."""
return self._POD
@POD.setter
def POD(self, POD):
if hasattr(self, "POD"): PODold = self.POD
else: PODold = -1
self._POD = POD
self._approxParameters["POD"] = self.POD
if PODold != self.POD:
self.resetSamples()
def solveHF(self, mu : complex = None):
"""
Find high fidelity solution with original parameters and arbitrary
parameter.
Args:
mu: Target parameter.
"""
if mu is None: mu = self.mu0
if (not hasattr(self, "lastSolvedHF")
or not np.isclose(self.lastSolvedHF, mu)):
- self.uHF = self.LAEngine.solveLS(mu)
+ self.uHF = self.samplingEngine.solveLS(mu)
self.lastSolvedHF = mu
- def buildLS(self, mu:complex, As : List[Np2D] = None,
- bs : List[Np1D] = None) -> Tuple[Np2D, Np1D]:
- """
- Build linear system from polynomial blocks.
-
- Args:
- mu: Parameter value.
- As: Matrix blocks. If None, set to self.As. Defaults to None.
- bs: RHS blocks. If None, set to self.bs. Defaults to None.
-
- Returns:
- Matrix and RHS of system.
- """
- return self.LAEngine.buildLS(mu, As, bs)
-
def resetSamples(self):
"""Reset samples."""
- if hasattr(self, "LAEngine"):
- self.LAEngine.resetHistory()
+ if hasattr(self, "samplingEngine"):
+ self.samplingEngine.resetHistory()
@abstractmethod
def setupApprox(self):
"""
Setup approximant. (ABSTRACT)
Any specialization should include something like
self.computeDerivatives()
if not self.checkComputedApprox():
...
self.lastApproxParameters = copy(self.approxParameters)
"""
pass
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return (hasattr(self, "lastApproxParameters")
and self.approxParameters == self.lastApproxParameters)
@abstractmethod
def evalApprox(self, mu:complex):
"""
Evaluate approximant at arbitrary parameter. (ABSTRACT)
Any specialization should include something like
self.setupApprox()
self.uApp = ...
Args:
mu: Target parameter.
"""
pass
- def getHF(self, mu:complex) -> Np1D:
+ def getHF(self, mu:complex) -> HS1D:
"""
Get HF solution at arbitrary wavenumber.
Args:
mu: Target parameter.
Returns:
HFsolution as numpy complex vector.
"""
self.solveHF(mu)
return self.uHF
- def getApp(self, mu:complex) -> Np1D:
+ def getApp(self, mu:complex) -> HS1D:
"""
Get approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
Returns:
Approximant as numpy complex vector.
"""
self.evalApprox(mu)
return self.uApp
@abstractmethod
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
pass
- def HFNorm(self, mu:complex, normType : N2FSExpr = None) -> float:
+ def HFNorm(self, mu:complex) -> float:
"""
Compute norm of HF solution at arbitrary wavenumber.
Args:
mu: Target parameter.
- normType(optional): Target norm identifier. Must be recognizable
- by HSEngine norm command. If None, set to
- self.energyNormMatrix. Defaults to None.
Returns:
Target norm of HFsolution.
"""
- if normType is None: normType = self.energyNormMatrix
uHF = self.getHF(mu)
- return self.HSEngine.norm(uHF, normType)
+ return self.HFEngine.norm(uHF)
- def approxNorm(self, mu:complex, normType : N2FSExpr = None) -> float:
+ def approxNorm(self, mu:complex) -> float:
"""
Compute norm of (translated) approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
- normType(optional): Target norm identifier. Must be recognizable
- by HSEngine norm command. If None, set to
- self.energyNormMatrix. Defaults to None.
Returns:
Target norm of approximant.
"""
- if normType is None: normType = self.energyNormMatrix
uApp = self.getApp(mu)
- return self.HSEngine.norm(uApp, normType)
+ return self.HFEngine.norm(uApp)
- def approxError(self, mu:complex, normType : N2FSExpr = None) -> float:
+ def approxError(self, mu:complex) -> float:
"""
Compute norm of approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
- normType(optional): Target norm identifier. Must be recognizable
- by HSEngine norm command. If None, set to
- self.energyNormMatrix. Defaults to None.
Returns:
Target norm of (approximant - HFsolution).
"""
- if normType is None: normType = self.energyNormMatrix
uApp = self.getApp(mu)
uHF = self.getHF(mu)
- return self.HSEngine.norm(uApp - uHF, normType)
+ return self.HFEngine.norm(uApp - uHF)
def plotHF(self, mu:complex, name : str = "u", **figspecs):
"""
Do some nice plots of the HF solution at arbitrary wavenumber.
Args:
mu: Target parameter.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
uHF = self.getHF(mu)
- self.HSEngine.plot(uHF, name = name, **figspecs)
+ self.HFEngine.plot(uHF, name = name, **figspecs)
def plotApp(self, mu:complex, name : str = "u", **figspecs):
"""
Do some nice plots of approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
uApp = self.getApp(mu)
- self.HSEngine.plot(uApp, name = name, **figspecs)
+ self.HFEngine.plot(uApp, name = name, **figspecs)
def plotErr(self, mu:complex, name : str = "u", **figspecs):
"""
Do some nice plots of approximation error at arbitrary wavenumber.
Args:
mu: Target parameter.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
uApp = self.getApp(mu)
uHF = self.getHF(mu)
- self.HSEngine.plot(uApp - uHF, name = name, **figspecs)
+ self.HFEngine.plot(uApp - uHF, name = name, **figspecs)
diff --git a/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py
index d3914d0..b584b68 100644
--- a/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py
+++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py
@@ -1,274 +1,255 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import copy
import warnings
import numpy as np
from rrompy.reduction_methods.lagrange.generic_approximant_lagrange import GenericApproximantLagrange
-from rrompy.utilities.base.types import Np1D, ListAny, DictAny, HSEng, HFEng
+from rrompy.utilities.base.types import Np1D, ListAny, DictAny, HFEng
from rrompy.utilities.base import purgeDict
__all__ = ['ApproximantLagrangePade']
class ApproximantLagrangePade(GenericApproximantLagrange):
"""
ROM Lagrange Pade' interpolant computation for parametric problems.
Args:
- HFEngine: HF problem solver. Should have members:
- - energyNormMatrix: assemble sparse matrix (in CSC format)
- associated to weighted H10 norm;
- - problemData: list of HF problem data (excluding mu);
- - setProblemData: set HF problem data and mu to prescribed values;
- - getLSBlocks: get algebraic system blocks.
- HSEngine: Hilbert space general purpose engine. Should have members:
- - norm: compute Hilbert norm of expression represented as complex
- numpy vector;
- - plot: plot Hilbert expression represented as complex numpy
- vector.
+ HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of samples current approximant relies upon;
defaults to 2;
- 'sampler': sample point generator; defaults to uniform sampler on
[0, 1];
- 'M': degree of Pade' interpolant numerator; defaults to 0;
- 'N': degree of Pade' interpolant denominator; defaults to 0.
Defaults to empty dict.
plotSnap(optional): What to plot of snapshots of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to [].
+ solution map. See plot method in HFEngine. Defaults to [].
plotSnapSpecs(optional): How to plot snapshots of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to {}.
+ solution map. See plot method in HFEngine. Defaults to {}.
Attributes:
HFEngine: HF problem solver.
- HSEngine: Hilbert space general purpose engine.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'M': degree of Pade' interpolant numerator;
- 'N': degree of Pade' interpolant denominator;
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
plotSnap: What to plot of snapshots of the Helmholtz solution map.
plotSnapSpecs: How to plot snapshots of the Helmholtz solution map.
solSnapshots: Array whose columns are FE dofs of snapshots of solution
map.
RPOD: Right factor of snapshots POD. If no POD, set to None.
POD: Whether to compute POD of snapshots.
- energyNormMatrix: Sparse matrix (in CSC format) assoociated to
- w-weighted H10 norm.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
- As: List of sparse matrices (in CSC format) representing coefficients
- of linear system matrix wrt mu.
- bs: List of numpy vectors representing coefficients of linear system
- RHS wrt mu.
"""
- def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0.,
+ def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
approxParameters : DictAny = {}, plotSnap : ListAny = [],
plotSnapSpecs : DictAny = {}):
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["M", "N"]
- super().__init__(HFEngine = HFEngine, HSEngine = HSEngine, mu0 = mu0,
+ super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
plotSnap = plotSnap, plotSnapSpecs = plotSnapSpecs)
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change M, N and S.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters, ["M", "N"],
True, True, baselevel = 1)
GenericApproximantLagrange.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
if "M" in keyList:
self.M = approxParameters["M"]
elif hasattr(self, "M"):
self.M = self.M
else:
self.M = 0
if "N" in keyList:
self.N = approxParameters["N"]
elif hasattr(self, "N"):
self.N = self.N
else:
self.N = 0
@property
def M(self):
"""Value of M. Its assignment may change S."""
return self._M
@M.setter
def M(self, M):
if M < 0: raise ArithmeticError("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
if hasattr(self, "S") and self.S < self.M + 1:
warnings.warn("Prescribed S is too small. Updating S to M + 1.",
stacklevel = 2)
self.S = self.M + 1
@property
def N(self):
"""Value of N. Its assignment may change S."""
return self._N
@N.setter
def N(self, N):
if N < 0: raise ArithmeticError("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
if hasattr(self, "S") and self.S < self.N + 1:
warnings.warn("Prescribed S is too small. Updating S to N + 1.",
stacklevel = 2)
self.S = self.N + 1
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
if S <= 0: raise ArithmeticError("S must be positive.")
if hasattr(self, "S"): Sold = self.S
else: Sold = -1
vals, label = [0] * 2, {0:"M", 1:"N"}
if hasattr(self, "M"): vals[0] = self.M
if hasattr(self, "N"): vals[1] = self.N
idxmax = np.argmax(vals)
if vals[idxmax] + 1 > S:
warnings.warn("Prescribed S is too small. Updating S to {} + 1."\
.format(label[idxmax]), stacklevel = 2)
self.S = vals[idxmax] + 1
else:
self._S = S
self._approxParameters["S"] = self.S
if Sold != self.S:
self.resetSamples()
- def approxRadius(self) -> float:
- """Sample radius."""
- return np.max(self.mu0 - self.mus)
-
def setupApprox(self):
"""
Compute Pade' interpolant.
SVD-based robust eigenvalue management.
"""
if not self.checkComputedApprox():
self.computeSnapshots()
TN = np.vander(self.radiusPade(self.mus), N = self.N + 1,
increasing = True)
if self.POD:
data = self.RPOD.T
else:
data = self.solSnapshots.T
RHSFull = np.empty((TN.shape[0], data.shape[1] * TN.shape[1]),
dtype = np.complex)
for j in range(TN.shape[0]):
RHSFull[j, :] = np.kron(TN[j, :], data[j, :])
G = np.polyfit(self.radiusPade(self.mus), RHSFull, deg = self.N,
w = self.ws)[0, :].reshape((TN.shape[1],
data.shape[1])).T
if self.POD:
_, _, V = np.linalg.svd(G, full_matrices = False)
self.Q = V[-1, :].conj()
else:
- G2 = G.conj().T.dot(self.energyNormMatrix.dot(G))
+ G2 = self.HFEngine.innerProduct(G, G)
_, eV = np.linalg.eigh(G2)
self.Q = eV[:, 0]
TNQ = TN.dot(self.Q).reshape((self.S, 1))
if self.POD:
data = self.solSnapshots.dot(self.RPOD)
else:
data = self.solSnapshots
RHS = np.multiply(data.T, TNQ)
self.P = np.polyfit(self.radiusPade(self.mus), RHS, deg = self.M,
w = self.ws)[::-1, :].T
self.approxParameters = {"N" : self.N, "M" : self.M, "S" : self.S}
self.lastApproxParameters = copy(self.approxParameters)
def radiusPade(self, mu:Np1D, mu0 : float = None) -> float:
"""
Compute translated radius to be plugged into Pade' approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Translated radius to be plugged into Pade' approximant.
"""
if mu0 is None: mu0 = self.mu0
- return (mu - mu0) / self.approxRadius()
+ return self.HFEngine.rescaling(mu) - self.HFEngine.rescaling(mu0)
def evalApprox(self, mu:complex):
"""
Evaluate Pade' approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
"""
self.setupApprox()
powerlist = np.power(self.radiusPade(mu),
range(max(self.M, self.N) + 1))
self.uApp = (self.P.dot(powerlist[:self.M + 1])
/ self.Q.dot(powerlist[:self.N + 1]))
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
- return np.roots(self.Q[::-1]) * self.approxRadius() + self.mu0
+ return self.HFEngine.rescalingInv(np.roots(self.Q[::-1])
+ + self.HFEngine.rescaling(self.mu0))
+
diff --git a/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py
index 8b146c4..7716612 100644
--- a/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py
+++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py
@@ -1,249 +1,244 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import copy
import warnings
import numpy as np
import scipy as sp
from rrompy.reduction_methods.lagrange.generic_approximant_lagrange import GenericApproximantLagrange
-from rrompy.utilities.base.types import Np1D, ListAny, DictAny, HFEng, HSEng
+from rrompy.utilities.base.types import Np1D, ListAny, DictAny, HFEng
from rrompy.utilities.base import purgeDict
__all__ = ['ApproximantLagrangeRB']
class ApproximantLagrangeRB(GenericApproximantLagrange):
"""
ROM RB approximant computation for parametric problems.
Args:
- HFEngine: HF problem solver. Should have members:
- - energyNormMatrix: assemble sparse matrix (in CSC format)
- associated to weighted H10 norm;
- - problemData: list of HF problem data (excluding mu);
- - setProblemData: set HF problem data and mu to prescribed values;
- - getLSBlocks: get algebraic system blocks;
- - liftDirichletData: perform lifting of Dirichlet BC as numpy
- complex vector.
- HSEngine: Hilbert space general purpose engine. Should have members:
- - norm: compute Hilbert norm of expression represented as complex
- numpy vector;
- - plot: plot Hilbert expression represented as complex numpy
- vector.
+ HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of samples current approximant relies upon;
defaults to 2;
- 'sampler': sample point generator; defaults to uniform sampler on
[0, 1];
- 'R': rank for Galerkin projection; defaults to S.
Defaults to empty dict.
plotSnap(optional): What to plot of snapshots of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to [].
+ solution map. See plot method in HFEngine. Defaults to [].
plotSnapSpecs(optional): How to plot snapshots of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to {}.
+ solution map. See plot method in HFEngine. Defaults to {}.
Attributes:
HFEngine: HF problem solver.
- HSEngine: Hilbert space general purpose engine.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths (unused).
approxRadius: Dummy radius of approximant (i.e. distance from mu0 to
farthest sample point).
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'R': rank for Galerkin projection.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
R: Rank for Galerkin projection.
plotSnap: What to plot of snapshots of the Helmholtz solution map.
plotSnapSpecs: How to plot snapshots of the Helmholtz solution map.
solSnapshots: Array whose columns are FE dofs of snapshots of solution
map.
RPOD: Right factor of snapshots POD. If no POD, set to None.
POD: Whether to compute POD of snapshots.
projMat: Projection matrix for RB system assembly.
- energyNormMatrix: Sparse matrix (in CSC format) assoociated to
- w-weighted H10 norm.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
As: List of sparse matrices (in CSC format) representing coefficients
- of linear system matrix wrt mu.
+ of linear system matrix wrt theta(mu).
bs: List of numpy vectors representing coefficients of linear system
- RHS wrt mu.
+ RHS wrt theta(mu).
+ thetaAs: List of callables representing coefficients of linear system
+ matrix wrt mu.
+ thetabs: List of callables representing coefficients of linear system
+ RHS wrt mu.
ARBs: List of sparse matrices (in CSC format) representing coefficients
- of compressed linear system matrix wrt mu.
+ of compressed linear system matrix wrt theta(mu).
bRBs: List of numpy vectors representing coefficients of compressed
- linear system RHS wrt mu.
+ linear system RHS wrt theta(mu).
"""
- def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0.,
+ def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
approxParameters : DictAny = {}, plotSnap : ListAny = [],
plotSnapSpecs : DictAny = {}):
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["R"]
- super().__init__(HFEngine = HFEngine, HSEngine = HSEngine, mu0 = mu0,
+ super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
plotSnap = plotSnap, plotSnapSpecs = plotSnapSpecs)
+ self.As, self.thetaAs = self.HFEngine.affineBlocksA(self.mu0)
+ self.bs, self.thetabs = self.HFEngine.affineBlocksb(self.mu0)
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self.projMat = None
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change M, N and S.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters, ["R"], True, True,
baselevel = 1)
GenericApproximantLagrange.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
if "R" in keyList:
self.R = approxParameters["R"]
elif hasattr(self, "R"):
self.R = self.R
else:
self.R = self.S
@property
def R(self):
"""Value of R. Its assignment may change S."""
return self._R
@R.setter
def R(self, R):
if R < 0: raise ArithmeticError("R must be non-negative.")
self._R = R
self._approxParameters["R"] = self.R
if hasattr(self, "S") and self.S < self.R:
warnings.warn("Prescribed S is too small. Updating S to R.",
stacklevel = 2)
self.S = self.R
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return (self.projMat is not None and super().checkComputedApprox())
def setupApprox(self):
"""Compute RB projection matrix."""
if not self.checkComputedApprox():
self.computeSnapshots()
if self.POD:
U, _, _ = np.linalg.svd(self.RPOD, full_matrices = False)
self.projMat = self.solSnapshots.dot(U[:, : self.R])
else:
self.projMat = self.solSnapshots[:, : self.R]
self.lastApproxParameters = copy(self.approxParameters)
self.assembleReducedSystem()
def assembleReducedSystem(self):
"""Build affine blocks of RB linear system through projections."""
self.setupApprox()
projMatH = self.projMat.T.conjugate()
self.ARBs = [None] * len(self.As)
self.bRBs = [None] * len(self.bs)
for j in range(len(self.As)):
self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat))
for j in range(len(self.bs)):
self.bRBs[j] = projMatH.dot(self.bs[j])
def solveReducedSystem(self, mu:complex) -> Np1D:
"""
Solve RB linear system.
Args:
mu: Target parameter.
Returns:
Solution of RB linear system.
"""
self.setupApprox()
- ARBmu = self.ARBs[0][:self.R,:self.R]
- bRBmu = self.bRBs[0][:self.R]
+ ARBmu = self.thetaAs[0](mu) * self.ARBs[0][:self.R,:self.R]
+ bRBmu = self.thetabs[0](mu) * self.bRBs[0][:self.R]
for j in range(1, len(self.ARBs)):
- ARBmu = ARBmu + np.power(mu, j) * self.ARBs[j][:self.R, :self.R]
+ ARBmu += self.thetaAs[j](mu) * self.ARBs[j][:self.R, :self.R]
for j in range(1, len(self.bRBs)):
- bRBmu = bRBmu + np.power(mu, j) * self.bRBs[j][:self.R]
+ bRBmu += self.thetabs[j](mu) * self.bRBs[j][:self.R]
return self.projMat[:, :self.R].dot(np.linalg.solve(ARBmu, bRBmu))
def evalApprox(self, mu:complex):
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
"""
self.setupApprox()
self.uApp = self.solveReducedSystem(mu)
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
+ warnings.warn(("Impossible to compute poles in general affine "
+ "parameter dependence. Results subject to "
+ "interpretation/rescaling, or possibly completely "
+ "wrong."), stacklevel = 2)
self.setupApprox()
A = np.eye(self.ARBs[0].shape[0] * (len(self.ARBs) - 1),
dtype = np.complex)
B = np.zeros_like(A)
A[: self.ARBs[0].shape[0], : self.ARBs[0].shape[0]] = - self.ARBs[0]
for j in range(len(self.ARBs) - 1):
Aj = self.ARBs[j + 1]
B[: Aj.shape[0], j * Aj.shape[0] : (j + 1) * Aj.shape[0]] = Aj
II = np.arange(self.ARBs[0].shape[0],
self.ARBs[0].shape[0] * (len(self.ARBs) - 1))
B[II, II - self.ARBs[0].shape[0]] = 1.
try:
return sp.linalg.eigvals(A, B)
except:
warnings.warn("Generalized eig algorithm did not converge.",
stacklevel = 2)
x = np.empty(A.shape[0])
x[:] = np.NaN
return x
diff --git a/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py
index b1ab829..deb9efe 100644
--- a/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py
+++ b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py
@@ -1,224 +1,207 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.reduction_methods.base.generic_approximant import GenericApproximant
-from rrompy.linalg.pod_engine import PODEngine
-from rrompy.utilities.base.types import Np1D, ListAny, DictAny, HFEng, HSEng
+from rrompy.sampling.scipy.pod_engine import PODEngine
+from rrompy.utilities.base.types import Np1D, ListAny, DictAny, HFEng
from rrompy.utilities.base import purgeDict
__all__ = ['GenericApproximantLagrange']
class GenericApproximantLagrange(GenericApproximant):
"""
ROM Lagrange interpolant computation for parametric problems (ABSTRACT).
Args:
- HFEngine: HF problem solver. Should have members:
- - energyNormMatrix: assemble sparse matrix (in CSC format)
- associated to weighted H10 norm;
- - problemData: list of HF problem data (excluding mu);
- - setProblemData: set HF problem data and k to prescribed values;
- - getLSBlocks: get algebraic system blocks.
- HSEngine: Hilbert space general purpose engine. Should have members:
- - norm: compute Hilbert norm of expression represented as complex
- numpy vector;
- - plot: plot Hilbert expression represented as complex numpy
- vector.
+ HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator; defaults to uniform sampler on
[0, 1].
Defaults to empty dict.
plotSnap(optional): What to plot of snapshots of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to [].
+ solution map. See plot method in HFEngine. Defaults to [].
plotSnapSpecs(optional): How to plot snapshots of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to {}.
+ solution map. See plot method in HFEngine. Defaults to {}.
Attributes:
HFEngine: HF problem solver.
- HSEngine: Hilbert space general purpose engine.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of snapshots current approximant relies upon;
- 'sampler': sample point generator.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
plotSnap: What to plot of snapshots of the Helmholtz solution map.
plotSnapSpecs: How to plot snapshots of the Helmholtz solution map.
solSnapshots: Array whose columns are FE dofs of snapshots of solution
map.
RPOD: Right factor of snapshots POD. If no POD, set to None.
POD: Whether to compute POD of snapshots.
- energyNormMatrix: Sparse matrix (in CSC format) assoociated to
- w-weighted H10 norm.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
- As: List of sparse matrices (in CSC format) representing coefficients
- of linear system matrix wrt mu.
- bs: List of numpy vectors representing coefficients of linear system
- RHS wrt mu.
"""
- def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0.,
+ def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
approxParameters : DictAny = {}, plotSnap : ListAny = [],
plotSnapSpecs : DictAny = {}):
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["S", "sampler"]
- super().__init__(HFEngine = HFEngine, HSEngine = HSEngine,
- mu0 = mu0, approxParameters = approxParameters)
+ super().__init__(HFEngine = HFEngine, mu0 = mu0,
+ approxParameters = approxParameters)
self.plotSnap = plotSnap
self.plotSnapSpecs = plotSnapSpecs
self.resetSamples()
@property
def mus(self):
"""Value of mus. Its assignment may reset snapshots."""
return self._mus
@mus.setter
def mus(self, mus):
musOld = self.mus if hasattr(self, 'mus') else None
self._mus = np.array(mus)
_, musCounts = np.unique(self._mus, return_counts = True)
if len(np.where(musCounts > 1)[0]) > 0:
raise Exception("Repeated sample points not allowed.")
if (musOld is None or len(self.mus) != len(musOld)
or not np.allclose(self.mus, musOld, 1e-14)):
self.resetSamples()
self.autoNode = None
@property
def approxParameters(self):
"""Value of approximant parameters. Its assignment may change S."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters, ["S", "sampler"],
True, True, baselevel = 1)
GenericApproximant.approxParameters.fset(self, approxParametersCopy)
keyList = list(approxParameters.keys())
if "S" in keyList:
self.S = approxParameters["S"]
elif hasattr(self, "S"):
self.S = self.S
else:
self.S = 2
if "sampler" in keyList:
self.sampler = approxParameters["sampler"]
elif not hasattr(self, "S"):
from rrompy.utilities.parameter_sampling import QuadratureSampler
self.sampler = QuadratureSampler([0., 1.], "UNIFORM")
del QuadratureSampler
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
if S <= 0: raise ArithmeticError("S must be positive.")
if hasattr(self, "S"): Sold = self.S
else: Sold = -1
self._S = S
self._approxParameters["S"] = self.S
if Sold != self.S:
self.resetSamples()
@property
def sampler(self):
"""Value of sampler."""
return self._sampler
@sampler.setter
def sampler(self, sampler):
if 'generatePoints' not in dir(sampler):
raise Exception("Sampler type not recognized.")
if hasattr(self, '_sampler'):
samplerOld = self.sampler
self._sampler = sampler
self._approxParameters["sampler"] = self.sampler
if not 'samplerOld' in locals() or samplerOld != self.sampler:
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
self.solSnapshots = None
self.RPOD = None
def computeSnapshots(self):
"""
Compute snapshots of solution map.
"""
if self.solSnapshots is None:
self.mus, self.ws = self.sampler.generatePoints(self.S)
self.mus = np.array([x[0] for x in self.mus])
for j, mu in enumerate(self.mus):
self.solveHF(mu)
- self.HSEngine.plot(self.uHF, name = "u({:.4f})".format(mu),
+ self.HFEngine.plot(self.uHF, name = "u({:.4f})".format(mu),
what = self.plotSnap, **self.plotSnapSpecs)
self.manageSnapshots(self.uHF, j)
def manageSnapshots(self, u:Np1D, pos:int):
"""
Store snapshots of solution map.
Args:
u: solution derivative as numpy complex vector;
pos: Derivative index.
"""
if pos == 0:
self.solSnapshots = np.empty((u.shape[0], self.S),
dtype = np.complex)
self.solSnapshots[:, pos] = u
if self.POD:
if pos == 0:
- self.PODEngine = PODEngine(self.energyNormMatrix)
+ self.PODEngine = PODEngine(self.HFEngine)
self.RPOD = np.eye(self.S, dtype = np.complex)
self.solSnapshots[:, pos], self.RPOD[: pos + 1, pos], _ = \
self.PODEngine.GS(self.solSnapshots[:, pos],
self.solSnapshots, pos)
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return self.solSnapshots is not None and super().checkComputedApprox()
diff --git a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
index 81d5565..a99da76 100644
--- a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
+++ b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
@@ -1,417 +1,422 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import copy
import warnings
import numpy as np
-from rrompy.reduction_methods.taylor.generic_approximant_taylor import GenericApproximantTaylor
-from rrompy.linalg.pod_engine import PODEngine
+from rrompy.reduction_methods.taylor.generic_approximant_taylor import (
+ GenericApproximantTaylor)
+from rrompy.sampling.scipy.pod_engine import PODEngine
from rrompy.utilities.base.types import Np1D, Np2D, ListAny, Tuple, DictAny
-from rrompy.utilities.base.types import HFEng, HSEng
+from rrompy.utilities.base.types import HFEng
from rrompy.utilities.base import purgeDict
__all__ = ['ApproximantTaylorPade']
class ApproximantTaylorPade(GenericApproximantTaylor):
"""
ROM single-point fast Pade' approximant computation for parametric
problems.
Args:
- HFEngine: HF problem solver. Should have members:
- - energyNormMatrix: assemble sparse matrix (in CSC format)
- associated to weighted H10 norm;
- - problemData: list of HF problem data (excluding mu);
- - setProblemData: set HF problem data and mu to prescribed values;
- - setupDerivativeComputation: setup HF problem data to compute j-th
- solution derivative at mu0;
- - solve: get HF solution as complex numpy vector.
- HSEngine: Hilbert space general purpose engine. Should have members:
- - norm: compute Hilbert norm of expression represented as complex
- numpy vector;
- - plot: plot Hilbert expression represented as complex numpy vector.
+ HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'rho': weight for computation of original Pade' approximant;
defaults to np.inf, i.e. fast approximant;
- 'M': degree of Pade' approximant numerator; defaults to 0;
- 'N': degree of Pade' approximant denominator; defaults to 0;
- 'E': total number of derivatives current approximant relies upon;
defaults to Emax;
- 'Emax': total number of derivatives of solution map to be
computed; defaults to E;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0;
- 'sampleType': label of sampling type; available values are:
- 'ARNOLDI': orthogonalization of solution derivatives through
Arnoldi algorithm;
- 'KRYLOV': standard computation of solution derivatives.
Defaults to 'KRYLOV'.
Defaults to empty dict.
plotDer(optional): What to plot of derivatives of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to [].
+ solution map. See plot method in HFEngine. Defaults to [].
plotDerSpecs(optional): How to plot derivatives of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to {}.
+ solution map. See plot method in HFEngine. Defaults to {}.
Attributes:
HFEngine: HF problem solver.
- HSEngine: Hilbert space general purpose engine.
mu0: Default parameter.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'rho': weight for computation of original Pade' approximant;
- 'M': degree of Pade' approximant numerator;
- 'N': degree of Pade' approximant denominator;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'robustTol': tolerance for robust Pade' denominator management;
- 'sampleType': label of sampling type.
POD: Whether to compute QR factorization of derivatives.
rho: Weight of approximant.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
robustTol: tolerance for robust Pade' denominator management.
sampleType: Label of sampling type.
plotDer: What to plot of derivatives of the Helmholtz solution map.
plotDerSpecs: How to plot derivatives of the Helmholtz solution map.
initialHFData: HF problem initial data.
- energyNormMatrix: Sparse matrix (in CSC format) assoociated to
- w-weighted H10 norm.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
G: Square Numpy 2D vector of size (N+1) corresponding to Pade'
denominator matrix (see paper).
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
- def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0,
+ def __init__(self, HFEngine:HFEng, mu0 : complex = 0,
approxParameters : DictAny = {}, plotDer : ListAny = [],
plotDerSpecs : DictAny = {}):
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["M", "N", "robustTol", "rho"]
- super().__init__(HFEngine = HFEngine, HSEngine = HSEngine, mu0 = mu0,
+ super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
plotDer = plotDer, plotDerSpecs = plotDerSpecs)
@property
def approxParameters(self):
"""Value of approximant parameters."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters,
["M", "N", "robustTol", "rho"],
True, True, baselevel = 1)
keyList = list(approxParameters.keys())
if "rho" in keyList:
self._rho = approxParameters["rho"]
elif hasattr(self, "rho"):
self._rho = self.rho
else:
self._rho = np.inf
GenericApproximantTaylor.approxParameters.fset(self,
approxParametersCopy)
self.rho = self._rho
if "robustTol" in keyList:
self.robustTol = approxParameters["robustTol"]
elif hasattr(self, "robustTol"):
self.robustTol = self.robustTol
else:
self.robustTol = 0
self._ignoreParWarnings = True
if "M" in keyList:
self.M = approxParameters["M"]
elif hasattr(self, "M"):
self.M = self.M
else:
self.M = 0
del self._ignoreParWarnings
if "N" in keyList:
self.N = approxParameters["N"]
elif hasattr(self, "N"):
self.N = self.N
else:
self.N = 0
@property
def rho(self):
"""Value of rho."""
return self._rho
@rho.setter
def rho(self, rho):
self._rho = np.abs(rho)
self._approxParameters["rho"] = self.rho
@property
def M(self):
"""Value of M. Its assignment may change Emax and E."""
return self._M
@M.setter
def M(self, M):
if M < 0: raise ArithmeticError("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
if not hasattr(self, "_ignoreParWarnings"):
self.checkMNEEmax()
@property
def N(self):
"""Value of N. Its assignment may change Emax."""
return self._N
@N.setter
def N(self, N):
if N < 0: raise ArithmeticError("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
self.checkMNEEmax()
def checkMNEEmax(self):
"""Check consistency of M, N, E, and Emax."""
M = self.M if hasattr(self, "_M") else 0
N = self.N if hasattr(self, "_N") else 0
E = self.E if hasattr(self, "_E") else M + N
Emax = self.Emax if hasattr(self, "_Emax") else M + N
if self.rho == np.inf:
if Emax < max(N, M):
warnings.warn(("Prescribed Emax is too small. Updating Emax "
"to max(M, N)."), stacklevel = 3)
self.Emax = max(N, M)
if E < max(N, M):
warnings.warn(("Prescribed E is too small. Updating E to "
"max(M, N)."), stacklevel = 3)
self.E = max(N, M)
else:
if Emax < N + M:
warnings.warn(("Prescribed Emax is too small. Updating "
"Emax to M + N."), stacklevel = 3)
self.Emax = self.N + M
if E < N + M:
warnings.warn(("Prescribed E is too small. Updating E to "
"M + N."), stacklevel = 3)
self.E = self.N + M
@property
def robustTol(self):
"""Value of tolerance for robust Pade' denominator management."""
return self._robustTol
@robustTol.setter
def robustTol(self, robustTol):
if robustTol < 0.:
warnings.warn(("Overriding prescribed negative robustness "
"tolerance to 0."), stacklevel = 2)
robustTol = 0.
self._robustTol = robustTol
self._approxParameters["robustTol"] = self.robustTol
@property
def E(self):
"""Value of E. Its assignment may change Emax."""
return self._E
@E.setter
def E(self, E):
if E < 0: raise ArithmeticError("E must be non-negative.")
self._E = E
self.checkMNEEmax()
self._approxParameters["E"] = self.E
if hasattr(self, "Emax") and self.Emax < self.E:
warnings.warn(("Prescribed Emax is too small. Updating Emax "
"to E."), stacklevel = 2)
self.Emax = self.E
@property
def Emax(self):
"""Value of Emax. Its assignment may reset computed derivatives."""
return self._Emax
@Emax.setter
def Emax(self, Emax):
if Emax < 0: raise ArithmeticError("Emax must be non-negative.")
if hasattr(self, "Emax"): EmaxOld = self.Emax
else: EmaxOld = -1
if hasattr(self, "_N"): N = self.N
else: N = 0
if hasattr(self, "_M"): M = self.M
else: M = 0
if hasattr(self, "_E"): E = self.E
else: E = 0
if self.rho == np.inf:
if max(N, M, E) > Emax:
warnings.warn(("Prescribed Emax is too small. Updating Emax "
"to max(N, M, E)."), stacklevel = 2)
Emax = max(N, M, E)
else:
if max(N + M, E) > Emax:
warnings.warn(("Prescribed Emax is too small. Updating Emax "
"to max(N + M, E)."), stacklevel = 2)
Emax = max(N + M, E)
self._Emax = Emax
self._approxParameters["Emax"] = Emax
- if EmaxOld >= self.Emax and self.LAEngine.samples is not None:
- self.LAEngine.samples = self.LAEngine.samples[:,
+ if EmaxOld >= self.Emax and self.samplingEngine.samples is not None:
+ self.samplingEngine.samples = self.samplingEngine.samples[:,
: self.Emax + 1]
if (self.sampleType == "ARNOLDI"
- and self.LAEngine.HArnoldi is not None):
- self.LAEngine.HArnoldi = self.LAEngine.HArnoldi[
+ and self.samplingEngine.HArnoldi is not None):
+ self.samplingEngine.HArnoldi = self.samplingEngine.HArnoldi[
: self.Emax + 1,
: self.Emax + 1]
- self.LAEngine.RArnoldi = self.LAEngine.RArnoldi[
+ self.samplingEngine.RArnoldi = self.samplingEngine.RArnoldi[
: self.Emax + 1,
: self.Emax + 1]
def setupApprox(self):
"""
Compute Pade' approximant. SVD-based robust eigenvalue management.
"""
if not self.checkComputedApprox():
self.computeDerivatives()
while True:
if self.POD:
ev, eV = self.findeveVGQR()
else:
ev, eV = self.findeveVGExplicit()
ts = self.robustTol * np.linalg.norm(ev)
Nnew = np.sum(np.abs(ev) >= ts)
diff = self.N - Nnew
if diff <= 0: break
Enew = self.E - diff
Mnew = min(self.M, Enew)
if Mnew == self.M:
strM = ""
else:
strM = ", M from {0} to {1},".format(self.M, Mnew)
print(("Smallest {0} eigenvalues below tolerance.\n"
"Reducing N from {1} to {2}{5} and E from {3} to {4}.")\
.format(diff + 1, self.N, Nnew, self.E, Enew, strM))
newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew}
self.approxParameters = newParameters
self.Q = eV[::-1, 0]
QToeplitz = np.zeros((self.Emax + 1, self.M + 1),
dtype = np.complex)
for i in range(len(self.Q)):
rng = np.arange(self.M + 1 - i)
QToeplitz[rng, rng + i] = self.Q[i]
if self.sampleType == "ARNOLDI":
- QToeplitz = self.LAEngine.RArnoldi.dot(QToeplitz)
- self.P = self.LAEngine.samples.dot(QToeplitz)
+ QToeplitz = self.samplingEngine.RArnoldi.dot(QToeplitz)
+ self.P = self.samplingEngine.samples.dot(QToeplitz)
self.lastApproxParameters = copy(self.approxParameters)
def buildG(self):
"""Assemble Pade' denominator matrix."""
self.computeDerivatives()
if self.rho == np.inf:
Nmin = self.E - self.N
else:
Nmin = self.M - self.N + 1
if self.sampleType == "KRYLOV":
- DerE = self.LAEngine.samples[:, Nmin : self.E + 1]
- G = DerE.conj().T.dot(self.energyNormMatrix.dot(DerE))
+ DerE = self.samplingEngine.samples[:, Nmin : self.E + 1]
+ G = self.HFEngine.innerProduct(DerE, DerE)
else:
- RArnE = self.LAEngine.RArnoldi[: self.E + 1, Nmin : self.E + 1]
+ RArnE = self.samplingEngine.RArnoldi[: self.E + 1,
+ Nmin : self.E + 1]
G = RArnE.conj().T.dot(RArnE)
if self.rho == np.inf:
self.G = G
else:
Gbig = G
- self.G = np.zeros((self.N + 1, self.N + 1), dtype=np.complex)
+ self.G = np.zeros((self.N + 1, self.N + 1), dtype = np.complex)
for k in range(self.E - self.M):
self.G += self.rho ** (2 * k) * Gbig[k : k + self.N + 1,
k : k + self.N + 1]
def findeveVGExplicit(self) -> Tuple[Np1D, Np2D]:
"""
Compute explicitly eigenvalues and eigenvectors of Pade' denominator
matrix.
"""
self.buildG()
ev, eV = np.linalg.eigh(self.G)
return ev, eV
def findeveVGQR(self) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of Pade' denominator matrix
through SVD of R factor. See ``Householder triangularization of a
quasimatrix'', L.Trefethen, 2008 for QR algorithm.
Returns:
Eigenvalues in ascending order and corresponding eigenvector
matrix.
"""
self.computeDerivatives()
if self.rho == np.inf:
Nmin = self.E - self.N
else:
Nmin = self.M - self.N + 1
if self.sampleType == "KRYLOV":
- A = copy(self.LAEngine.samples[:, Nmin : self.E + 1])
- self.PODEngine = PODEngine(self.energyNormMatrix)
+ A = copy(self.samplingEngine.samples[:, Nmin : self.E + 1])
+ self.PODEngine = PODEngine(self.HFEngine)
R = self.PODEngine.QRHouseholder(A, only_R = True)
else:
- R = self.LAEngine.RArnoldi[: self.E + 1, Nmin : self.E + 1]
+ R = self.samplingEngine.RArnoldi[: self.E + 1, Nmin : self.E + 1]
if self.rho == np.inf:
_, s, V = np.linalg.svd(R, full_matrices = False)
else:
Rtower = np.zeros((R.shape[0] * (self.E - self.M), self.N + 1),
dtype = np.complex)
for k in range(self.E - self.M):
Rtower[k * R.shape[0] : (k + 1) * R.shape[0], :] = (
self.rho ** k
* R[:, self.M - self.N + 1 + k : self.M + 1 + k])
_, s, V = np.linalg.svd(Rtower, full_matrices = False)
return s[::-1], V.conj().T[:, ::-1]
+ def radiusPade(self, mu:Np1D, mu0 : float = None) -> float:
+ """
+ Compute translated radius to be plugged into Pade' approximant.
+
+ Args:
+ mu: Parameter(s) 1.
+ mu0: Parameter(s) 2. If None, set to self.mu0.
+
+ Returns:
+ Translated radius to be plugged into Pade' approximant.
+ """
+ if mu0 is None: mu0 = self.mu0
+ return self.HFEngine.rescaling(mu) - self.HFEngine.rescaling(mu0)
+
def evalApprox(self, mu:complex):
"""
Evaluate Pade' approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
"""
self.setupApprox()
- powerlist = np.power(mu - self.mu0, range(max(self.M, self.N) + 1))
+ powerlist = np.power(self.radiusPade(mu),
+ range(max(self.M, self.N) + 1))
self.uApp = (self.P.dot(powerlist[:self.M + 1])
/ self.Q.dot(powerlist[:self.N + 1]))
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
- return np.roots(self.Q[::-1]) + self.mu0
+ return self.HFEngine.rescalingInv(np.roots(self.Q[::-1])
+ + self.HFEngine.rescaling(self.mu0))
+
diff --git a/rrompy/reduction_methods/taylor/approximant_taylor_rb.py b/rrompy/reduction_methods/taylor/approximant_taylor_rb.py
index a3b28ef..decdd1d 100644
--- a/rrompy/reduction_methods/taylor/approximant_taylor_rb.py
+++ b/rrompy/reduction_methods/taylor/approximant_taylor_rb.py
@@ -1,277 +1,266 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import copy
import warnings
import numpy as np
import scipy as sp
from rrompy.reduction_methods.taylor.generic_approximant_taylor import GenericApproximantTaylor
-from rrompy.linalg.pod_engine import PODEngine
-from rrompy.utilities.base.types import Np1D, ListAny, DictAny, HFEng, HSEng
+from rrompy.sampling.scipy.pod_engine import PODEngine
+from rrompy.utilities.base.types import Np1D, ListAny, DictAny, HFEng
from rrompy.utilities.base import purgeDict
__all__ = ['ApproximantTaylorRB']
class ApproximantTaylorRB(GenericApproximantTaylor):
"""
ROM single-point fast RB approximant computation for parametric problems
with polynomial dependence up to degree 2.
Args:
- HFEngine: HF problem solver. Should have members:
- - energyNormMatrix: sparse matrix (in CSC format) associated to
- w-weighted H10 norm;
- - problemData: list of HF problem data (excluding mu);
- - setProblemData: set HF problem data and mu to prescribed values;
- - getLSBlocks: get blocks of HF linear system as sparse matrices in
- CSC format;
- - liftDirichletData: perform lifting of Dirichlet BC as numpy
- complex vector;
- - setupDerivativeComputation: setup HF problem data to compute j-th
- solution derivative at mu0;
- - solve: get HF solution as complex numpy vector.
- HSEngine: Hilbert space general purpose engine. Should have members:
- - norm: compute Hilbert norm of expression represented as complex
- numpy vector;
- - plot: plot Hilbert expression represented as complex numpy vector.
+ HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'R': rank for Galerkin projection; defaults to E + 1;
- 'E': total number of derivatives current approximant relies upon;
defaults to Emax;
- 'Emax': total number of derivatives of solution map to be
computed; defaults to E;
- 'sampleType': label of sampling type; available values are:
- 'ARNOLDI': orthogonalization of solution derivatives through
Arnoldi algorithm;
- 'KRYLOV': standard computation of solution derivatives.
Defaults to 'KRYLOV'.
Defaults to empty dict.
plotDer(optional): What to plot of derivatives of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to [].
+ solution map. See plot method in HFEngine. Defaults to [].
plotDerSpecs(optional): How to plot derivatives of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to {}.
+ solution map. See plot method in HFEngine. Defaults to {}.
Attributes:
HFEngine: HF problem solver.
- HSEngine: Hilbert space general purpose engine.
mu0: Default parameter.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'R': rank for Galerkin projection;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'sampleType': label of sampling type.
POD: Whether to compute QR factorization of derivatives.
R: Rank for Galerkin projection.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
sampleType: Label of sampling type, i.e. 'KRYLOV'.
plotDer: What to plot of derivatives of the Helmholtz solution map.
plotDerSpecs: How to plot derivatives of the Helmholtz solution map.
- energyNormMatrix: Sparse matrix (in CSC format) associated to
- w-weighted H10 norm.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
solLifting: Numpy complex vector with lifting of real part of Dirichlet
boundary datum.
projMat: Numpy matrix representing projection onto RB space.
projMat: Numpy matrix representing projection onto RB space.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix wrt mu.
bs: List of numpy vectors representing coefficients of linear system
RHS wrt mu.
ARBs: List of sparse matrices (in CSC format) representing RB
coefficients of linear system matrix wrt mu.
bRBs: List of numpy vectors representing RB coefficients of linear
system RHS wrt mu.
"""
- def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0,
+ def __init__(self, HFEngine:HFEng, mu0 : complex = 0,
approxParameters : DictAny = {}, plotDer : ListAny = [],
plotDerSpecs : DictAny = {}):
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["R"]
- super().__init__(HFEngine = HFEngine, HSEngine = HSEngine, mu0 = mu0,
+ super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
plotDer = plotDer, plotDerSpecs = plotDerSpecs)
+ self.As, self.thetaAs = self.HFEngine.affineBlocksA(self.mu0)
+ self.bs, self.thetabs = self.HFEngine.affineBlocksb(self.mu0)
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self.projMat = None
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change M, N and S.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters, ["R"],
True, True, baselevel = 1)
GenericApproximantTaylor.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
if "R" in keyList:
self.R = approxParameters["R"]
elif hasattr(self, "R"):
self.R = self.R
else:
self.R = self.E + 1
@property
def POD(self):
"""Value of POD."""
return self._POD
@POD.setter
def POD(self, POD):
GenericApproximantTaylor.POD.fset(self, POD)
if (hasattr(self, "sampleType") and self.sampleType == "ARNOLDI"
and not self.POD):
warnings.warn(("Arnoldi sampling implicitly forces POD-type "
"derivative management."), stacklevel = 2)
@property
def sampleType(self):
"""Value of sampleType."""
return self._sampleType
@sampleType.setter
def sampleType(self, sampleType):
GenericApproximantTaylor.sampleType.fset(self, sampleType)
if (hasattr(self, "POD") and not self.POD
and self.sampleType == "ARNOLDI"):
warnings.warn(("Arnoldi sampling implicitly forces POD-type "
"derivative management."), stacklevel = 2)
@property
def R(self):
"""Value of R. Its assignment may change S."""
return self._R
@R.setter
def R(self, R):
if R < 0: raise ArithmeticError("R must be non-negative.")
self._R = R
self._approxParameters["R"] = self.R
if hasattr(self, "E") and self.E + 1 < self.R:
warnings.warn("Prescribed E is too small. Updating E to R - 1.",
stacklevel = 2)
self.E = self.R - 1
def setupApprox(self):
"""Setup RB system."""
if not self.checkComputedApprox():
self.computeDerivatives()
if self.POD and not self.sampleType == "ARNOLDI":
- self.PODEngine = PODEngine(self.energyNormMatrix)
+ self.PODEngine = PODEngine(self.HFEngine)
self.projMatQ, self.projMatR = self.PODEngine.QRHouseholder(
- self.LAEngine.samples)
+ self.samplingEngine.samples)
if self.POD:
if self.sampleType == "ARNOLDI":
- self.projMatR = self.LAEngine.RArnoldi
- self.projMatQ = self.LAEngine.samples
+ self.projMatR = self.samplingEngine.RArnoldi
+ self.projMatQ = self.samplingEngine.samples
U, _, _ = np.linalg.svd(self.projMatR[: self.E + 1, : self.E + 1])
self.projMat = self.projMatQ[:, : self.E + 1].dot(U[:, : self.R])
else:
- self.projMat = self.LAEngine.samples[:, : self.R + 1]
+ self.projMat = self.samplingEngine.samples[:, : self.R + 1]
self.lastApproxParameters = copy(self.approxParameters)
self.assembleReducedSystem()
def assembleReducedSystem(self):
"""Build affine blocks of RB linear system through projections."""
if not self.checkComputedApprox():
self.setupApprox()
projMatH = self.projMat.T.conj()
self.ARBs = [None] * len(self.As)
self.bRBs = [None] * len(self.bs)
for j in range(len(self.As)):
self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat))
for j in range(len(self.bs)):
self.bRBs[j] = projMatH.dot(self.bs[j])
def solveReducedSystem(self, mu:complex) -> Np1D:
"""
Solve RB linear system.
Args:
mu: Target parameter.
Returns:
Solution of RB linear system.
"""
self.setupApprox()
- ARBmu = self.ARBs[0][: self.R, : self.R]
- bRBmu = self.bRBs[0][: self.R]
+ ARBmu = self.thetaAs[0](mu) * self.ARBs[0][:self.R,:self.R]
+ bRBmu = self.thetabs[0](mu) * self.bRBs[0][:self.R]
for j in range(1, len(self.ARBs)):
- ARBmu = ARBmu + np.power(mu, j) * self.ARBs[j][:self.R, :self.R]
+ ARBmu += self.thetaAs[j](mu) * self.ARBs[j][:self.R, :self.R]
for j in range(1, len(self.bRBs)):
- bRBmu = bRBmu + np.power(mu, j) * self.bRBs[j][:self.R]
+ bRBmu += self.thetabs[j](mu) * self.bRBs[j][:self.R]
return self.projMat[:, :self.R].dot(np.linalg.solve(ARBmu, bRBmu))
def evalApprox(self, mu:complex):
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
"""
self.setupApprox()
self.uApp = self.solveReducedSystem(mu)
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
+ warnings.warn(("Impossible to compute poles in general affine "
+ "parameter dependence. Results subject to "
+ "interpretation/rescaling, or possibly completely "
+ "wrong."), stacklevel = 2)
self.setupApprox()
A = np.eye(self.ARBs[0].shape[0] * (len(self.ARBs) - 1),
dtype = np.complex)
B = np.zeros_like(A)
A[: self.ARBs[0].shape[0], : self.ARBs[0].shape[0]] = - self.ARBs[0]
for j in range(len(self.ARBs) - 1):
Aj = self.ARBs[j + 1]
B[: Aj.shape[0], j * Aj.shape[0] : (j + 1) * Aj.shape[0]] = Aj
II = np.arange(self.ARBs[0].shape[0],
self.ARBs[0].shape[0] * (len(self.ARBs) - 1))
B[II, II - self.ARBs[0].shape[0]] = 1.
try:
return sp.linalg.eigvals(A, B)
except:
warnings.warn("Generalized eig algorithm did not converge.",
stacklevel = 2)
x = np.empty(A.shape[0])
x[:] = np.NaN
return x
+
diff --git a/rrompy/reduction_methods/taylor/generic_approximant_taylor.py b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py
index 88faaeb..190d37e 100644
--- a/rrompy/reduction_methods/taylor/generic_approximant_taylor.py
+++ b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py
@@ -1,252 +1,231 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import warnings
-import numpy as np
-from scipy.special import comb
-import scipy.sparse.linalg as spla
from rrompy.reduction_methods.base.generic_approximant import GenericApproximant
-from rrompy.linalg.pod_engine import PODEngine
-from rrompy.utilities.base.types import Np1D, ListAny, DictAny, HFEng, HSEng, LAEng
+from rrompy.utilities.base.types import ListAny, DictAny, HFEng
from rrompy.utilities.base import purgeDict
__all__ = ['GenericApproximantTaylor']
class GenericApproximantTaylor(GenericApproximant):
"""
ROM single-point approximant computation for parametric problems
(ABSTRACT).
Args:
- HFEngine: HF problem solver. Should have members:
- - energyNormMatrix: assemble sparse matrix (in CSC format)
- associated to weighted H10 norm;
- - problemData: list of HF problem data (excluding mu);
- - setProblemData: set HF problem data and mu to prescribed values;
- - getLSBlocks: get algebraic system blocks.
- HSEngine: Hilbert space general purpose engine. Should have members:
- - norm: compute Hilbert norm of expression represented as complex
- numpy vector;
- - plot: plot Hilbert expression represented as complex numpy vector.
+ HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'E': total number of derivatives current approximant relies upon;
defaults to Emax;
- 'Emax': total number of derivatives of solution map to be
computed; defaults to E;
- 'sampleType': label of sampling type; available values are:
- 'ARNOLDI': orthogonalization of solution derivatives through
Arnoldi algorithm;
- 'KRYLOV': standard computation of solution derivatives.
Defaults to 'KRYLOV'.
Defaults to empty dict.
plotDer(optional): What to plot of derivatives of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to [].
+ solution map. See plot method in HFEngine. Defaults to [].
plotDerSpecs(optional): How to plot derivatives of the Helmholtz
- solution map. See plot method in HSEngine. Defaults to {}.
+ solution map. See plot method in HFEngine. Defaults to {}.
Attributes:
HFEngine: HF problem solver.
- HSEngine: Hilbert space general purpose engine.
mu0: Default parameter.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'sampleType': label of sampling type.
POD: Whether to compute QR factorization of derivatives.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
sampleType: Label of sampling type.
plotDer: What to plot of derivatives of the Helmholtz solution map.
plotDerSpecs: How to plot derivatives of the Helmholtz solution map.
initialHFData: HF problem initial data.
- energyNormMatrix: Sparse matrix (in CSC format) assoociated to
- w-weighted H10 norm.
- LAEngine: Linear algebra engine.
+ samplingEngine: Sampling engine.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
solLifting: Lifting of Dirichlet boundary data as numpy vector.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
- As: List of sparse matrices (in CSC format) representing coefficients
- of linear system matrix wrt mu.
- bs: List of numpy vectors representing coefficients of linear system
- RHS wrt mu.
"""
- def __init__(self, HFEngine:HFEng, HSEngine:HSEng, mu0 : complex = 0,
+ def __init__(self, HFEngine:HFEng, mu0 : complex = 0,
approxParameters : DictAny = {}, plotDer : ListAny = [],
plotDerSpecs : DictAny = {}):
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["E", "Emax", "sampleType"]
- super().__init__(HFEngine = HFEngine, HSEngine = HSEngine, mu0 = mu0,
+ super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters)
self.plotDer = plotDer
self.plotDerSpecs = plotDerSpecs
self.resetSamples()
- def setupLAEngine(self) -> LAEng:
- """Setup linear algebra engine."""
+ def setupSampling(self):
+ """Setup sampling engine."""
if self.sampleType == "ARNOLDI":
- from rrompy.linalg.linalg_arnoldi_engine import LinAlgArnoldiEngine
- self.LAEngine = LinAlgArnoldiEngine(self.As, self.bs,
- self.energyNormMatrix,
- mu0 = self.HFEngine.mu0Blocks)
+ from rrompy.sampling.scipy.sampling_engine_arnoldi import (
+ SamplingEngineArnoldi)
+ self.samplingEngine = SamplingEngineArnoldi(self.HFEngine)
elif self.sampleType == "KRYLOV":
- from rrompy.linalg.linalg_krylov_engine import LinAlgKrylovEngine
- self.LAEngine = LinAlgKrylovEngine(self.As, self.bs,
- mu0 = self.HFEngine.mu0Blocks)
+ from rrompy.sampling.scipy.sampling_engine_krylov import (
+ SamplingEngineKrylov)
+ self.samplingEngine = SamplingEngineKrylov(self.HFEngine)
else:
raise Exception("Sample type not recognized.")
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change E and Emax.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters,
["E", "Emax", "sampleType"],
True, True, baselevel = 1)
GenericApproximant.approxParameters.fset(self, approxParametersCopy)
keyList = list(approxParameters.keys())
if "E" in keyList:
self._E = approxParameters["E"]
self._approxParameters["E"] = self.E
if "Emax" in keyList:
self.Emax = approxParameters["Emax"]
else:
if not hasattr(self, "Emax"):
self.Emax = self.E
else:
self.Emax = self.Emax
else:
if "Emax" in keyList:
self._E = approxParameters["Emax"]
self._approxParameters["E"] = self.E
self.Emax = self.E
else:
if not (hasattr(self, "Emax") and hasattr(self, "E")):
raise Exception("At least one of E and Emax must be set.")
if "sampleType" in keyList:
self.sampleType = approxParameters["sampleType"]
elif hasattr(self, "sampleType"):
self.sampleType = self.sampleType
else:
self.sampleType = "KRYLOV"
@property
def E(self):
"""Value of E. Its assignment may change Emax."""
return self._E
@E.setter
def E(self, E):
if E < 0: raise ArithmeticError("E must be non-negative.")
self._E = E
self._approxParameters["E"] = self.E
if hasattr(self, "Emax") and self.Emax < self.E:
warnings.warn("Prescribed E is too large. Updating Emax to E.",
stacklevel = 2)
self.Emax = self.E
@property
def Emax(self):
"""Value of Emax. Its assignment may reset computed derivatives."""
return self._Emax
@Emax.setter
def Emax(self, Emax):
if Emax < 0: raise ArithmeticError("Emax must be non-negative.")
if hasattr(self, "Emax"): EmaxOld = self.Emax
else: EmaxOld = -1
self._Emax = Emax
if hasattr(self, "E") and self.Emax < self.E:
warnings.warn("Prescribed Emax is too small. Updating Emax to E.",
stacklevel = 2)
self.Emax = self.E
else:
self._approxParameters["Emax"] = self.Emax
- if EmaxOld >= self.Emax and self.LAEngine.samples is not None:
- self.LAEngine.samples = self.LAEngine.samples[:,
- : self.Emax + 1]
+ if EmaxOld >= self.Emax and self.samplingEngine.samples is not None:
+ self.samplingEngine.samples = self.samplingEngine.samples[:,
+ : self.Emax + 1]
if (self.sampleType == "ARNOLDI"
- and self.LAEngine.HArnoldi is not None):
- self.LAEngine.HArnoldi = self.LAEngine.HArnoldi[
+ and self.samplingEngine.HArnoldi is not None):
+ self.samplingEngine.HArnoldi= self.samplingEngine.HArnoldi[
: self.Emax + 1,
: self.Emax + 1]
- self.LAEngine.RArnoldi = self.LAEngine.RArnoldi[
+ self.samplingEngine.RArnoldi= self.samplingEngine.RArnoldi[
: self.Emax + 1,
: self.Emax + 1]
else:
self.resetSamples()
@property
def sampleType(self):
"""Value of sampleType."""
return self._sampleType
@sampleType.setter
def sampleType(self, sampleType):
if hasattr(self, "sampleType"): sampleTypeOld = self.sampleType
else: sampleTypeOld = -1
try:
sampleType = sampleType.upper().strip().replace(" ","")
if sampleType not in ["ARNOLDI", "KRYLOV"]:
raise Exception("Sample type not recognized.")
self._sampleType = sampleType
except:
warnings.warn(("Prescribed sampleType not recognized. Overriding "
"to 'KRYLOV'."), stacklevel = 2)
self._sampleType = "KRYLOV"
self._approxParameters["sampleType"] = self.sampleType
if sampleTypeOld != self.sampleType:
self.resetSamples()
def computeDerivatives(self):
"""Compute derivatives of solution map starting from order 0."""
- if self.LAEngine.samples is None:
- self.LAEngine.iterSample(self.mu0, self.Emax + 1)
+ if self.samplingEngine.samples is None:
+ self.samplingEngine.iterSample(self.mu0, self.Emax + 1)
for j in range(self.Emax + 1):
- self.HSEngine.plot(self.LAEngine.samples[:, j],
+ self.HFEngine.plot(self.samplingEngine.samples[:, j],
name = "u_{0}".format(j),
what = self.plotDer, **self.plotDerSpecs)
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
- return (self.LAEngine.samples is not None
+ return (self.samplingEngine.samples is not None
and super().checkComputedApprox())
diff --git a/rrompy/utilities/base/types.py b/rrompy/utilities/base/types.py
index 40c14d7..a0f8497 100644
--- a/rrompy/utilities/base/types.py
+++ b/rrompy/utilities/base/types.py
@@ -1,50 +1,53 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from typing import TypeVar, List, Tuple, Dict, Any
-__all__ = ['TupleAny','ListAny','DictAny','Np1D','Np2D','Np1DLst','N2FSExpr',
- 'FenSpace','FenExpr','FenFunc','HFEng','HSEng','ROMEng','GenExpr',
- 'strLst','BfSExpr']
+__all__ = ['TupleAny','ListAny','DictAny','HS1D','HS2D','HSOp','Np1D','Np2D',
+ 'Np1DLst','N2FSExpr','FenSpace','FenExpr','FenFunc','HFEng',
+ 'ROMEng','GenExpr','strLst','BfSExpr']
# ANY
TupleAny = Tuple[Any]
ListAny = List[Any]
DictAny = Dict[Any, Any]
+# GENERIC
+HS1D = TypeVar("Hilbert space element")
+HS2D = TypeVar("Hilbert space quasi-matrix")
+HSOp = TypeVar("Hilbert space operator")
+
# NUMPY
Np1D = TypeVar("NumPy 1D array")
Np2D = TypeVar("NumPy 2D array")
Np1DLst = TypeVar("NumPy 1D array or list of NumPy 1D array")
N2FSExpr = TypeVar("NumPy 2D array, float or str")
# FENICS
FenSpace = TypeVar("FEniCS FESpace")
FenExpr = TypeVar("FEniCS expression")
FenFunc = TypeVar("FEniCS function")
# ENGINES
HFEng = TypeVar("High fidelity engine")
-HSEng = TypeVar("Hilbert space engine")
ROMEng = TypeVar("ROM engine")
-LAEng = TypeVar("Linear algebra engine")
# OTHERS
GenExpr = TypeVar("Generic expression")
strLst = TypeVar("str or list of str")
BfSExpr = TypeVar("Boolean function or string")
diff --git a/rrompy/utilities/parameter_sampling/fft_sampler.py b/rrompy/utilities/parameter_sampling/fft_sampler.py
index 9cc17fb..b028726 100644
--- a/rrompy/utilities/parameter_sampling/fft_sampler.py
+++ b/rrompy/utilities/parameter_sampling/fft_sampler.py
@@ -1,56 +1,61 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import Np1D, List
__all__ = ['FFTSampler']
class FFTSampler(GenericSampler):
"""Generator of FFT-type sample points on scaled roots of unity."""
def generatePoints(self, n:Np1D) -> List[Np1D]:
"""Array of sample points and array of weights."""
+ super().generatePoints(n)
d = len(self.lims[0])
try:
len(n)
except:
n = np.array([n])
if len(n) != d:
raise Exception(("Numbers of points must have same dimension as "
"limits."))
for j in range(d):
a, b = self.lims[0][j], self.lims[1][j]
+ if self.scaling is not None:
+ a, b = self.scaling[j](a), self.scaling[j](b)
c, r = (a + b) / 2., np.abs(a - b) / 2.
xj = c + r * np.exp(1.j *
np.linspace(0, 2 * np.pi, n[j] + 1)[:-1, None])
wj = r / n[j] * np.ones(n[j])
+ if self.scalingInv is not None:
+ xj = self.scalingInv[j](xj)
if j == 0:
x = xj
w = wj
xsize = n[0]
else:
x = np.concatenate((np.kron(np.ones(n[j])[:, None], x),
np.kron(xj, np.ones(xsize)[:, None])),
axis = 1)
w = np.multiply(np.kron(np.ones(n[j]), w),
np.kron(wj, np.ones(xsize)))
xsize = xsize * n[j]
return [y.flatten() for y in np.split(x, xsize)], w
diff --git a/rrompy/utilities/parameter_sampling/generic_sampler.py b/rrompy/utilities/parameter_sampling/generic_sampler.py
index 4d65876..f035001 100644
--- a/rrompy/utilities/parameter_sampling/generic_sampler.py
+++ b/rrompy/utilities/parameter_sampling/generic_sampler.py
@@ -1,66 +1,94 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
import numpy as np
-from rrompy.utilities.base.types import Np1D, Tuple, GenExpr
+from rrompy.utilities.base.types import Np1D, Tuple, GenExpr, List
__all__ = ['GenericSampler']
class GenericSampler:
"""ABSTRACT. Generic generator of sample points."""
- def __init__(self, lims:Tuple[Np1D, Np1D]):
+ def __init__(self, lims:Tuple[Np1D, Np1D], scaling : List[callable] = None,
+ scalingInv : List[callable] = None):
self.lims = lims
+ self.scaling = scaling
+ self.scalingInv = scalingInv
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return "{}[{}{}]".format(self.name(),
np.array2string(self.lims[0], separator = '_'),
np.array2string(self.lims[1], separator = '_'))
def __eq__(self, other) -> bool:
return self.__dict__ == other.__dict__
@property
def lims(self):
"""Value of lims."""
return self._lims
@lims.setter
def lims(self, lims):
if len(lims) != 2:
raise Exception("2 limits must be specified.")
try:
lims = lims.tolist()
except:
lims = list(lims)
for j in range(2):
try:
len(lims[j])
except:
lims[j] = np.array([lims[j]])
if len(lims[0]) != len(lims[1]):
raise Exception("The limits must have the same length.")
self._lims = lims
+ @property
+ def scaling(self):
+ """Value of scaling."""
+ return self._scaling
+ @scaling.setter
+ def scaling(self, scaling):
+ if scaling is not None and not isinstance(scaling, list):
+ scaling = [scaling]
+ self._scaling = scaling
+
+ @property
+ def scalingInv(self):
+ """Value of scalingInv."""
+ return self._scalingInv
+ @scalingInv.setter
+ def scalingInv(self, scalingInv):
+ if scalingInv is not None and not isinstance(scalingInv, list):
+ scalingInv = [scalingInv]
+ self._scalingInv = scalingInv
+
@abstractmethod
def generatePoints(self, n:GenExpr) -> Tuple[Np1D, Np1D]:
"""Array of points and array of weights."""
+ assert ((self.scaling is None
+ or len(self.scaling) == len(self.lims[0]))
+ and (self.scalingInv is None
+ or len(self.scalingInv) == len(self.lims[0])))
pass
+
diff --git a/rrompy/utilities/parameter_sampling/quadrature_sampler.py b/rrompy/utilities/parameter_sampling/quadrature_sampler.py
index f69b887..ee92e22 100644
--- a/rrompy/utilities/parameter_sampling/quadrature_sampler.py
+++ b/rrompy/utilities/parameter_sampling/quadrature_sampler.py
@@ -1,79 +1,87 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import Np1D, Tuple, List
__all__ = ['QuadratureSampler']
class QuadratureSampler(GenericSampler):
"""Generator of quadrature sample points."""
allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"]
- def __init__(self, lims:Tuple[Np1D, Np1D], kind : str = "UNIFORM"):
- super().__init__(lims = lims)
+ def __init__(self, lims:Tuple[Np1D, Np1D], kind : str = "UNIFORM",
+ scaling : List[callable] = None,
+ scalingInv : List[callable] = None):
+ super().__init__(lims = lims, scaling = scaling,
+ scalingInv = scalingInv)
self.kind = kind
@property
def kind(self):
"""Value of kind."""
return self._kind
@kind.setter
def kind(self, kind):
if kind.upper() not in self.allowedKinds:
raise Exception("Generator kind not recognized.")
self._kind = kind.upper()
def generatePoints(self, n:Np1D) -> List[Np1D]:
"""Array of quadrature points and array of weights."""
+ super().generatePoints(n)
d = len(self.lims[0])
try:
len(n)
except:
n = np.array([n])
if len(n) != d:
raise Exception(("Numbers of points must have same dimension as"
"limits."))
for j in range(d):
a, b = self.lims[0][j], self.lims[1][j]
+ if self.scaling is not None:
+ a, b = self.scaling[j](a), self.scaling[j](b)
if self.kind == "UNIFORM":
xj = np.linspace(a, b, n[j])[:, None]
wj = np.abs(a - b) / n[j] * np.ones(n[j])
elif self.kind == "CHEBYSHEV":
nodes, weights = np.polynomial.chebyshev.chebgauss(n[j])
xj = (a + b) / 2 + (a - b) / 2 * nodes[:, None]
wj = np.abs(a - b) / np.pi * weights
elif self.kind == "GAUSSLEGENDRE":
nodes, weights = np.polynomial.legendre.leggauss(n[j])
xj = (a + b) / 2 + (a - b) / 2 * nodes[:, None]
wj = np.abs(a - b) * weights
+ if self.scalingInv is not None:
+ xj = self.scalingInv[j](xj)
if j == 0:
x = xj
w = wj
xsize = n[0]
else:
x = np.concatenate((np.kron(np.ones(n[j])[:, None], x),
np.kron(xj, np.ones(xsize)[:, None])),
axis = 1)
w = np.multiply(np.kron(np.ones(n[j]), w),
np.kron(wj, np.ones(xsize)))
xsize = xsize * n[j]
return [y.flatten() for y in np.split(x, xsize)], w
diff --git a/rrompy/utilities/parameter_sampling/random_sampler.py b/rrompy/utilities/parameter_sampling/random_sampler.py
index 4073f4a..13cc5e4 100644
--- a/rrompy/utilities/parameter_sampling/random_sampler.py
+++ b/rrompy/utilities/parameter_sampling/random_sampler.py
@@ -1,57 +1,64 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.base.sobol import sobolGenerate
from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import Np1D, Tuple, List
__all__ = ['RandomSampler']
class RandomSampler(GenericSampler):
"""Generator of quadrature sample points."""
allowedKinds = ["UNIFORM", "SOBOL"]
- def __init__(self, lims:Tuple[Np1D, Np1D], kind : str = "UNIFORM"):
- super().__init__(lims = lims)
+ def __init__(self, lims:Tuple[Np1D, Np1D], kind : str = "UNIFORM",
+ scaling : callable = None, scalingInv : callable = None):
+ super().__init__(lims = lims, scaling = scaling,
+ scalingInv = scalingInv)
self.kind = kind
@property
def kind(self):
"""Value of kind."""
return self._kind
@kind.setter
def kind(self, kind):
if kind.upper() not in self.allowedKinds:
raise Exception("Generator kind not recognized.")
self._kind = kind.upper()
def generatePoints(self, n:int, seed : int = 0) -> List[Np1D]:
"""Array of quadrature points and array of weights."""
+ super().generatePoints(n)
d = len(self.lims[0])
if self.kind == "UNIFORM":
np.random.seed(seed)
x = np.random.uniform(size = (n, d))
else:
x = sobolGenerate(d, n, seed)
for j in range(d):
a, b = self.lims[0][j], self.lims[1][j]
+ if self.scaling is not None:
+ a, b = self.scaling[j](a), self.scaling[j](b)
x[:, j] = a + (b - a) * x[:, j]
+ if self.scalingInv is not None:
+ x[:, j] = self.scalingInv[j](x[:, j])
return [y.flatten() for y in np.split(x, n)], np.ones(n)
diff --git a/rrompy/utilities/parameter_sweeper/parameter_sweeper.py b/rrompy/utilities/parameter_sweeper/parameter_sweeper.py
index 8d2b015..981b5f9 100644
--- a/rrompy/utilities/parameter_sweeper/parameter_sweeper.py
+++ b/rrompy/utilities/parameter_sweeper/parameter_sweeper.py
@@ -1,328 +1,324 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import copy
import itertools
import csv
import warnings
import numpy as np
from matplotlib import pyplot as plt
from rrompy.utilities.base.types import Np1D, N2FSExpr, DictAny, List, ROMEng
from rrompy.utilities.base import purgeList, getNewFilename
__all__ = ['ParameterSweeper']
def C2R2csv(x):
x = np.ravel(x)
y = np.concatenate((np.real(x), np.imag(x)))
z = np.ravel(np.reshape(y, [2, np.size(x)]).T)
return np.array2string(z, separator = '_', suppress_small = False,
max_line_width = np.inf, sign = '+',
formatter = {'all' : lambda x : "{:.15E}".format(x)}
)[1 : -1]
class ParameterSweeper:
"""
ROM approximant parameter sweeper.
Args:
ROMEngine(optional): Generic approximant class. Defaults to None.
mutars(optional): Array of parameter values to sweep. Defaults to empty
array.
params(optional): List of parameter settings (each as a dict) to
explore. Defaults to single empty set.
mostExpensive(optional): String containing label of most expensive
step, to be executed fewer times. Allowed options are 'HF' and
'Approx'. Defaults to 'HF'.
- normType(optional): Target norm identifier. Must be recognizable by
- HSEngine norm command. Defaults to None.
Attributes:
ROMEngine: Generic approximant class.
mutars: Array of parameter values to sweep.
params: List of parameter settings (each as a dict) to explore.
mostExpensive: String containing label of most expensive step, to be
executed fewer times.
"""
allowedOutputsStandard = ["HFNorm", "AppNorm", "ErrNorm"]
allowedOutputs = allowedOutputsStandard + ["HFFunc", "AppFunc", "ErrFunc"]
allowedOutputsFull = allowedOutputs + ["poles"]
def __init__(self, ROMEngine : ROMEng = None, mutars : Np1D = np.array([]),
- params : List[DictAny] = [{}], mostExpensive : str = "HF",
- normType : N2FSExpr = None):
+ params : List[DictAny] = [{}], mostExpensive : str = "HF"):
self.ROMEngine = ROMEngine
self.mutars = mutars
self.params = params
self.mostExpensive = mostExpensive
- self.normType = normType
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
@property
def mostExpensive(self):
"""Value of mostExpensive."""
return self._mostExpensive
@mostExpensive.setter
def mostExpensive(self, mostExpensive:str):
mostExpensive = mostExpensive.upper()
if mostExpensive not in ["HF", "APPROX"]:
warnings.warn(("Value of mostExpensive not recognized. Overriding "
"to 'APPROX'."), stacklevel = 2)
mostExpensive = "APPROX"
self._mostExpensive = mostExpensive
def checkValues(self) -> bool:
"""Check if sweep can be performed."""
if self.ROMEngine is None:
warnings.warn("ROMEngine is missing. Aborting.", stacklevel = 2)
return False
if len(self.mutars) == 0:
warnings.warn("Empty target parameter vector. Aborting.",
stacklevel = 2)
return False
if len(self.params) == 0:
warnings.warn("Empty method parameters vector. Aborting.",
stacklevel = 2)
return False
return True
def sweep(self, filename : str = "out.dat", outputs : List[str] = [],
verbose : int = 1):
if not self.checkValues(): return
try:
if outputs.upper() == "ALL":
outputs = self.allowedOutputsFull
except:
if len(outputs) == 0:
outputs = self.allowedOutputsStandard
outputs = purgeList(outputs, self.allowedOutputsFull,
listname = self.name() + ".outputs",
baselevel = 1)
poles = ("poles" in outputs)
if len(outputs) == 0:
warnings.warn("Empty outputs. Aborting.", stacklevel = 2)
return
outParList = self.ROMEngine.parameterList
Nparams = len(self.params)
if poles: polesCheckList = []
allowedParams = self.ROMEngine.parameterList
dotPos = filename.rfind('.')
if dotPos in [-1, len(filename) - 1]:
filename = getNewFilename(filename[:dotPos])
else:
filename = getNewFilename(filename[:dotPos], filename[dotPos + 1:])
append_write = "w"
initial_row = (outParList + ["muRe", "muIm"]
+ [x for x in self.allowedOutputs if x in outputs]
+ ["type"] + ["poles"] * poles)
with open(filename, append_write, buffering = 1) as fout:
writer = csv.writer(fout, delimiter=",")
writer.writerow(initial_row)
if self.mostExpensive == "HF":
outerSet = self.mutars
innerSet = self.params
elif self.mostExpensive == "APPROX":
outerSet = self.params
innerSet = self.mutars
for outerIdx, outerPar in enumerate(outerSet):
if self.mostExpensive == "HF":
i, mutar = outerIdx, outerPar
elif self.mostExpensive == "APPROX":
j, par = outerIdx, outerPar
self.ROMEngine.approxParameters = {k: par[k] for k in\
par.keys() & allowedParams}
self.ROMEngine.setupApprox()
for innerIdx, innerPar in enumerate(innerSet):
if self.mostExpensive == "APPROX":
i, mutar = innerIdx, innerPar
elif self.mostExpensive == "HF":
j, par = innerIdx, innerPar
self.ROMEngine.approxParameters = {k: par[k] for k in\
par.keys() & allowedParams}
self.ROMEngine.setupApprox()
if verbose >= 1:
print("Set {}/{}\tmu_{} = {:.10f}".format(j+1, Nparams,
i, mutar))
outData = []
if "HFNorm" in outputs:
- val = self.ROMEngine.HFNorm(mutar, self.normType)
+ val = self.ROMEngine.HFNorm(mutar)
if isinstance(val, (list,)): val = val[0]
outData = outData + [val]
if "AppNorm" in outputs:
- val = self.ROMEngine.approxNorm(mutar, self.normType)
+ val = self.ROMEngine.approxNorm(mutar)
if isinstance(val, (list,)): val = val[0]
outData = outData + [val]
if "ErrNorm" in outputs:
- val = self.ROMEngine.approxError(mutar, self.normType)
+ val = self.ROMEngine.approxError(mutar)
if isinstance(val, (list,)): val = val[0]
outData = outData + [val]
if "HFFunc" in outputs:
outData = outData +[self.ROMEngine.HFEngine.functional(
self.ROMEngine.getHF(mutar))]
if "AppFunc" in outputs:
outData = outData +[self.ROMEngine.HFEngine.functional(
self.ROMEngine.getApp(mutar))]
if "ErrFunc" in outputs:
outData = outData +[self.ROMEngine.HFEngine.functional(
self.ROMEngine.getApp(mutar))
- self.ROMEngine.HFEngine.functional(
self.ROMEngine.getHF(mutar))]
writeData = []
for parn in outParList:
writeData = (writeData
+ [self.ROMEngine.approxParameters[parn]])
writeData = (writeData + [mutar.real, mutar.imag]
+ outData + [self.ROMEngine.name()])
if poles:
if j not in polesCheckList:
polesCheckList += [j]
writeData = writeData + [C2R2csv(
self.ROMEngine.getPoles())]
else:
writeData = writeData + [""]
writer.writerow(str(x) for x in writeData)
if verbose >= 1:
if self.mostExpensive == "APPROX":
print("Set {}/{}\tdone".format(j+1, Nparams))
elif self.mostExpensive == "HF":
print("Point mu_{} = {:.10f}\tdone".format(i, mutar))
self.filename = filename
return self.filename
def read(self, filename:str, restrictions : DictAny = {},
outputs : List[str] = []) -> DictAny:
"""
Execute a query on a custom format CSV.
Args:
filename: CSV filename.
restrictions(optional): Parameter configurations to output.
Defaults to empty dictionary, i.e. output all.
outputs(optional): Values to output. Defaults to empty list, i.e.
no output.
Returns:
Dictionary of desired results, with a key for each entry of
outputs, and a numpy 1D array as corresponding value.
"""
with open(filename, 'r') as f:
reader = csv.reader(f, delimiter=',')
header = next(reader)
restrIndices, outputIndices, outputData = {}, {}, {}
for key in restrictions.keys():
try:
restrIndices[key] = header.index(key)
if not isinstance(restrictions[key], list):
restrictions[key] = [restrictions[key]]
restrictions[key] = copy(restrictions[key])
except:
warnings.warn("Ignoring key {} from restrictions"\
.format(key), stacklevel = 2)
for key in outputs:
try:
outputIndices[key] = header.index(key)
outputData[key] = np.array([])
except:
warnings.warn("Ignoring key {} from outputs".format(key),
stacklevel = 2)
for row in reader:
restrTrue = True
for key in restrictions.keys():
if row[restrIndices[key]] == restrictions[key]:
continue
try:
if np.any(np.isclose(float(row[restrIndices[key]]),
[float(x) for x in restrictions[key]])):
continue
except: pass
restrTrue = False
if restrTrue:
for key in outputIndices.keys():
try:
val = row[outputIndices[key]]
val = float(val)
finally:
outputData[key] = np.append(outputData[key], val)
return outputData
def plot(self, filename:str, xs:List[str], ys:List[str], zs:List[str],
onePlot : bool = False):
"""
Perform plots from data in filename.
Args:
filename: CSV filename.
xs: Values to put on x axes.
ys: Values to put on y axes.
zs: Meta-values for constraints.
onePlot: Whether to create a single figure per x. Defaults to
False.
"""
zsVals = self.read(filename, outputs = zs)
zs = list(zsVals.keys())
zss = None
for key in zs:
vals = np.unique(zsVals[key])
if zss is None:
zss = copy(vals)
else:
zss = list(itertools.product(zss, vals))
lzs = len(zs)
for z in zss:
if lzs <= 1:
constr = {zs[0] : z}
else:
constr = {zs[j] : z[j] for j in range(len(zs))}
data = self.read(filename, restrictions = constr, outputs = xs+ys)
if onePlot:
for x in xs:
xVals = data[x]
plt.figure()
for y in ys:
yVals = data[y]
label = '{} vs {} for {}'.format(x, y, constr)
plt.semilogy(xVals, yVals, label = label)
plt.legend()
plt.grid()
plt.show()
plt.close()
else:
for x, y in itertools.product(xs, ys):
xVals, yVals = data[x], data[y]
label = '{} vs {} for {}'.format(x, y, constr)
plt.figure()
plt.semilogy(xVals, yVals, label = label)
plt.legend()
plt.grid()
plt.show()
plt.close()