diff --git a/examples/airfoil/airfoil_engine.py b/examples/airfoil/airfoil_engine.py
index 17bcce9..f63b57d 100644
--- a/examples/airfoil/airfoil_engine.py
+++ b/examples/airfoil/airfoil_engine.py
@@ -1,50 +1,50 @@
import numpy as np
import fenics as fen
import ufl
from rrompy.hfengines.linear_problem import \
HelmholtzBoxScatteringProblemEngine as HSP
-from rrompy.utilities.fenics import fenONE
+from rrompy.solver.fenics import fenONE
PI = np.pi
class AirfoilScatteringEngine(HSP):
def __init__(self, kappa:float, theta:float,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(R = 5, kappa = kappa, theta = theta, n = 1,
degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
mesh = fen.Mesh('../data/mesh/airfoil2412_1.xml')
self.V = fen.FunctionSpace(mesh, "P", 1)
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
u0R = - fen.cos(kappa * (c * x + s * y))
u0I = - fen.sin(kappa * (c * x + s * y))
self.DirichletDatum = [u0R, u0I]
mu = 1.1
epsilon = .1
checkReal = x**2-x+y**2
rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25
phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2
phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2
kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/
((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.)
)**.5 - mu
kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/
((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.)
)**.5 - mu
Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1
Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1
cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE)
c_F = fen.Constant(.1)
cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE)
c_F = fen.Constant(.1)
cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F)
cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F)
a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF)
self.diffusivity = a
diff --git a/examples/airfoil/greedy.py b/examples/airfoil/greedy.py
index df59726..373d19c 100644
--- a/examples/airfoil/greedy.py
+++ b/examples/airfoil/greedy.py
@@ -1,106 +1,106 @@
import numpy as np
from airfoil_engine import AirfoilScatteringEngine
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
-from rrompy.utilities.fenics import L2NormMatrix
+from rrompy.solver.fenics import L2NormMatrix
verb = 2
timed = False
algo = "Pade"
#algo = "RB"
polyBasis = "LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
homog = True
homog = False
k0s = np.linspace(5, 20, 100)
k0 = np.mean(k0s)
kl, kr = min(k0s), max(k0s)
params = {'muBounds':[kl, kr], 'nTestPoints':500, 'Delta':0,
'greedyTol':1e-2, 'S':2, 'basis':polyBasis}
#########
kappa = 10
theta = np.pi * - 45 / 180.
solver = AirfoilScatteringEngine(kappa, theta, verbosity = verb,
degree_threshold = 8)
#########
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb, homogeneized = homog)
else:
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb,
homogeneized = homog)
approx.initEstimatorNormEngine(L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
kl, kr = np.real(kl), np.real(kr)
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(
approx.getRes(k0s[j], homogeneized=homog))
/ approx.estimatorNormEngine.norm(
approx.getRHS(k0s[j], homogeneized=homog)))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.plot(k0s, norm)
plt.plot(k0s, normApp, '--')
plt.plot(np.real(approx.mus),
1.05*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/examples/diapason/greedy.py b/examples/diapason/greedy.py
index 473d5aa..cd24a8c 100644
--- a/examples/diapason/greedy.py
+++ b/examples/diapason/greedy.py
@@ -1,179 +1,179 @@
import numpy as np
import fenics as fen
import ufl
from rrompy.hfengines.vector_linear_problem import \
LinearElasticityHelmholtzProblemEngine as LEHPE
from rrompy.hfengines.vector_linear_problem import \
LinearElasticityHelmholtzProblemEngineDamped as LEHPED
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
-from rrompy.utilities.fenics import L2NormMatrix
+from rrompy.solver.fenics import L2NormMatrix
verb = 2
timed = False
algo = "Pade"
#algo = "RB"
polyBasis = "LEGENDRE"
polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
if timed: verb = 0
dampingEta = 0 * 1e4 / 2. / np.pi
k0s = np.linspace(2.5e2, 7.5e3, 100)
k0s = np.linspace(2.5e3, 1.5e4, 100)
k0s = np.linspace(5.0e4, 1.0e5, 100)
k0s = np.linspace(2.0e5, 2.5e5, 100)
k0 = np.mean(np.power(k0s, 2.)) ** .5 # [Hz]
kl, kr = min(k0s), max(k0s)
params = {'muBounds':[kl, kr], 'nTestPoints':500, 'Delta':0,
'greedyTol':1e-2, 'S':2, 'polybasis':polyBasis,
'robustTol':2e-16, 'interpRcond':None, 'errorEstimatorKind':'EXACT'}
theta = 20. * np.pi / 180.
phi = 10. * np.pi / 180.
mesh = fen.Mesh("../data/mesh/diapason_1.xml")
subdomains = fen.MeshFunction("size_t", mesh,
"../data/mesh/diapason_1_physical_region.xml")
meshBall = fen.SubMesh(mesh, subdomains, 2)
meshFork = fen.SubMesh(mesh, subdomains, 1)
Hball = np.max(meshBall.coordinates()[:, 1]) #.00257
Ltot = np.max(mesh.coordinates()[:, 1]) #.1022
Lhandle = np.max(meshFork.coordinates()[:, 1]) #.026
Lrod = Ltot - Lhandle #.0762
L2var = (Lrod / 4.) ** 2.
Ehandle_ratio = 3.
rhohandle_ratio = 1.5
c = 3.e2
rho = 8e3 * (2. * np.pi) ** 2.
E = 1.93e11
nu = .3
T = 1e6
lambda_ = E * nu / (1. + nu) / (1. - 2. * nu)
mu_ = E / (1. + nu)
kWave = (np.cos(theta) * np.cos(phi), np.sin(phi), np.sin(theta) * np.cos(phi))
x, y, z = fen.SpatialCoordinate(mesh)[:]
yCorr = y - Ltot
compPlane = kWave[0] * x + kWave[1] * y + kWave[2] * z
xPlane, yPlane, zPlane = (xx - compPlane * xx for xx in (x, y, z))
xOrtho, yOrtho, zOrtho = (compPlane * xx for xx in (x, y, z))
forcingBase = (T / (2. * np.pi * L2var)**.5
* fen.exp(- (xPlane**2. + yPlane**2. + zPlane**2.) / (2.*L2var)))
forcingWeight = np.real(k0) / c * (xOrtho + yOrtho + zOrtho)
neumannDatum = [ufl.as_vector(
tuple(forcingBase * fen.cos(forcingWeight) * kWavedir for kWavedir in kWave)),
ufl.as_vector(
tuple(forcingBase * fen.sin(forcingWeight) * kWavedir for kWavedir in kWave))]
lambda_eff = ufl.conditional(ufl.ge(y, Lhandle), fen.Constant(lambda_),
fen.Constant(Ehandle_ratio * lambda_))
mu_eff = ufl.conditional(ufl.ge(y, Lhandle), fen.Constant(mu_),
fen.Constant(Ehandle_ratio * mu_))
rho_eff = ufl.conditional(ufl.ge(y, Lhandle), fen.Constant(rho),
fen.Constant(rhohandle_ratio * rho))
###
if dampingEta > 0:
solver = LEHPED(degree_threshold = 8, verbosity = 0)
solver.eta = dampingEta
else:
solver = LEHPE(degree_threshold = 8, verbosity = 0)
solver.omega = np.real(k0)
solver.lambda_ = lambda_eff
solver.mu_ = mu_eff
solver.rho_ = rho_eff
solver.V = fen.VectorFunctionSpace(mesh, "P", 1)
solver.DirichletBoundary = lambda x, on_b: on_b and x[1] < Hball
solver.NeumannBoundary = "REST"
solver.forcingTerm = neumannDatum
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
else:
params.pop("Delta")
params.pop("polybasis")
params.pop("robustTol")
params.pop("interpRcond")
params.pop("errorEstimatorKind")
approx = RB(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.initEstimatorNormEngine(L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
polesApp = approx.getPoles()
print("Poles:\n", polesApp)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
kl, kr = np.real(kl), np.real(kr)
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j]))
/ approx.estimatorNormEngine.norm(approx.getRHS(k0s[j])))
err[j] = approx.normErr(k0s[j]) / norm[j]
resApp = approx.errorEstimator(k0s)
res[res < 1e-5 * approx.greedyTol] = np.nan
resApp[resApp < 1e-5 * approx.greedyTol] = np.nan
err[err < 1e-8 * approx.greedyTol] = np.nan
plt.figure()
plt.semilogy(k0s, norm)
plt.semilogy(k0s, normApp, '--')
plt.semilogy(np.real(approx.mus),
1.05*np.max(norm)*np.ones_like(approx.mus, dtype = float),
'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus),
approx.greedyTol*np.ones_like(approx.mus, dtype = float),
'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesAppEff = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesAppEff), np.imag(polesAppEff), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/examples/from_papers/greedy_internalBox.py b/examples/from_papers/greedy_internalBox.py
index 9532548..72c798a 100644
--- a/examples/from_papers/greedy_internalBox.py
+++ b/examples/from_papers/greedy_internalBox.py
@@ -1,112 +1,112 @@
import numpy as np
import fenics as fen
from rrompy.hfengines.linear_problem import HelmholtzProblemEngine as HPE
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
-from rrompy.utilities.fenics import L2NormMatrix
+from rrompy.solver.fenics import L2NormMatrix
dim = 3
verb = 2
timed = False
algo = "Pade"
#algo = "RB"
polyBasis = "LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
k0s = np.power(np.linspace(500 ** 2., 2250 ** 2., 200), .5)
k0 = np.mean(np.power(k0s, 2.)) ** .5
kl, kr = min(k0s), max(k0s)
params = {'muBounds':[kl, kr], 'nTestPoints':500, 'Delta':0,
'greedyTol':1e-2, 'S':2, 'basis':polyBasis}
if dim == 2:
mesh = fen.RectangleMesh(fen.Point(0., 0.), fen.Point(.1, .15), 10, 15)
x, y = fen.SpatialCoordinate(mesh)[:]
f = fen.exp(- 1e2 * (x + y))
else:#if dim == 3:
mesh = fen.BoxMesh(fen.Point(0., 0., 0.), fen.Point(.1, .15, .25), 4, 6,10)
x, y, z = fen.SpatialCoordinate(mesh)[:]
f = fen.exp(- 1e2 * (x + y + z))
solver = HPE(verbosity = verb)
solver.omega = np.real(k0)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.refractionIndex = fen.Constant(1. / 54.6)
solver.forcingTerm = f
solver.NeumannBoundary = "ALL"
#########
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
else:
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approx.initEstimatorNormEngine(L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
kl, kr = np.real(kl), np.real(kr)
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j]))
/ approx.estimatorNormEngine.norm(approx.getRHS(k0s[j])))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.plot(k0s, norm)
plt.plot(k0s, normApp, '--')
plt.plot(np.real(approx.mus),
1.05*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/examples/from_papers/greedy_scatteringAirfoil.py b/examples/from_papers/greedy_scatteringAirfoil.py
index c9ea8b7..e261c9a 100644
--- a/examples/from_papers/greedy_scatteringAirfoil.py
+++ b/examples/from_papers/greedy_scatteringAirfoil.py
@@ -1,150 +1,149 @@
import numpy as np
import fenics as fen
import ufl
from rrompy.hfengines.linear_problem import \
HelmholtzBoxScatteringProblemEngine as HSP
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
-from rrompy.utilities.fenics import fenONE
-from rrompy.utilities.fenics import L2NormMatrix
+from rrompy.solver.fenics import fenONE, L2NormMatrix
verb = 2
timed = False
algo = "Pade"
#algo = "RB"
polyBasis = "LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
homog = True
homog = False
k0s = np.linspace(5, 20, 25)
k0 = np.mean(k0s)
kl, kr = min(k0s), max(k0s)
params = {'muBounds':[kl, kr], 'nTestPoints':500, 'Delta':0,
'greedyTol':1e-2, 'S':2, 'basis':polyBasis}
#########
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
mesh = fen.Mesh('../data/mesh/airfoil.xml')
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
u0R = - fen.cos(kappa * (c * x + s * y))
u0I = - fen.sin(kappa * (c * x + s * y))
checkReal = x**2-x+y**2
rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25
phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2
phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2
kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/
((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.)
)**.5 - mu
kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/
((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.)
)**.5 - mu
Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1
Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1
cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE)
c_F = fen.Constant(.1)
cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE)
c_F = fen.Constant(.1)
cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F)
cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F)
a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF)
solver = HSP(R, np.real(k0), theta, n = 1, verbosity = verb,
degree_threshold = 8)
solver.omega = np.real(k0)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.diffusivity = a
solver.DirichletBoundary = Dboundary
solver.RobinBoundary = "REST"
solver.DirichletDatum = [u0R, u0I]
#########
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb, homogeneized = homog)
else:
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb,
homogeneized = homog)
approx.initEstimatorNormEngine(L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
kl, kr = np.real(kl), np.real(kr)
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(
approx.getRes(k0s[j], homogeneized=homog))
/ approx.estimatorNormEngine.norm(
approx.getRHS(k0s[j], homogeneized=homog)))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.plot(k0s, norm)
plt.plot(k0s, normApp, '--')
plt.plot(np.real(approx.mus),
1.05*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/examples/from_papers/pod_scatteringAirfoil.py b/examples/from_papers/pod_scatteringAirfoil.py
index 2f11c01..29d9a57 100644
--- a/examples/from_papers/pod_scatteringAirfoil.py
+++ b/examples/from_papers/pod_scatteringAirfoil.py
@@ -1,144 +1,143 @@
-from copy import deepcopy as copy
import numpy as np
import fenics as fen
import ufl
from rrompy.hfengines.linear_problem import \
HelmholtzBoxScatteringProblemEngine as HSP
from rrompy.reduction_methods.centered import RationalPade as PC
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RationalInterpolant as PD
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
-from rrompy.utilities.fenics import fenONE
+from rrompy.solver.fenics import fenONE
from operator import itemgetter
def subdict(d, ks):
return dict(zip(ks, itemgetter(*ks)(d)))
verb = 0
####################
homog = True
#homog = False
####################
test = "solve"
test = "Centered"
test = "Distributed"
plotSamples = True
k0 = 10
kLeft, kRight = 8 + 0.j, 12 + 0.j
ktar = 11
ktars = np.linspace(8, 12, 21) + 0.j
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
mesh = fen.Mesh('../data/mesh/airfoil.xml')
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
u0R = - fen.cos(kappa * (c * x + s * y))
u0I = - fen.sin(kappa * (c * x + s * y))
checkReal = x**2-x+y**2
rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25
phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2
phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2
kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/
((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.)
)**.5 - mu
kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/
((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.)
)**.5 - mu
Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1
Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1
cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE)
c_F = fen.Constant(.1)
cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE)
c_F = fen.Constant(.1)
cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F)
cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F)
a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF)
###
solver = HSP(R, np.abs(k0), theta, n = 1, verbosity = verb,
degree_threshold = 8)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.diffusivity = a
solver.DirichletBoundary = Dboundary
solver.RobinBoundary = "REST"
solver.DirichletDatum = [u0R, u0I]
###
if test == "solve":
uinc = solver.liftDirichletData(k0)
if homog:
uhtot = solver.solve(k0, homogeneized = homog)
uh = uhtot + uinc
else:
uh = solver.solve(k0, homogeneized = homog)
uhtot = uh - uinc
print(solver.norm(uh))
print(solver.norm(uhtot))
solver.plot(fen.project(a, solver.V).vector(), what = 'Real',
name = 'a')
solver.plot(uinc, what = 'Real', name = 'u_inc')
solver.plot(uh, what = 'ABS')
solver.plot(uhtot, what = 'ABS', name = 'u + u_inc')
elif test in ["Centered", "Distributed"]:
if test == "Centered":
params = {'N':8, 'M':8, 'R':8, 'E':8, 'sampleType':'Arnoldi',
'POD':True}
parPade = subdict(params, ['N', 'M', 'E', 'sampleType', 'POD'])
parRB = subdict(params, ['R', 'E', 'sampleType', 'POD'])
approxPade = PC(solver, mu0 = k0, approxParameters = parPade,
verbosity = verb, homogeneized = homog)
approxRB = RBC(solver, mu0 = k0, approxParameters = parRB,
verbosity = verb, homogeneized = homog)
else:
params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True, 'basis':"CHEBYSHEV",
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
parPade = subdict(params, ['N', 'M', 'S', 'POD', 'basis', 'sampler'])
parRB = subdict(params, ['R', 'S', 'POD', 'sampler'])
approxPade = PD(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parPade, verbosity = verb,
homogeneized = homog)
approxRB = RBD(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parRB, verbosity = verb,
homogeneized = homog)
approxPade.setupApprox()
approxRB.setupApprox()
if plotSamples:
approxPade.plotSamples()
approxPade.plotHF(ktar, name = 'u_HF')
approxPade.plotApprox(ktar, name = 'u_Pade''')
approxPade.plotErr(ktar, name = 'err_Pade''')
approxPade.plotRes(ktar, name = 'res_Pade''')
approxRB.plotApprox(ktar, name = 'u_RB')
approxRB.plotErr(ktar, name = 'err_RB')
approxRB.plotRes(ktar, name = 'res_RB')
HFNorm, RHSNorm = approxPade.normHF(ktar), approxPade.normRHS(ktar)
PadeRes, PadeErr = approxPade.normRes(ktar), approxPade.normErr(ktar)
RBRes, RBErr = approxRB.normRes(ktar), approxRB.normErr(ktar)
print('HFNorm:\t{}\nRHSNorm:\t{}'.format(HFNorm, RHSNorm))
print('PadeRes:\t{}\nPadeErr:\t{}'.format(PadeRes, PadeErr))
print('RBRes:\t{}\nRBErr:\t{}'.format(RBRes, RBErr))
print('\nPoles Pade'':')
print(approxPade.getPoles())
diff --git a/examples/greedy/matrix_greedy.py b/examples/greedy/matrix_greedy.py
index 7754260..7419cd8 100644
--- a/examples/greedy/matrix_greedy.py
+++ b/examples/greedy/matrix_greedy.py
@@ -1,111 +1,112 @@
import numpy as np
import scipy.sparse as sp
from matplotlib import pyplot as plt
from rrompy.hfengines.base import MatrixEngineBase as MEB
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
test = 1
timed = False
method = "Pade"
#method = "RB"
verb = 2
errorEstimatorKind = "BARE"
#errorEstimatorKind = "BASIC"
#errorEstimatorKind = "EXACT"
N = 100
solver = MEB(verbosity = verb)
solver.nAs = 2
if test == 1:
solver.As = [sp.spdiags([np.arange(1, 1 + N)], [0], N, N),
- sp.eye(N)]
elif test == 2:
solver.setSolver("SOLVE")
fftB = np.fft.fft(np.eye(N)) * N**-.5
solver.As = [fftB.dot(np.multiply(np.arange(1, 1 + N), fftB.conj()).T),
- np.eye(N)]
np.random.seed(420)
solver.nbs = 1
solver.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
mu0 = 10.25
murange = [1.25, 19.25]
mutars = np.linspace(murange[0], murange[1], 500)
if method == "Pade":
params = {'muBounds':murange, 'nTestPoints':200, 'Delta':0, 'S':5,
- 'greedyTol':1e-2, 'polybasis':"CHEBYSHEV",
- 'errorEstimatorKind':errorEstimatorKind}
+ 'greedyTol':1e-2, 'polybasis':"CHEBYSHEV", "maxIter":17,
+ 'coarseTolerance': .01, 'errorEstimatorKind':errorEstimatorKind}
approx = Pade(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
elif method == "RB":
- params = {'muBounds':murange, 'nTestPoints':500, 'greedyTol':1e-2, 'S':5}
+ params = {'muBounds':murange, 'nTestPoints':500, 'greedyTol':1e-2, 'S':5,
+ 'coarseTolerance': .01}
approx = RB(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
normApp = np.zeros(len(mutars))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(mutars)):
normApp[j] = approx.normApprox(mutars[j])
norm[j] = approx.normHF(mutars[j])
res[j] = (approx.estimatorNormEngine.norm(approx.getRes(mutars[j]))
/ approx.estimatorNormEngine.norm(approx.getRHS(mutars[j])))
err[j] = approx.normErr(mutars[j]) / approx.normHF(mutars[j])
resApp = approx.errorEstimator(mutars)
plt.figure()
plt.semilogy(mutars, norm)
plt.semilogy(mutars, normApp, '--')
plt.semilogy(np.real(approx.mus),
1.25*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim(murange)
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(mutars, res)
plt.semilogy(mutars, resApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim(murange)
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(mutars, err)
plt.xlim(murange)
plt.grid()
plt.show()
plt.close()
polesTrue = np.arange(1, 1 + N)
polesTrue = polesTrue[polesTrue >= murange[0]]
polesTrue = polesTrue[polesTrue <= murange[1]]
polesApp = approx.getPoles()
mask = (np.real(polesApp) < murange[0]) | (np.real(polesApp) > murange[1])
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.plot(polesTrue, np.zeros_like(polesTrue), 'm.')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/examples/greedy/parametricDomain.py b/examples/greedy/parametricDomain.py
index 8b5d640..d1d5910 100644
--- a/examples/greedy/parametricDomain.py
+++ b/examples/greedy/parametricDomain.py
@@ -1,102 +1,102 @@
import numpy as np
from rrompy.hfengines.linear_problem import \
HelmholtzSquareBubbleDomainProblemEngine as HSBDPE
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
-from rrompy.utilities.fenics import L2NormMatrix
+from rrompy.solver.fenics import L2NormMatrix
verb = 2
timed = False
algo = "Pade"
#algo = "RB"
polyBasis = "LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
errorEstimatorKind = "BARE"
#errorEstimatorKind = "BASIC"
#errorEstimatorKind = "EXACT"
k0s = np.power(np.linspace(9, 19, 100), .5)
k0 = np.mean(np.power(k0s, 2.)) ** .5
kl, kr = min(k0s), max(k0s)
params = {'muBounds':[kl, kr], 'nTestPoints':500, 'Delta':0,
'greedyTol':1e-2, 'S':2, 'polybasis':polyBasis,
'errorEstimatorKind':errorEstimatorKind}
if timed:
verb = 0
solver = HSBDPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20, mu0 = k0,
degree_threshold = 15, verbosity = verb)
solver.omega = np.real(k0)
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
else:
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approx.initEstimatorNormEngine(np.abs(k0) ** 2. * L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j]))
/ approx.estimatorNormEngine.norm(approx.getRHS(k0s[j])))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.semilogy(k0s, norm)
plt.semilogy(k0s, normApp, '--')
plt.semilogy(np.real(approx.mus),
1.25*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/examples/greedy/squareBubbleHomog.py b/examples/greedy/squareBubbleHomog.py
index 7a463b5..800fbf0 100644
--- a/examples/greedy/squareBubbleHomog.py
+++ b/examples/greedy/squareBubbleHomog.py
@@ -1,114 +1,114 @@
import numpy as np
from rrompy.hfengines.linear_problem import \
HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
from rrompy.utilities.base import squareResonances
-from rrompy.utilities.fenics import L2NormMatrix
+from rrompy.solver.fenics import L2NormMatrix
verb = 2
timed = False
algo = "Pade"
#algo = "RB"
polyBasis = "LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
errorEstimatorKind = "BARE"
#errorEstimatorKind = "BASIC"
#errorEstimatorKind = "EXACT"
k0s = np.power(np.linspace(95, 149, 250), .5)
#k0s = np.power(np.linspace(95, 129, 100), .5)
#k0s = np.power(np.linspace(95, 109, 100), .5)
k0 = np.mean(np.power(k0s, 2.)) ** .5
kl, kr = min(k0s), max(k0s)
polesexact = np.unique(np.power(squareResonances(kl**2., kr**2., False), .5))
params = {'muBounds':[kl, kr], 'nTestPoints':500, 'Delta':0,
'greedyTol':1e-2, 'S':10, 'polybasis':polyBasis,
'errorEstimatorKind':errorEstimatorKind, 'interactive':False}
if timed:
verb = 0
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20,
verbosity = verb)
solver.omega = np.real(k0)
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
else:
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approx.initEstimatorNormEngine(L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.trainedModel.verbosity = 0
approx.verbosity = 0
from matplotlib import pyplot as plt
normApp = np.zeros_like(k0s)
norm = np.zeros_like(k0s)
res = np.zeros_like(k0s)
err = np.zeros_like(k0s)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j]))
/ approx.estimatorNormEngine.norm(approx.getRHS(k0s[j])))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.semilogy(k0s, norm)
plt.semilogy(k0s, normApp, '--')
plt.semilogy(polesexact,
2.*np.max(norm)*np.ones_like(polesexact, dtype = float), 'm.')
plt.semilogy(np.real(approx.mus),
4.*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(polesexact,
2.*np.max(resApp)*np.ones_like(polesexact, dtype = float), 'm.')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.semilogy(polesexact,
2.*np.max(err)*np.ones_like(polesexact, dtype = float), 'm.')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.plot(np.real(polesexact), np.imag(polesexact), 'm.')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/examples/greedy/squareScatteringHomog.py b/examples/greedy/squareScatteringHomog.py
index aebfed4..f0f2ca7 100644
--- a/examples/greedy/squareScatteringHomog.py
+++ b/examples/greedy/squareScatteringHomog.py
@@ -1,104 +1,104 @@
import numpy as np
from rrompy.hfengines.linear_problem import \
HelmholtzCavityScatteringProblemEngine as HCSPE
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
-from rrompy.utilities.fenics import L2NormMatrix
+from rrompy.solver.fenics import L2NormMatrix
verb = 2
timed = False
algo = "Pade"
#algo = "RB"
polyBasis = "LEGENDRE"
polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
errorEstimatorKind = "BARE"
errorEstimatorKind = "BASIC"
#errorEstimatorKind = "EXACT"
k0s = np.linspace(10, 15, 100)
k0 = np.mean(k0s)
kl, kr = min(k0s), max(k0s)
params = {'muBounds':[kl, kr], 'nTestPoints':500, 'Delta':0,
'greedyTol':5e-3, 'S':2, 'polybasis':polyBasis,
'errorEstimatorKind':errorEstimatorKind}
if timed:
verb = 0
solver = HCSPE(kappa = 5, n = 20, verbosity = verb)
solver.omega = np.real(k0)
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
else:
params.pop('Delta')
params.pop('polybasis')
params.pop('errorEstimatorKind')
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
approx.initEstimatorNormEngine(L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.trainedModel.verbosity = 0
approx.verbosity = 0
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j]))
/ approx.estimatorNormEngine.norm(approx.getRHS(k0s[j])))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
plt.figure()
plt.plot(k0s, norm)
plt.plot(k0s, normApp, '--')
plt.plot(np.real(approx.mus),
1.25*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/examples/greedy/squareTransmissionNonHomog.py b/examples/greedy/squareTransmissionNonHomog.py
index 5abfc82..a6f9804 100644
--- a/examples/greedy/squareTransmissionNonHomog.py
+++ b/examples/greedy/squareTransmissionNonHomog.py
@@ -1,113 +1,113 @@
import numpy as np
from rrompy.hfengines.linear_problem import \
HelmholtzSquareTransmissionProblemEngine as HSTPE
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
-from rrompy.utilities.fenics import L2NormMatrix
+from rrompy.solver.fenics import L2NormMatrix
timed = False
verb = 2
algo = "Pade"
#algo = "RB"
polyBasis = "LEGENDRE"
#polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
homog = True
#homog = False
errorEstimatorKind = "BARE"
errorEstimatorKind = "BASIC"
#errorEstimatorKind = "EXACT"
k0s = np.power(np.linspace(4, 15, 100), .5)
k0 = np.mean(np.power(k0s, 2.)) ** .5
kl, kr = min(k0s), max(k0s)
rescaling = lambda x: np.power(x, 2.)
rescalingInv = lambda x: np.power(x, .5)
params = {'muBounds':[kl, kr], 'nTestPoints':500, 'Delta':0,
'greedyTol':1e-2, 'S':5, 'polybasis':polyBasis,
'errorEstimatorKind':errorEstimatorKind}
solver = HSTPE(nT = 1, nB = 2, theta = np.pi * 20 / 180, kappa = 4.,
n = 20, verbosity = verb)
solver.omega = np.real(k0)
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb, homogeneized = homog)
else:
approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb,
homogeneized = homog)
approx.initEstimatorNormEngine(L2NormMatrix(solver.V))
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.verbosity = 0
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
res[j] = (approx.estimatorNormEngine.norm(
approx.getRes(k0s[j], homogeneized=homog))
/ approx.estimatorNormEngine.norm(
approx.getRHS(k0s[j], homogeneized=homog)))
err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
resApp = approx.errorEstimator(k0s)
polesApp = approx.getPoles()
polesApp = polesApp[np.abs(np.imag(polesApp)) < 1e-3]
plt.figure()
plt.semilogy(k0s, norm)
plt.semilogy(k0s, normApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.semilogy(np.real(polesApp),
2.*np.max(norm)*np.ones_like(polesApp, dtype = float), 'k.')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.semilogy(k0s, resApp, '--')
plt.semilogy(np.real(approx.mus),
4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
plt.semilogy(np.real(polesApp),
2.*np.max(resApp)*np.ones_like(polesApp, dtype = float), 'k.')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
plt.semilogy(np.real(polesApp),
2.*np.max(err)*np.ones_like(polesApp, dtype = float), 'k.')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/rrompy/hfengines/base/boundary_conditions.py b/rrompy/hfengines/base/boundary_conditions.py
index 7ca334f..1c4e302 100644
--- a/rrompy/hfengines/base/boundary_conditions.py
+++ b/rrompy/hfengines/base/boundary_conditions.py
@@ -1,126 +1,126 @@
# 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 fenics import SubDomain, AutoSubDomain
from rrompy.utilities.base.types import GenExpr
-from rrompy.utilities.fenics import bdrFalse
+from rrompy.solver.fenics import bdrFalse
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['BoundaryConditions']
class BoundaryConditions:
"""
Boundary conditions manager.
Attributes:
DirichletBoundary: Callable returning True when on Dirichlet boundary.
NeumannBoundary: Callable returning True when on Neumann boundary.
RobinBoundary: Callable returning True when on Robin boundary.
"""
allowedKinds = ["Dirichlet", "Neumann", "Robin"]
def __init__(self, kind : str = None):
if kind is None: return
kind = kind[0].upper() + kind[1:].lower()
if kind in self.allowedKinds:
getattr(self.__class__, kind + "Boundary", None).fset(self, "ALL")
else:
raise RROMPyException("BC kind not recognized.")
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def _generalManagement(self, kind:str, value:GenExpr):
if isinstance(value, (str,)):
value = value.upper()
if value.upper() == "ALL":
self._complementaryManagementAll(kind)
elif value.upper() == "REST":
self._complementaryManagementRest(kind)
else:
raise RROMPyException("Wildcard not recognized.")
elif callable(value):
self._standardManagementCallable(kind, value)
elif isinstance(value, (SubDomain,)):
self._standardManagement(kind, value)
else:
raise RROMPyException(kind + "Boundary type not recognized.")
def _complementaryManagementAll(self, kind:str):
if kind not in self.allowedKinds:
raise RROMPyException("BC kind not recognized.")
for k in self.allowedKinds:
if k != kind:
self._standardManagementCallable(k, bdrFalse)
self._complementaryManagementRest(kind)
def _complementaryManagementRest(self, kind:str):
if kind not in self.allowedKinds:
raise RROMPyException("BC kind not recognized.")
otherBCs = []
for k in self.allowedKinds:
if k != kind:
if hasattr(self, "_" + k + "Rest"):
self._standardManagement(k, bdrFalse)
otherBCs += [getattr(self, k + "Boundary")]
def restCall(x, on_boundary):
return (on_boundary
and not any([bc.inside(x, on_boundary) for bc in otherBCs]))
self._standardManagementCallable(kind, restCall)
super().__setattr__("_" + kind + "Rest", 1)
def _standardManagementCallable(self, kind:str, bc:callable):
bcSD = AutoSubDomain(bc)
self._standardManagement(kind, bcSD)
def _standardManagement(self, kind:str, bc:SubDomain):
super().__setattr__("_" + kind + "Boundary", bc)
if hasattr(self, "_" + kind + "Rest"):
super().__delattr__("_" + kind + "Rest")
@property
def DirichletBoundary(self):
"""Function handle to DirichletBoundary."""
return self._DirichletBoundary
@DirichletBoundary.setter
def DirichletBoundary(self, DirichletBoundary):
self._generalManagement("Dirichlet", DirichletBoundary)
@property
def NeumannBoundary(self):
"""Function handle to NeumannBoundary."""
return self._NeumannBoundary
@NeumannBoundary.setter
def NeumannBoundary(self, NeumannBoundary):
self._generalManagement("Neumann", NeumannBoundary)
@property
def RobinBoundary(self):
"""Function handle to RobinBoundary."""
return self._RobinBoundary
@RobinBoundary.setter
def RobinBoundary(self, RobinBoundary):
self._generalManagement("Robin", RobinBoundary)
diff --git a/rrompy/hfengines/base/problem_engine_base.py b/rrompy/hfengines/base/problem_engine_base.py
index 8658551..e5acc29 100644
--- a/rrompy/hfengines/base/problem_engine_base.py
+++ b/rrompy/hfengines/base/problem_engine_base.py
@@ -1,361 +1,361 @@
# 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
from os import path, mkdir
import fenics as fen
import numpy as np
from matplotlib import pyplot as plt
from rrompy.utilities.base.types import (Np1D, ScOp, strLst, FenFunc, Tuple,
List)
from rrompy.utilities.base import purgeList, getNewFilename, verbosityDepth
-from rrompy.utilities.fenics import L2NormMatrix
+from rrompy.solver.fenics import L2NormMatrix
from .boundary_conditions import BoundaryConditions
from .matrix_engine_base import MatrixEngineBase
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['ProblemEngineBase']
class ProblemEngineBase(MatrixEngineBase):
"""
Generic solver for parametric problems.
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
"""
def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(verbosity = verbosity, timestamp = timestamp)
self.BCManager = BoundaryConditions("Dirichlet")
self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
self.bsmu = np.nan
self.liftDirichletDatamu = np.nan
self.mu0BC = np.nan
self.degree_threshold = degree_threshold
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
self.resetAs()
self.resetbs()
if not type(V).__name__ == 'FunctionSpace':
raise RROMPyException("V type not recognized.")
self._V = V
self.u = fen.TrialFunction(V)
self.v = fen.TestFunction(V)
def spacedim(self):
return self.V.dim()
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
self.energyNormMatrix = L2NormMatrix(self.V)
def setDirichletDatum(self, mu:complex):
"""Set Dirichlet datum if parametric."""
if hasattr(self, "liftedDirichletDatum"):
self.liftDirichletDatamu = mu
def liftDirichletData(self, mu:complex) -> Np1D:
"""Lift Dirichlet datum."""
self.setDirichletDatum(mu)
if not np.isclose(self.liftDirichletDatamu, mu):
try:
liftRe = fen.interpolate(self.DirichletDatum[0], self.V)
except:
liftRe = fen.project(self.DirichletDatum[0], self.V)
try:
liftIm = fen.interpolate(self.DirichletDatum[1], self.V)
except:
liftIm = fen.project(self.DirichletDatum[1], self.V)
self.liftedDirichletDatum = (np.array(liftRe.vector())
+ 1.j * np.array(liftIm.vector()))
return self.liftedDirichletDatum
def reduceQuadratureDegree(self, fun:FenFunc, name:str):
"""Check whether to reduce compiler parameters to degree threshold."""
if not np.isinf(self.degree_threshold):
from ufl.algorithms.estimate_degrees import (
estimate_total_polynomial_degree as ETPD)
try:
deg = ETPD(fun)
except:
return False
if deg > self.degree_threshold:
if self.verbosity >= 15:
verbosityDepth("MAIN", ("Reducing quadrature degree from "
"{} to {} for {}.").format(
deg,
self.degree_threshold,
name),
timestamp = self.timestamp)
return True
return False
def iterReduceQuadratureDegree(self, funsNames:List[Tuple[FenFunc, str]]):
"""
Iterate reduceQuadratureDegree over list and define reduce compiler
parameters.
"""
if funsNames is not None:
for fun, name in funsNames:
if self.reduceQuadratureDegree(fun, name):
return {"quadrature_degree" : self.degree_threshold}
return {}
@abstractmethod
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
if self.As[der] is None:
self.As[der] = 0.
return self.As[der]
@abstractmethod
def b(self, mu:complex, der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
bnull = self.checkbInBounds(der, homogeneized)
if bnull is not None: return bnull
b = self.bsH[der] if homogeneized else self.bs[der]
if b is None:
if homogeneized:
self.bsH[der] = 0.
else:
self.bs[der] = 0.
b = 0.
return b
def plot(self, u:Np1D, name : str = "u", save : str = None,
what : strLst = 'all', saveFormat : str = "eps",
saveDPI : int = 100, show : bool = True, **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'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
show(optional): Whether to show figure. Defaults to True.
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
plt.figure(**figspecs)
plt.jet()
if 'ABS' in what:
uAb = fen.Function(self.V)
uAb.vector().set_local(np.abs(u))
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().set_local(np.angle(u))
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().set_local(np.real(u))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fen.plot(uRe, title = "Re({0})".format(name))
plt.colorbar(p)
if 'IMAG' in what:
uIm = fen.Function(self.V)
uIm.vector().set_local(np.imag(u))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fen.plot(uIm, title = "Im({0})".format(name))
plt.colorbar(p)
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_fig_".format(save), saveFormat),
format = saveFormat, dpi = saveDPI)
if show:
plt.show()
plt.close()
def plotmesh(self, name : str = "Mesh", save : str = None,
saveFormat : str = "eps", saveDPI : int = 100,
show : bool = True, **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): Where to save plot(s). Defaults to None, i.e. no
saving.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
show(optional): Whether to show figure. Defaults to True.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
plt.figure(**figspecs)
fen.plot(self.V.mesh())
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_msh_".format(save), saveFormat),
format = saveFormat, dpi = saveDPI)
if show:
plt.show()
plt.close()
def outParaview(self, u:Np1D, name : str = "u", filename : str = "out",
time : float = 0., what : strLst = 'all',
forceNewFile : bool = True, folder : bool = False,
filePW = None):
"""
Output complex-valued function with given dofs to ParaView file.
Args:
u: numpy complex array with function dofs.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
time(optional): Timestamp.
what(optional): Which plots to do. If list, can contain 'MESH',
'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard
'ALL'. Defaults to 'ALL'.
forceNewFile(optional): Whether to create new output file.
folder(optional): Whether to create an additional folder layer.
filePW(optional): Fenics File entity (for time series).
"""
if isinstance(what, (str,)):
if what.upper() == 'ALL':
what = ['MESH', 'ABS', 'PHASE', 'REAL', 'IMAG']
else:
what = [what]
what = purgeList(what, ['MESH', 'ABS', 'PHASE', 'REAL', 'IMAG'],
listname = self.name() + ".what", baselevel = 1)
if len(what) == 0: return
if filePW is None:
if folder:
if not path.exists(filename + "/"):
mkdir(filename)
idxpath = filename.rfind("/")
filename += "/" + filename[idxpath + 1 :]
if forceNewFile:
filePW = fen.File(getNewFilename(filename, "pvd"))
else:
filePW = fen.File("{}.pvd".format(filename))
if what == ['MESH']:
filePW << (self.V.mesh(), time)
if 'ABS' in what:
uAb = fen.Function(self.V, name = "{}_ABS".format(name))
uAb.vector().set_local(np.abs(u))
filePW << (uAb, time)
if 'PHASE' in what:
uPh = fen.Function(self.V, name = "{}_PHASE".format(name))
uPh.vector().set_local(np.angle(u))
filePW << (uPh, time)
if 'REAL' in what:
uRe = fen.Function(self.V, name = "{}_REAL".format(name))
uRe.vector().set_local(np.real(u))
filePW << (uRe, time)
if 'IMAG' in what:
uIm = fen.Function(self.V, name = "{}_IMAG".format(name))
uIm.vector().set_local(np.imag(u))
filePW << (uIm, time)
return filePW
def outParaviewTimeDomain(self, u:Np1D, omega:float,
timeFinal : float = None,
periodResolution : int = 20, name : str = "u",
filename : str = "out",
forceNewFile : bool = True,
folder : bool = False):
"""
Output complex-valued function with given dofs to ParaView file,
converted to time domain.
Args:
u: numpy complex array with function dofs.
omega: frequency.
timeFinal(optional): final time of simulation.
periodResolution(optional): number of time steps per period.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
forceNewFile(optional): Whether to create new output file.
folder(optional): Whether to create an additional folder layer.
"""
if folder:
if not path.exists(filename + "/"):
mkdir(filename)
idxpath = filename.rfind("/")
filename += "/" + filename[idxpath + 1 :]
if forceNewFile:
filePW = fen.File(getNewFilename(filename, "pvd"))
else:
filePW = fen.File("{}.pvd".format(filename))
omega = np.abs(omega)
t = 0.
dt = 2. * np.pi / omega / periodResolution
if timeFinal is None: timeFinal = 2. * np.pi / omega - dt
for j in range(int(np.ceil(timeFinal / dt)) + 1):
ut = fen.Function(self.V, name = name)
ut.vector().set_local(np.real(u) * np.cos(omega * t)
+ np.imag(u) * np.sin(omega * t))
filePW << (ut, t)
t += dt
return filePW
diff --git a/rrompy/hfengines/linear_problem/helmholtz_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_problem_engine.py
index a515bc7..af81a00 100644
--- a/rrompy/hfengines/linear_problem/helmholtz_problem_engine.py
+++ b/rrompy/hfengines/linear_problem/helmholtz_problem_engine.py
@@ -1,163 +1,163 @@
# 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 .laplace_base_problem_engine import LaplaceBaseProblemEngine
from rrompy.utilities.base.types import ScOp
-from rrompy.utilities.fenics import fenZERO, fenONE
+from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.utilities.base import verbosityDepth
__all__ = ['HelmholtzProblemEngine']
class HelmholtzProblemEngine(LaplaceBaseProblemEngine):
"""
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:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
bsH: Numpy array representation of homogeneized bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
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.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
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.
dsToBeSet: Whether ds needs to be set.
"""
nAs = 2
rescalingExp = 2.
def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.omega = 1.
self.refractionIndex = fenONE
@property
def refractionIndex(self):
"""Value of n."""
return self._refractionIndex
@refractionIndex.setter
def refractionIndex(self, refractionIndex):
self.resetAs()
if not isinstance(refractionIndex, (list, tuple,)):
refractionIndex = [refractionIndex, fenZERO]
self._refractionIndex = refractionIndex
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
self.autoSetDS()
if der <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
aRe, aIm = self.diffusivity
hRe, hIm = self.RobinDatumH
termNames = ["diffusivity", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip(
[aRe, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[aIm, hIm],
[x + "Imag" for x in termNames]))
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, form_compiler_parameters = parsRe)
A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
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.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size)
+ 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
shape = A0ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a1Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a1Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
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.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
shape = A1ReMat.size)
+ 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
shape = A1ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der == 0:
return self.As[0] + mu**2 * self.As[1]
return self.As[1]
diff --git a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py
index 679bc9c..2060185 100644
--- a/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py
+++ b/rrompy/hfengines/linear_problem/helmholtz_square_bubble_domain_problem_engine.py
@@ -1,245 +1,245 @@
# 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.utilities.base.types import Np1D, ScOp, Tuple, FenExpr
-from rrompy.utilities.fenics import fenZERO
+from rrompy.solver.fenics import fenZERO
from .helmholtz_problem_engine import HelmholtzProblemEngine
from rrompy.utilities.base import verbosityDepth
__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.
"""
nAs, nbs = 3, 20
rescalingExp = 1.
def __init__(self, kappa:float, theta:float, n:int, mu0 : np.complex = 1.,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.omega = kappa
self.kappa = kappa
self.theta = theta
self.mu0 = mu0
self.forcingTermMu = np.nan
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(np.pi,np.pi), n, n)
self.V = fen.FunctionSpace(mesh, "P", 3)
def buildEnergyNormForm(self): # H1
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
mudx = np.abs(self.mu0) * fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
muM = np.abs(self.mu0) * fen.dot(self.u, self.v) * fen.dx
imudy = 1. / np.abs(self.mu0) * (fen.dot(self.u.dx(1), self.v.dx(1))
* fen.dx)
normMatFen = fen.assemble(mudx + imudy + muM)
normMat = fen.as_backend_type(normMatFen).mat()
self.energyNormMatrix = scsp.csr_matrix(normMat.getValuesCSR()[::-1],
shape = normMat.size)
def getForcingTerm(self, mu:complex) -> Tuple[FenExpr, FenExpr]:
"""Compute forcing term."""
if not np.isclose(mu, self.forcingTermMu):
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling base expression for "
"forcing term."),
timestamp = self.timestamp)
pi = np.pi
c, s = np.cos(self.theta), np.sin(self.theta)
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
muR, muI = np.real(mu), np.imag(mu)
mu2R, mu2I = np.real(mu ** 2.), np.imag(mu ** 2.)
C = 16. / pi ** 4.
bR = C * (2 * (x * (pi - x) + y * (pi - y))
+ (self.kappa * s) ** 2. * (mu2R - 1.)
* x * (pi - x) * y * (pi - y))
bI = C * (2 * self.kappa * (c * (pi - 2 * x) * y * (pi - y)
+ s * x * (pi - x) * (pi - 2 * y))
+ (self.kappa * s) ** 2. * mu2I
* x * (pi - x) * y * (pi - y))
wR = (fen.cos(self.kappa * (c * x + s * muR * y))
* fen.exp(self.kappa * s * muI * y))
wI = (fen.sin(self.kappa * (c * x + s * muR * y))
* fen.exp(self.kappa * s * muI * y))
self.forcingTerm = [bR * wR + bI * wI, bI * wR - bR * wI]
self.forcingTermMu = mu
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling base expression.",
timestamp = self.timestamp)
return self.forcingTerm
def getExtraFactorB(self, mu:complex, der:int) -> Tuple[FenExpr, FenExpr]:
"""Compute extra expression in RHS."""
def getPowMinusj(x, power):
powR = x ** power
powI = fenZERO
if power % 2 == 1:
powR, powI = powI, powR
if (power + 3) % 4 < 2:
powR, powI = - powR, - powI
return powR, powI
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling auxiliary expression for "
"forcing term derivative."),
timestamp = self.timestamp)
from scipy.special import factorial as fact
y = fen.SpatialCoordinate(self.V.mesh())[1]
powR, powI = [(self.kappa * np.sin(self.theta)) ** der * k\
for k in getPowMinusj(y, der)]
mu2R, mu2I = np.real(mu ** 2.), np.imag(mu ** 2.)
exprR = mu2R * powR - mu2I * powI
exprI = mu2I * powR + mu2R * powI
if der >= 1:
muR, muI = np.real(2. * mu), np.imag(2. * mu)
powR, powI = [(self.kappa * np.sin(self.theta)) ** (der - 1) * k\
* der for k in getPowMinusj(y, der - 1)]
exprR += muR * powR - muI * powI
exprI += muI * powR + muR * powI
if der >= 2:
powR, powI = [(self.kappa * np.sin(self.theta)) ** (der - 2) * k\
* der * (der - 1) for k in getPowMinusj(y, der - 2)]
exprR += powR
exprI += powI
fac = fact(der)
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling auxiliary expression.",
timestamp = self.timestamp)
return [exprR / fac, exprI / fac]
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
self.autoSetDS()
if der <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
A0Re = fen.assemble(a0Re)
DirichletBC0.apply(A0Re)
A0ReMat = fen.as_backend_type(A0Re).mat()
A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size,
dtype = np.complex)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der <= 2 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
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
parsRe = self.iterReduceQuadratureDegree(zip([k2n2Re],
["kappaSquaredRefractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([k2n2Im],
["kappaSquaredRefractionIndexSquaredImag"]))
a2Re = (fen.dot(self.u.dx(0), self.v.dx(0))
- k2n2Re * fen.dot(self.u, self.v)) * fen.dx
a2Im = - k2n2Im * fen.dot(self.u, self.v) * fen.dx
A2Re = fen.assemble(a2Re, form_compiler_parameters = parsRe)
A2Im = fen.assemble(a2Im, form_compiler_parameters = parsIm)
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.As[2] = (scsp.csr_matrix((A2Rev, A2Rec, A2Rer),
shape = A2ReMat.size)
+ 1.j * scsp.csr_matrix((A2Imv, A2Imc, A2Imr),
shape = A2ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der == 0:
return self.As[0] + mu ** 2 * self.As[2]
if der == 1:
return 2. * mu * self.As[2]
return self.As[2]
def b(self, mu:complex, der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
bnull = self.checkbInBounds(der, homogeneized)
if bnull is not None: return bnull
if homogeneized and not np.isclose(self.mu0BC, mu):
self.u0BC = self.liftDirichletData(mu)
if not np.isclose(self.bsmu, mu):
self.bsmu = mu
self.resetbs()
b = self.bsH[der] if homogeneized else self.bs[der]
if b is None:
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(der),
timestamp = self.timestamp)
if der < self.nbs:
fRe, fIm = self.getForcingTerm(mu)
cRe, cIm = self.getExtraFactorB(mu, der)
cfRe = cRe * fRe - cIm * fIm
cfIm = cRe * fIm + cIm * fRe
else:
cfRe, cfIm = fenZERO, fenZERO
parsRe = self.iterReduceQuadratureDegree(zip([cfRe],
["forcingTermDer{}Real".format(der)]))
parsIm = self.iterReduceQuadratureDegree(zip([cfIm],
["forcingTermDer{}Imag".format(der)]))
L0Re = fen.dot(cfRe, self.v) * fen.dx
L0Im = fen.dot(cfIm, self.v) * fen.dx
b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
if homogeneized:
Ader = self.A(mu, der)
b0Re[:] -= np.real(Ader.dot(self.u0BC))
b0Im[:] -= np.imag(Ader.dot(self.u0BC))
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
DirichletBC0.apply(b0Re)
DirichletBC0.apply(b0Im)
b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
if homogeneized:
self.bsH[der] = b
else:
self.bs[der] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return b
diff --git a/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py b/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py
index 3196353..a42eba7 100644
--- a/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py
+++ b/rrompy/hfengines/linear_problem/laplace_base_problem_engine.py
@@ -1,318 +1,318 @@
# 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.base.problem_engine_base import ProblemEngineBase
from rrompy.utilities.base.types import Np1D, ScOp
-from rrompy.utilities.fenics import fenZERO, fenONE, H1NormMatrix
+from rrompy.solver.fenics import fenZERO, fenONE, H1NormMatrix
from rrompy.utilities.base import verbosityDepth
__all__ = ['LaplaceBaseProblemEngine']
class LaplaceBaseProblemEngine(ProblemEngineBase):
"""
Solver for generic Laplace problems.
- \nabla \cdot (a \nabla u) = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
bsH: Numpy array representation of homogeneized bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
diffusivity: Value of a.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
A0: Scipy sparse array representation (in CSC format) of A0.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.omega = 0.
self.diffusivity = fenONE
self.forcingTerm = fenZERO
self.DirichletDatum = fenZERO
self.NeumannDatum = fenZERO
self.RobinDatumG = fenZERO
self.RobinDatumH = fenZERO
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
ProblemEngineBase.V.fset(self, V)
self.dsToBeSet = True
@property
def diffusivity(self):
"""Value of a."""
return self._diffusivity
@diffusivity.setter
def diffusivity(self, diffusivity):
self.resetAs()
if not isinstance(diffusivity, (list, tuple,)):
diffusivity = [diffusivity, fenZERO]
self._diffusivity = diffusivity
@property
def forcingTerm(self):
"""Value of f."""
return self._forcingTerm
@forcingTerm.setter
def forcingTerm(self, forcingTerm):
self.resetbs()
if not isinstance(forcingTerm, (list, tuple,)):
forcingTerm = [forcingTerm, fenZERO]
self._forcingTerm = forcingTerm
@property
def DirichletDatum(self):
"""Value of u0."""
return self._DirichletDatum
@DirichletDatum.setter
def DirichletDatum(self, DirichletDatum):
self.resetbs()
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):
self.resetbs()
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):
self.resetbs()
if not isinstance(RobinDatumG, (list, tuple,)):
RobinDatumG = [RobinDatumG, fenZERO]
self._RobinDatumG = RobinDatumG
@property
def RobinDatumH(self):
"""Value of h."""
return self._RobinDatumH
@RobinDatumH.setter
def RobinDatumH(self, RobinDatumH):
self.resetAs()
if not isinstance(RobinDatumH, (list, tuple,)):
RobinDatumH = [RobinDatumH, fenZERO]
self._RobinDatumH = RobinDatumH
@property
def DirichletBoundary(self):
"""Function handle to DirichletBoundary."""
return self.BCManager.DirichletBoundary
@DirichletBoundary.setter
def DirichletBoundary(self, DirichletBoundary):
self.resetAs()
self.resetbs()
self.BCManager.DirichletBoundary = DirichletBoundary
@property
def NeumannBoundary(self):
"""Function handle to NeumannBoundary."""
return self.BCManager.NeumannBoundary
@NeumannBoundary.setter
def NeumannBoundary(self, NeumannBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.NeumannBoundary = NeumannBoundary
@property
def RobinBoundary(self):
"""Function handle to RobinBoundary."""
return self.BCManager.RobinBoundary
@RobinBoundary.setter
def RobinBoundary(self, RobinBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.RobinBoundary = RobinBoundary
def autoSetDS(self):
"""Set FEniCS boundary measure based on boundary function handles."""
if self.dsToBeSet:
if self.verbosity >= 20:
verbosityDepth("INIT", "Initializing boundary measures.",
timestamp = self.timestamp)
mesh = self.V.mesh()
NB = self.NeumannBoundary
RB = self.RobinBoundary
boundary_markers = fen.MeshFunction("size_t", mesh,
mesh.topology().dim() - 1)
NB.mark(boundary_markers, 0)
RB.mark(boundary_markers, 1)
self.ds = fen.Measure("ds", domain = mesh,
subdomain_data = boundary_markers)
self.dsToBeSet = False
if self.verbosity >= 20:
verbosityDepth("DEL", "Done initializing boundary measures.",
timestamp = self.timestamp)
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
self.energyNormMatrix = H1NormMatrix(self.V, np.abs(self.omega)**2)
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
self.autoSetDS()
if self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
aRe, aIm = self.diffusivity
hRe, hIm = self.RobinDatumH
termNames = ["diffusivity", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip(
[aRe, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[aIm, hIm],
[x + "Imag" for x in termNames]))
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, form_compiler_parameters = parsRe)
A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
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.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size)
+ 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
shape = A0ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self.As[0]
def b(self, mu:complex, der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
bnull = self.checkbInBounds(der, homogeneized)
if bnull is not None: return bnull
if homogeneized and not np.isclose(self.mu0BC, mu):
self.u0BC = self.liftDirichletData(mu)
if (max(self.nbs, self.nAs * homogeneized) > 1
and not np.isclose(self.bsmu, mu)):
self.bsmu = mu
self.resetbs()
b = self.bsH[der] if homogeneized else self.bs[der]
if b is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(der),
timestamp = self.timestamp)
if der == 0:
fRe, fIm = self.forcingTerm
g1Re, g1Im = self.NeumannDatum
g2Re, g2Im = self.RobinDatumG
else:
fRe, fIm = fenZERO, fenZERO
g1Re, g1Im = fenZERO, fenZERO
g2Re, g2Im = fenZERO, fenZERO
termNames = ["forcingTerm", "NeumannDatum", "RobinDatumG"]
parsRe = self.iterReduceQuadratureDegree(zip(
[fRe, g1Re, g2Re],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[fIm, g1Im, g2Im],
[x + "Imag" for x in termNames]))
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, form_compiler_parameters = parsRe)
b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
if homogeneized:
Ader = self.A(mu, der)
b0Re[:] -= np.real(Ader.dot(self.u0BC))
b0Im[:] -= np.imag(Ader.dot(self.u0BC))
DBCR = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary)
else:
DBCR = fen.DirichletBC(self.V, self.DirichletDatum[0],
self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, self.DirichletDatum[1],
self.DirichletBoundary)
DBCR.apply(b0Re)
DBCI.apply(b0Im)
b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
if homogeneized:
self.bsH[der] = b
else:
self.bs[der] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return b
diff --git a/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py b/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py
index e65f308..f1146ef 100644
--- a/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py
+++ b/rrompy/hfengines/linear_problem/laplace_disk_gaussian.py
@@ -1,158 +1,158 @@
# 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.utilities.base.types import Np1D, Tuple, FenExpr
from .laplace_base_problem_engine import LaplaceBaseProblemEngine
-from rrompy.utilities.fenics import fenZERO, fenONE
+from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.utilities.base import verbosityDepth
__all__ = ['LaplaceDiskGaussian']
class LaplaceDiskGaussian(LaplaceBaseProblemEngine):
"""
Solver for disk Laplace problems with parametric forcing term center.
- \Delta u = C exp(-.5 * ||\cdot - (mu, 0)||^2) in \Omega = B(0, 5)
u = 0 on \partial\Omega.
"""
nbs = 20
def __init__(self, n:int, degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.computebsFactors()
self.forcingTermMu = np.nan
import mshr
mesh = mshr.generate_mesh(mshr.Circle(fen.Point(0., 0.), 5.), n)
self.V = fen.FunctionSpace(mesh, "P", 3)
def getForcingTerm(self, mu:complex) -> Tuple[FenExpr, FenExpr]:
"""Compute forcing term."""
if not np.isclose(mu, self.forcingTermMu):
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling base expression for "
"forcing term."),
timestamp = self.timestamp)
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
C = np.exp(-.5 * mu ** 2.)
CR, CI = np.real(C), np.imag(C)
f0 = (2 * np.pi) ** -.5 * fen.exp(-.5 * (x ** 2. + y ** 2.))
muR, muI = np.real(mu), np.imag(mu)
f1R = fen.exp(muR * x) * fen.cos(muI * x)
f1I = fen.exp(muR * x) * fen.sin(muI * x)
self.forcingTerm = [f0 * (CR * f1R - CI * f1I),
f0 * (CR * f1I + CI * f1R)]
self.forcingTermMu = mu
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling base expression.",
timestamp = self.timestamp)
return self.forcingTerm
def computebsFactors(self):
self.bsFactors = np.zeros((self.nbs, self.nbs), dtype = float)
self.bsFactors[0, 0] = 1.
self.bsFactors[1, 1] = 1.
for j in range(2, self.nbs):
l = (j + 1) % 2 + 1
J = np.arange(l, j + 1, 2)
self.bsFactors[j, J] = self.bsFactors[j - 1, J - 1]
if l == 2:
l = 0
J = np.arange(l, j, 2)
self.bsFactors[j, J] += np.multiply(- 1 - J,
self.bsFactors[j - 1, J + 1])
self.bsFactors[j, l : j + 2 : 2] /= j
def getExtraFactorB(self, mu:complex, der:int) -> Tuple[FenExpr, FenExpr]:
"""Compute extra expression in RHS."""
if self.verbosity >= 25:
verbosityDepth("INIT", ("Assembling auxiliary expression for "
"forcing term derivative."),
timestamp = self.timestamp)
muR, muI = np.real(mu), np.imag(mu)
x = fen.SpatialCoordinate(self.V.mesh())[0]
l = der % 2
if l == 0:
powR, powI = fenONE, fenZERO
else:
powR, powI = x - muR, fen.Constant(muI)
exprR, exprI = [self.bsFactors[der, l] * k for k in [powR, powI]]
for j in range(l + 2, der + 1, 2):
for _ in range(2):
powR, powI = (powR * (x - muR) - powI * muI,
powR * muI + powI * (x - muR))
exprR += self.bsFactors[der, j] * powR
exprI += self.bsFactors[der, j] * powI
if self.verbosity >= 25:
verbosityDepth("DEL", "Done assembling auxiliary expression.",
timestamp = self.timestamp)
return[exprR, exprI]
def b(self, mu:complex, der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
bnull = self.checkbInBounds(der, homogeneized)
if bnull is not None: return bnull
if homogeneized and not np.isclose(self.mu0BC, mu):
self.u0BC = self.liftDirichletData(mu)
if not np.isclose(self.bsmu, mu):
self.bsmu = mu
self.resetbs()
b = self.bsH[der] if homogeneized else self.bs[der]
if b is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling forcing term b{}.".format(
der),
timestamp = self.timestamp)
if der < self.nbs:
fRe, fIm = self.getForcingTerm(mu)
cRe, cIm = self.getExtraFactorB(mu, der)
cfRe = cRe * fRe - cIm * fIm
cfIm = cRe * fIm + cIm * fRe
else:
cfRe, cfIm = fenZERO, fenZERO
parsRe = self.iterReduceQuadratureDegree(zip([cfRe],
["forcingTermDer{}Real".format(der)]))
parsIm = self.iterReduceQuadratureDegree(zip([cfIm],
["forcingTermDer{}Imag".format(der)]))
L0Re = fen.dot(cfRe, self.v) * fen.dx
L0Im = fen.dot(cfIm, self.v) * fen.dx
b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
if homogeneized:
Ader = self.A(mu, der)
b0Re[:] -= np.real(Ader.dot(self.u0BC))
b0Im[:] -= np.imag(Ader.dot(self.u0BC))
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
DirichletBC0.apply(b0Re)
DirichletBC0.apply(b0Im)
b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
if homogeneized:
self.bsH[der] = b
else:
self.bs[der] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return b
diff --git a/rrompy/hfengines/linear_problem/scattering_problem_engine.py b/rrompy/hfengines/linear_problem/scattering_problem_engine.py
index 88f82d4..8345ac3 100644
--- a/rrompy/hfengines/linear_problem/scattering_problem_engine.py
+++ b/rrompy/hfengines/linear_problem/scattering_problem_engine.py
@@ -1,177 +1,177 @@
# 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 numpy import inf
import scipy.sparse as scsp
import fenics as fen
from rrompy.utilities.base.types import ScOp
-from rrompy.utilities.fenics import fenZERO
+from rrompy.solver.fenics import fenZERO
from rrompy.utilities.base import verbosityDepth
from .helmholtz_problem_engine import HelmholtzProblemEngine
from rrompy.utilities.exception_manager import RROMPyWarning
__all__ = ['ScatteringProblemEngine']
class ScatteringProblemEngine(HelmholtzProblemEngine):
"""
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 omega u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
bsH: Numpy array representation of homogeneized bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
signR: Sign in ABC.
omega: Value of omega.
diffusivity: Value of a.
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.
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 A2.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
nAs = 3
rescalingExp = 1.
signR = - 1.
def __init__(self, degree_threshold : int = inf, verbosity : int = 10,
timestamp : bool = True):
self.silenceWarnings = True
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
del self.silenceWarnings
@property
def RobinDatumH(self):
"""Value of h."""
return self.signR * self.omega
@RobinDatumH.setter
def RobinDatumH(self, RobinDatumH):
if not hasattr(self, "silenceWarnings"):
RROMPyWarning(("Scattering problems do not allow changes of h. "
"Ignoring assignment."))
return
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
self.autoSetDS()
if der <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
aRe, aIm = self.diffusivity
parsRe = self.iterReduceQuadratureDegree(zip([aRe],
["diffusivityReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([aIm],
["diffusivityImag"]))
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, form_compiler_parameters = parsRe)
A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
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.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size)
+ 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
shape = A0ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
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.As[1] = self.signR * 1.j * scsp.csr_matrix((A1v, A1c, A1r),
shape = A1Mat.size)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der <= 2 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a2Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a2Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
A2Re = fen.assemble(a2Re, form_compiler_parameters = parsRe)
A2Im = fen.assemble(a2Im, form_compiler_parameters = parsIm)
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.As[2] = (scsp.csr_matrix((A2Rev, A2Rec, A2Rer),
shape = A2ReMat.size)
+ 1.j * scsp.csr_matrix((A2Imv, A2Imc, A2Imr),
shape = A2ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der == 0:
return self.As[0] + mu * self.As[1] + mu**2. * self.As[2]
if der == 1:
return self.As[1] + 2 * mu * self.As[2]
return self.As[2]
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py
index de82104..96d7701 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py
+++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_beam_poisson_ratio.py
@@ -1,151 +1,151 @@
# 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 .linear_elasticity_problem_engine import LinearElasticityProblemEngine
-from rrompy.utilities.fenics import fenZEROS
+from rrompy.solver.fenics import fenZEROS
from rrompy.utilities.base.types import Np1D, ScOp
from rrompy.utilities.base import verbosityDepth
__all__ = ['LinearElasticityBeamPoissonRatio']
class LinearElasticityBeamPoissonRatio(LinearElasticityProblemEngine):
"""
Solver for linear elasticity problem of a beam subject to its own weight,
with parametric Poisson's ratio.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u)) = rho_ * g in \Omega
u = 0 on \Gamma_D
\partial_nu = 0 on \Gamma_N
"""
nAs = 2
nbs = 3
def __init__(self, n:int, rho_:float, g:float, E:float, nu0:float,
length:float, degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.lambda_ = E * nu0 / (1. + nu0) / (1. - 2 * nu0)
self.mu_ = E / (1. + nu0)
mesh = fen.RectangleMesh(fen.Point(0., 0.), fen.Point(length, 1),
n, max(int(n / length), 1))
self.V = fen.VectorFunctionSpace(mesh, "P", 1)
self.forcingTerm = [fen.Constant((0., - rho_ * g / E)), fenZEROS(2)]
self.DirichletBoundary = lambda x, on_b: on_b and fen.near(x[0], 0.)
self.NeumannBoundary = "REST"
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
self.autoSetDS()
if der <= 1 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
epsilon = lambda u: .5 * (fen.grad(u) + fen.nabla_grad(u))
a0Re = 2 * fen.inner(epsilon(self.u), epsilon(self.v)) * fen.dx
A0Re = fen.assemble(a0Re)
DirichletBC0.apply(A0Re)
A0ReMat = fen.as_backend_type(A0Re).mat()
A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
self.As[0] = scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size,
dtype = np.complex)
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V, fenZEROS(2),
self.DirichletBoundary)
a1Re = fen.div(self.u) * fen.div(self.v) * fen.dx
A1Re = fen.assemble(a1Re)
DirichletBC0.apply(A1Re)
A1ReMat = fen.as_backend_type(A1Re).mat()
A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
self.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
shape = A1ReMat.size,
dtype = np.complex)
- 2. * self.As[0])
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der == 0:
return self.As[0] + mu * self.As[1]
return self.As[1]
def b(self, mu:complex, der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
assert homogeneized == False
bnull = self.checkbInBounds(der)
if bnull is not None: return bnull
if (self.nbs > 1 and not np.isclose(self.bsmu, mu)):
self.bsmu = mu
self.resetbs()
b = self.bsH[der] if homogeneized else self.bs[der]
if b is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(der),
timestamp = self.timestamp)
if self.bs[0] is None and der > 0: self.b(mu, 0)
if der == 0:
fRe, fIm = self.forcingTerm
parsRe = self.iterReduceQuadratureDegree(zip(
[fRe],
["forcingTermReal"]))
parsIm = self.iterReduceQuadratureDegree(zip(
[fIm],
["forcingTermImag"]))
L0Re = fen.inner(fRe, self.v) * fen.dx
L0Im = fen.inner(fIm, self.v) * fen.dx
b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
DBCR = fen.DirichletBC(self.V, self.DirichletDatum[0],
self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, self.DirichletDatum[1],
self.DirichletBoundary)
DBCR.apply(b0Re)
DBCI.apply(b0Im)
b = (1. - mu - 2 * mu ** 2) * np.array(b0Re[:] + 1.j * b0Im[:],
dtype = np.complex)
elif der == 1:
b = (- 1. - 4 * mu) / (1. - mu - 2 * mu ** 2) * self.bs[0]
elif der == 2:
b = - 2. / (1. - mu - 2 * mu ** 2) * self.bs[0]
if homogeneized:
self.bsH[der] = b
else:
self.bs[der] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return b
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_archway_frequency.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_archway_frequency.py
index a143faa..fa8ffef 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_archway_frequency.py
+++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_archway_frequency.py
@@ -1,67 +1,67 @@
# 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 .linear_elasticity_helmholtz_problem_engine import \
LinearElasticityHelmholtzProblemEngine
-from rrompy.utilities.fenics import fenZEROS
+from rrompy.solver.fenics import fenZEROS
__all__ = ['LinearElasticityHelmholtzArchwayFrequency']
class LinearElasticityHelmholtzArchwayFrequency(
LinearElasticityHelmholtzProblemEngine):
"""
Solver for archway linear elasticity Helmholtz problem with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * omega^2 * u = rho_ * g / omega in \Omega
u = 0 on \Gamma_D
\partial_nu = 0 on \Gamma_N
"""
def __init__(self, kappa:float, n:int, rho_:float, T:float, lambda_:float,
mu_:float, R:float, r:float, degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.omega = kappa
self.lambda_ = lambda_
self.mu_ = mu_
self.rho_ = rho_
import mshr
domain = (mshr.Circle(fen.Point(0, 0), R)
- mshr.Circle(fen.Point(0, 0), r)
- mshr.Rectangle(fen.Point(-1.05*R, -1.05*R),
fen.Point(1.05*R, 0)))
mesh = mshr.generate_mesh(domain, n)
self.V = fen.VectorFunctionSpace(mesh, "P", 1)
import ufl
x, y = fen.SpatialCoordinate(mesh)[:]
NeumannNonZero = ufl.And(ufl.gt(y, r),
ufl.And(ufl.ge(x, -.25 * R), ufl.le(x, .25 * R)))
self.NeumannDatum = [ufl.as_vector((0.,
ufl.conditional(NeumannNonZero,
fen.Constant(T),
0.))),
fenZEROS(2)]
self.DirichletBoundary = lambda x, on_b: on_b and fen.near(x[1], 0.)
self.NeumannBoundary = "REST"
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine.py
index a31d6ec..0e543d8 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine.py
+++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine.py
@@ -1,183 +1,182 @@
# 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 scipy.sparse import csr_matrix
import fenics as fen
from .linear_elasticity_problem_engine import LinearElasticityProblemEngine
from rrompy.utilities.base.types import ScOp
-from rrompy.utilities.fenics import (fenZERO, fenZEROS, fenONE,
- elasticNormMatrix)
+from rrompy.solver.fenics import fenZERO, fenZEROS, fenONE, elasticNormMatrix
from rrompy.utilities.base import verbosityDepth
__all__ = ['LinearElasticityHelmholtzProblemEngine']
class LinearElasticityHelmholtzProblemEngine(LinearElasticityProblemEngine):
"""
Solver for generic linear elasticity Helmholtz problems with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * mu^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:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real vector FE space.
u: Generic vector trial functions for variational form evaluation.
v: Generic vector test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
lambda_: Value of lambda_.
mu_: Value of mu_.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
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.
dsToBeSet: Whether ds needs to be set.
"""
nAs = 2
rescalingExp = 2.
def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.omega = 1.
self.rho_ = fenONE
@property
def rho_(self):
"""Value of rho_."""
return self._rho_
@rho_.setter
def rho_(self, rho_):
self.resetAs()
if not isinstance(rho_, (list, tuple,)):
rho_ = [rho_, fenZERO]
self._rho_ = rho_
def buildEnergyNormForm(self): # energy + omega norm
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
lambda_Re, _ = self.lambda_
mu_Re, _ = self.mu_
r_Re, _ = self.rho_
self.energyNormMatrix = elasticNormMatrix(self.V, lambda_Re, mu_Re,
np.abs(self.omega)**2 * r_Re)
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
self.autoSetDS()
if der <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
lambda_Re, lambda_Im = self.lambda_
mu_Re, mu_Im = self.mu_
hRe, hIm = self.RobinDatumH
termNames = ["lambda_", "mu_", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip(
[lambda_Re, mu_Re, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[lambda_Im, mu_Re, hIm],
[x + "Imag" for x in termNames]))
epsilon = lambda u: 0.5 * (fen.grad(u) + fen.nabla_grad(u))
sigma = lambda u, l_, m_: (
l_ * fen.div(u) * fen.Identity(u.geometric_dimension())
+ 2. * m_ * epsilon(u))
a0Re = (fen.inner(sigma(self.u, lambda_Re, mu_Re),
epsilon(self.v)) * fen.dx
+ hRe * fen.inner(self.u, self.v) * self.ds(1))
a0Im = (fen.inner(sigma(self.u, lambda_Im, mu_Im),
epsilon(self.v)) * fen.dx
+ hIm * fen.inner(self.u, self.v) * self.ds(1))
A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
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.As[0] = (csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size)
+ 1.j * csr_matrix((A0Imv, A0Imc, A0Imr),
shape = A0ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
rho_Re, rho_Im = self.rho_
parsRe = self.iterReduceQuadratureDegree(zip([rho_Re],
["rho_Real"]))
parsIm = self.iterReduceQuadratureDegree(zip([rho_Im],
["rho_Imag"]))
a1Re = - rho_Re * fen.inner(self.u, self.v) * fen.dx
a1Im = - rho_Im * fen.inner(self.u, self.v) * fen.dx
A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
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.As[1] = (csr_matrix((A1Rev, A1Rec, A1Rer),
shape = A1ReMat.size)
+ 1.j * csr_matrix((A1Imv, A1Imc, A1Imr),
shape = A1ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der == 0:
return self.As[0] + mu**2 * self.As[1]
return self.As[1]
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine_damped.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine_damped.py
index fcc6bf2..a1b3603 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine_damped.py
+++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_helmholtz_problem_engine_damped.py
@@ -1,209 +1,209 @@
# 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 scipy.sparse import csr_matrix
import fenics as fen
from .linear_elasticity_helmholtz_problem_engine import \
LinearElasticityHelmholtzProblemEngine
from rrompy.utilities.base.types import ScOp
-from rrompy.utilities.fenics import fenZERO, fenZEROS
+from rrompy.solver.fenics import fenZERO, fenZEROS
from rrompy.utilities.base import verbosityDepth
__all__ = ['LinearElasticityHelmholtzProblemEngineDamped']
class LinearElasticityHelmholtzProblemEngineDamped(
LinearElasticityHelmholtzProblemEngine):
"""
Solver for generic linear elasticity Helmholtz problems with parametric
wavenumber.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u))
- rho_ * (mu^2 - i * eta * mu) * u = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real vector FE space.
u: Generic vector trial functions for variational form evaluation.
v: Generic vector test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
omega: Value of omega.
lambda_: Value of lambda_.
mu_: Value of mu_.
eta: Value of eta.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
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.
dsToBeSet: Whether ds needs to be set.
"""
nAs = 3
rescalingExp = 1.
def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.eta = fenZERO
@property
def eta(self):
"""Value of eta."""
return self._eta
@eta.setter
def eta(self, eta):
self.resetAs()
if not isinstance(eta, (list, tuple,)):
eta = [eta, fenZERO]
self._eta = eta
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
self.autoSetDS()
if der <= 0 and self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
lambda_Re, lambda_Im = self.lambda_
mu_Re, mu_Im = self.mu_
hRe, hIm = self.RobinDatumH
termNames = ["lambda_", "mu_", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip(
[lambda_Re, mu_Re, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[lambda_Im, mu_Re, hIm],
[x + "Imag" for x in termNames]))
epsilon = lambda u: 0.5 * (fen.grad(u) + fen.nabla_grad(u))
sigma = lambda u, l_, m_: (
l_ * fen.div(u) * fen.Identity(u.geometric_dimension())
+ 2. * m_ * epsilon(u))
a0Re = (fen.inner(sigma(self.u, lambda_Re, mu_Re),
epsilon(self.v)) * fen.dx
+ hRe * fen.inner(self.u, self.v) * self.ds(1))
a0Im = (fen.inner(sigma(self.u, lambda_Im, mu_Im),
epsilon(self.v)) * fen.dx
+ hIm * fen.inner(self.u, self.v) * self.ds(1))
A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
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.As[0] = (csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size)
+ 1.j * csr_matrix((A0Imv, A0Imc, A0Imr),
shape = A0ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der <= 1 and self.As[1] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A1.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
rho_Re, rho_Im = self.rho_
eta_Re, eta_Im = self.eta
termNames = ["rho_", "eta"]
parsRe = self.iterReduceQuadratureDegree(zip([rho_Re, eta_Re],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip([rho_Im, eta_Im],
[x + "Imag" for x in termNames]))
a1Re = - ((eta_Re * rho_Im + eta_Im * rho_Re)
* fen.inner(self.u, self.v)) * fen.dx
a1Im = ((eta_Re * rho_Re - eta_Im * rho_Im)
* fen.inner(self.u, self.v)) * fen.dx
A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
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.As[1] = (csr_matrix((A1Rev, A1Rec, A1Rer),
shape = A1ReMat.size)
+ 1.j * csr_matrix((A1Imv, A1Imc, A1Imr),
shape = A1ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der <= 2 and self.As[2] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A2.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
rho_Re, rho_Im = self.rho_
parsRe = self.iterReduceQuadratureDegree(zip([rho_Re],
["rho_Real"]))
parsIm = self.iterReduceQuadratureDegree(zip([rho_Im],
["rho_Imag"]))
a2Re = - rho_Re * fen.inner(self.u, self.v) * fen.dx
a2Im = - rho_Im * fen.inner(self.u, self.v) * fen.dx
A2Re = fen.assemble(a2Re, form_compiler_parameters = parsRe)
A2Im = fen.assemble(a2Im, form_compiler_parameters = parsIm)
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.As[2] = (csr_matrix((A2Rev, A2Rec, A2Rer),
shape = A2ReMat.size)
+ 1.j * csr_matrix((A2Imv, A2Imc, A2Imr),
shape = A2ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
if der == 0:
return self.As[0] + mu * self.As[1] + mu**2 * self.As[2]
if der == 1:
return self.As[1] + 2 * mu * self.As[2]
return self.As[2]
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py
index 1efc734..be8d5b1 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py
+++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py
@@ -1,354 +1,353 @@
# 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 scipy.sparse import csr_matrix
import fenics as fen
from rrompy.hfengines.base.vector_problem_engine_base import \
VectorProblemEngineBase
from rrompy.utilities.base.types import Np1D, ScOp
-from rrompy.utilities.fenics import (fenZERO, fenZEROS, fenONE,
- elasticNormMatrix)
+from rrompy.solver.fenics import fenZERO, fenZEROS, fenONE, elasticNormMatrix
from rrompy.utilities.base import verbosityDepth
__all__ = ['LinearElasticityProblemEngine']
class LinearElasticityProblemEngine(VectorProblemEngineBase):
"""
Solver for generic linear elasticity problems.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u)) = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real vector FE space.
u: Generic vector trial functions for variational form evaluation.
v: Generic vector test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
bsmu: Mu value of last bs evaluation.
liftDirichletDatamu: Mu value of last Dirichlet datum evaluation.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
lambda_: Value of lambda_.
mu_: Value of mu_.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
A0: Scipy sparse array representation (in CSC format) of A0.
b0: Numpy array representation of b0.
dsToBeSet: Whether ds needs to be set.
"""
def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.lambda_ = fenONE
self.mu_ = fenONE
self.forcingTerm = fenZEROS(self.V.mesh().topology().dim())
self.DirichletDatum = fenZEROS(self.V.mesh().topology().dim())
self.NeumannDatum = fenZEROS(self.V.mesh().topology().dim())
self.RobinDatumG = fenZEROS(self.V.mesh().topology().dim())
self.RobinDatumH = fenZERO
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
VectorProblemEngineBase.V.fset(self, V)
self.forcingTerm = fenZEROS(self.V.mesh().topology().dim())
self.DirichletDatum = fenZEROS(self.V.mesh().topology().dim())
self.NeumannDatum = fenZEROS(self.V.mesh().topology().dim())
self.RobinDatumG = fenZEROS(self.V.mesh().topology().dim())
self.dsToBeSet = True
@property
def lambda_(self):
"""Value of lambda_."""
return self._lambda_
@lambda_.setter
def lambda_(self, lambda_):
self.resetAs()
if not isinstance(lambda_, (list, tuple,)):
lambda_ = [lambda_, fenZERO]
self._lambda_ = lambda_
@property
def mu_(self):
"""Value of mu_."""
return self._mu_
@mu_.setter
def mu_(self, mu_):
self.resetAs()
if not isinstance(mu_, (list, tuple,)):
mu_ = [mu_, fenZERO]
self._mu_ = mu_
@property
def forcingTerm(self):
"""Value of f."""
return self._forcingTerm
@forcingTerm.setter
def forcingTerm(self, forcingTerm):
self.resetbs()
if not isinstance(forcingTerm, (list, tuple,)):
forcingTerm = [forcingTerm,
fenZEROS(self.V.mesh().topology().dim())]
self._forcingTerm = forcingTerm
@property
def DirichletDatum(self):
"""Value of u0."""
return self._DirichletDatum
@DirichletDatum.setter
def DirichletDatum(self, DirichletDatum):
self.resetbs()
if not isinstance(DirichletDatum, (list, tuple,)):
DirichletDatum = [DirichletDatum,
fenZEROS(self.V.mesh().topology().dim())]
self._DirichletDatum = DirichletDatum
@property
def NeumannDatum(self):
"""Value of g1."""
return self._NeumannDatum
@NeumannDatum.setter
def NeumannDatum(self, NeumannDatum):
self.resetbs()
if not isinstance(NeumannDatum, (list, tuple,)):
NeumannDatum = [NeumannDatum,
fenZEROS(self.V.mesh().topology().dim())]
self._NeumannDatum = NeumannDatum
@property
def RobinDatumG(self):
"""Value of g2."""
return self._RobinDatumG
@RobinDatumG.setter
def RobinDatumG(self, RobinDatumG):
self.resetbs()
if not isinstance(RobinDatumG, (list, tuple,)):
RobinDatumG = [RobinDatumG,
fenZEROS(self.V.mesh().topology().dim())]
self._RobinDatumG = RobinDatumG
@property
def RobinDatumH(self):
"""Value of h."""
return self._RobinDatumH
@RobinDatumH.setter
def RobinDatumH(self, RobinDatumH):
self.resetAs()
if not isinstance(RobinDatumH, (list, tuple,)):
RobinDatumH = [RobinDatumH, fenZERO]
self._RobinDatumH = RobinDatumH
@property
def DirichletBoundary(self):
"""Function handle to DirichletBoundary."""
return self.BCManager.DirichletBoundary
@DirichletBoundary.setter
def DirichletBoundary(self, DirichletBoundary):
self.resetAs()
self.resetbs()
self.BCManager.DirichletBoundary = DirichletBoundary
@property
def NeumannBoundary(self):
"""Function handle to NeumannBoundary."""
return self.BCManager.NeumannBoundary
@NeumannBoundary.setter
def NeumannBoundary(self, NeumannBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.NeumannBoundary = NeumannBoundary
@property
def RobinBoundary(self):
"""Function handle to RobinBoundary."""
return self.BCManager.RobinBoundary
@RobinBoundary.setter
def RobinBoundary(self, RobinBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.RobinBoundary = RobinBoundary
def autoSetDS(self):
"""Set FEniCS boundary measure based on boundary function handles."""
if self.dsToBeSet:
if self.verbosity >= 20:
verbosityDepth("INIT", "Initializing boundary measures.",
timestamp = self.timestamp)
NB = self.NeumannBoundary
RB = self.RobinBoundary
boundary_markers = fen.MeshFunction("size_t", self.V.mesh(),
self.V.mesh().topology().dim() - 1)
NB.mark(boundary_markers, 0)
RB.mark(boundary_markers, 1)
self.ds = fen.Measure("ds", domain = self.V.mesh(),
subdomain_data = boundary_markers)
self.dsToBeSet = False
if self.verbosity >= 20:
verbosityDepth("DEL", "Done initializing boundary measures.",
timestamp = self.timestamp)
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
lambda_Re, _ = self.lambda_
mu_Re, _ = self.mu_
self.energyNormMatrix = elasticNormMatrix(self.V, lambda_Re, mu_Re)
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
self.autoSetDS()
if self.As[0] is None:
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling operator term A0.",
timestamp = self.timestamp)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
lambda_Re, lambda_Im = self.lambda_
mu_Re, mu_Im = self.mu_
hRe, hIm = self.RobinDatumH
termNames = ["lambda_", "mu_", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip(
[lambda_Re, mu_Re, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[lambda_Im, mu_Re, hIm],
[x + "Imag" for x in termNames]))
epsilon = lambda u: 0.5 * (fen.grad(u) + fen.nabla_grad(u))
sigma = lambda u, l_, m_: (
l_ * fen.div(u) * fen.Identity(u.geometric_dimension())
+ 2. * m_ * epsilon(u))
a0Re = (fen.inner(sigma(self.u, lambda_Re, mu_Re),
epsilon(self.v)) * fen.dx
+ hRe * fen.inner(self.u, self.v) * self.ds(1))
a0Im = (fen.inner(sigma(self.u, lambda_Im, mu_Im),
epsilon(self.v)) * fen.dx
+ hIm * fen.inner(self.u, self.v) * self.ds(1))
A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
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.As[0] = (csr_matrix((A0Rev, A0Rec, A0Rer),
shape = A0ReMat.size)
+ 1.j * csr_matrix((A0Imv, A0Imc, A0Imr),
shape = A0ImMat.size))
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling operator term.",
timestamp = self.timestamp)
return self.As[0]
def b(self, mu:complex, der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
bnull = self.checkbInBounds(der, homogeneized)
if bnull is not None: return bnull
if homogeneized and not np.isclose(self.mu0BC, mu):
self.u0BC = self.liftDirichletData(mu)
if (max(self.nbs, self.nAs * homogeneized) > 1
and not np.isclose(self.bsmu, mu)):
self.bsmu = mu
self.resetbs()
b = self.bsH[der] if homogeneized else self.bs[der]
if b is None:
self.autoSetDS()
if self.verbosity >= 20:
verbosityDepth("INIT", ("Assembling forcing term "
"b{}.").format(der),
timestamp = self.timestamp)
if der == 0:
fRe, fIm = self.forcingTerm
g1Re, g1Im = self.NeumannDatum
g2Re, g2Im = self.RobinDatumG
else:
fRe = fenZEROS(self.V.mesh().topology().dim())
fIm = fenZEROS(self.V.mesh().topology().dim())
g1Re = fenZEROS(self.V.mesh().topology().dim())
g1Im = fenZEROS(self.V.mesh().topology().dim())
g2Re = fenZEROS(self.V.mesh().topology().dim())
g2Im = fenZEROS(self.V.mesh().topology().dim())
termNames = ["forcingTerm", "NeumannDatum", "RobinDatumG"]
parsRe = self.iterReduceQuadratureDegree(zip(
[fRe, g1Re, g2Re],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[fIm, g1Im, g2Im],
[x + "Imag" for x in termNames]))
L0Re = (fen.inner(fRe, self.v) * fen.dx
+ fen.inner(g1Re, self.v) * self.ds(0)
+ fen.inner(g2Re, self.v) * self.ds(1))
L0Im = (fen.inner(fIm, self.v) * fen.dx
+ fen.inner(g1Im, self.v) * self.ds(0)
+ fen.inner(g2Im, self.v) * self.ds(1))
b0Re = fen.assemble(L0Re, form_compiler_parameters = parsRe)
b0Im = fen.assemble(L0Im, form_compiler_parameters = parsIm)
if homogeneized:
Ader = self.A(mu, der)
b0Re[:] -= np.real(Ader.dot(self.u0BC))
b0Im[:] -= np.imag(Ader.dot(self.u0BC))
DBCR = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
else:
DBCR = fen.DirichletBC(self.V, self.DirichletDatum[0],
self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, self.DirichletDatum[1],
self.DirichletBoundary)
DBCR.apply(b0Re)
DBCI.apply(b0Im)
b = np.array(b0Re[:] + 1.j * b0Im[:], dtype = np.complex)
if homogeneized:
self.bsH[der] = b
else:
self.bs[der] = b
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling forcing term.",
timestamp = self.timestamp)
return b
diff --git a/rrompy/reduction_methods/distributed_greedy/generic_distributed_greedy_approximant.py b/rrompy/reduction_methods/distributed_greedy/generic_distributed_greedy_approximant.py
index 0d938bd..f7d6d4a 100644
--- a/rrompy/reduction_methods/distributed_greedy/generic_distributed_greedy_approximant.py
+++ b/rrompy/reduction_methods/distributed_greedy/generic_distributed_greedy_approximant.py
@@ -1,651 +1,641 @@
# 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 deepcopy as copy
import numpy as np
from rrompy.reduction_methods.distributed.generic_distributed_approximant \
import GenericDistributedApproximant
-from rrompy.utilities.base.types import Np1D, Np2D, DictAny, HFEng, Tuple, List
+from rrompy.utilities.base.types import (Np1D, Np2D, DictAny, HFEng, Tuple,
+ List, normEng)
from rrompy.utilities.base import purgeDict, verbosityDepth
+from rrompy.solver import normEngine
from rrompy.utilities.exception_manager import (RROMPyException, modeAssert,
RROMPyWarning)
__all__ = ['GenericDistributedGreedyApproximant']
-class normEngine:
- def __init__(self, energyNormMatrix:Np2D):
- self.energyNormMatrix = copy(energyNormMatrix)
-
- def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False) -> Np2D:
- if onlyDiag:
- return np.sum(self.energyNormMatrix.dot(u) * v.conj(), axis = 0)
- return v.T.conj().dot(self.energyNormMatrix.dot(u))
-
- def norm(self, u:Np2D) -> Np1D:
- return np.abs(self.innerProduct(u, u, onlyDiag = True)) ** .5
-
def pruneSamples(mus:Np1D, badmus:Np1D, tol : float = 1e-8) -> Np1D:
"""Remove from mus all the elements which are too close to badmus."""
proximity = np.min(np.abs(mus.reshape(-1, 1)
- np.tile(badmus.reshape(1, -1), [len(mus), 1])),
axis = 1).flatten()
proxMask = proximity > tol
return mus[proxMask]
class GenericDistributedGreedyApproximant(GenericDistributedApproximant):
"""
ROM greedy interpolant computation for parametric problems
(ABSTRACT).
Args:
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;
- 'muBounds': list of bounds for parameter values; defaults to
[0, 1];
- 'S': number of starting training points; defaults to 2;
- 'sampler': sample point generator; defaults to uniform sampler on
muBounds;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'interactive': whether to interactively terminate greedy
algorithm; defaults to False;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement; defaults to 0.2;
- 'nTestPoints': number of test points; defaults to maxIter /
refinementRatio;
- 'trainSetGenerator': training sample points generator; defaults
to Chebyshev sampler within muBounds.
Defaults to empty dict.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
homogeneized: Whether to homogeneize Dirichlet BCs.
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;
- 'muBounds': list of bounds for parameter values;
- 'S': number of starting training points;
- 'sampler': sample point generator;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'interactive': whether to interactively terminate greedy
algorithm;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
POD: whether to compute POD of snapshots.
muBounds: list of bounds for parameter values.
S: number of test points.
sampler: Sample point generator.
greedyTol: uniform error tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
robustTol: tolerance for robust Pade' denominator management.
samplingEngine: Sampling engine.
estimatorNormEngine: Engine for estimator norm computation.
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.
"""
TOL_INSTABILITY = 1e-6
def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._addParametersToList(["greedyTol", "interactive", "maxIter",
"refinementRatio", "nTestPoints",
"trainSetGenerator"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneized = homogeneized, verbosity = verbosity,
timestamp = timestamp)
self._postInit()
@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,
["greedyTol", "interactive",
"maxIter", "refinementRatio",
"nTestPoints", "trainSetGenerator"],
True, True, baselevel = 1)
GenericDistributedApproximant.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
if "greedyTol" in keyList:
self.greedyTol = approxParameters["greedyTol"]
elif not hasattr(self, "_greedyTol") or self.greedyTol is None:
self.greedyTol = 1e-2
if "interactive" in keyList:
self.interactive = approxParameters["interactive"]
elif not hasattr(self, "interactive") or self.interactive is None:
self.interactive = False
if "maxIter" in keyList:
self.maxIter = approxParameters["maxIter"]
elif not hasattr(self, "_maxIter") or self.maxIter is None:
self.maxIter = 1e2
if "refinementRatio" in keyList:
self.refinementRatio = approxParameters["refinementRatio"]
elif (not hasattr(self, "_refinementRatio")
or self.refinementRatio is None):
self.refinementRatio = 0.2
if "nTestPoints" in keyList:
self.nTestPoints = approxParameters["nTestPoints"]
elif (not hasattr(self, "_nTestPoints")
or self.nTestPoints is None):
self.nTestPoints = np.int(np.ceil(self.maxIter
/ self.refinementRatio))
if "trainSetGenerator" in keyList:
self.trainSetGenerator = approxParameters["trainSetGenerator"]
elif (not hasattr(self, "_trainSetGenerator")
or self.trainSetGenerator is None):
from rrompy.utilities.parameter_sampling import QuadratureSampler
self.trainSetGenerator = QuadratureSampler(self.muBounds,
"CHEBYSHEV")
del QuadratureSampler
@property
def mus(self):
"""Value of mus."""
return self._mus
@mus.setter
def mus(self, mus):
self._mus = np.array(mus, dtype = np.complex)
@property
def greedyTol(self):
"""Value of greedyTol."""
return self._greedyTol
@greedyTol.setter
def greedyTol(self, greedyTol):
if greedyTol < 0:
raise RROMPyException("greedyTol must be non-negative.")
if hasattr(self, "_greedyTol") and self.greedyTol is not None:
greedyTolold = self.greedyTol
else:
greedyTolold = -1
self._greedyTol = greedyTol
self._approxParameters["greedyTol"] = self.greedyTol
if greedyTolold != self.greedyTol:
self.resetSamples()
@property
def maxIter(self):
"""Value of maxIter."""
return self._maxIter
@maxIter.setter
def maxIter(self, maxIter):
if maxIter <= 0: raise RROMPyException("maxIter must be positive.")
if hasattr(self, "_maxIter") and self.maxIter is not None:
maxIterold = self.maxIter
else:
maxIterold = -1
self._maxIter = maxIter
self._approxParameters["maxIter"] = self.maxIter
if maxIterold != self.maxIter:
self.resetSamples()
@property
def refinementRatio(self):
"""Value of refinementRatio."""
return self._refinementRatio
@refinementRatio.setter
def refinementRatio(self, refinementRatio):
if refinementRatio <= 0. or refinementRatio > 1.:
raise RROMPyException(("refinementRatio must be between 0 "
"(excluded) and 1."))
if (hasattr(self, "_refinementRatio")
and self.refinementRatio is not None):
refinementRatioold = self.refinementRatio
else:
refinementRatioold = -1
self._refinementRatio = refinementRatio
self._approxParameters["refinementRatio"] = self.refinementRatio
if refinementRatioold != self.refinementRatio:
self.resetSamples()
@property
def nTestPoints(self):
"""Value of nTestPoints."""
return self._nTestPoints
@nTestPoints.setter
def nTestPoints(self, nTestPoints):
if nTestPoints <= 0:
raise RROMPyException("nTestPoints must be positive.")
if not np.isclose(nTestPoints, np.int(nTestPoints)):
raise RROMPyException("nTestPoints must be an integer.")
nTestPoints = np.int(nTestPoints)
if hasattr(self, "_nTestPoints") and self.nTestPoints is not None:
nTestPointsold = self.nTestPoints
else:
nTestPointsold = -1
self._nTestPoints = nTestPoints
self._approxParameters["nTestPoints"] = self.nTestPoints
if nTestPointsold != self.nTestPoints:
self.resetSamples()
@property
def trainSetGenerator(self):
"""Value of trainSetGenerator."""
return self._trainSetGenerator
@trainSetGenerator.setter
def trainSetGenerator(self, trainSetGenerator):
if 'generatePoints' not in dir(trainSetGenerator):
raise RROMPyException("trainSetGenerator type not recognized.")
if (hasattr(self, '_trainSetGenerator')
and self.trainSetGenerator is not None):
trainSetGeneratorOld = self.trainSetGenerator
self._trainSetGenerator = trainSetGenerator
self._approxParameters["trainSetGenerator"] = self.trainSetGenerator
if (not 'trainSetGeneratorOld' in locals()
or trainSetGeneratorOld != self.trainSetGenerator):
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._mus = []
- def initEstimatorNormEngine(self, normEng = None):
+ def initEstimatorNormEngine(self, normEngn : normEng = None):
"""Initialize estimator norm engine."""
- if (normEng is not None or not hasattr(self, "estimatorNormEngine")
- or self.estimatorNormEngine is None):
- if normEng is None:
+ if (normEngn is not None or not hasattr(self, "estimatorNormEngine")
+ or self.estimatorNormEngine is None):
+ if normEngn is None:
if not hasattr(self.HFEngine, "energyNormMatrix"):
self.HFEngine.buildEnergyNormForm()
estimatorEnergyMatrix = self.HFEngine.energyNormMatrix
else:
- if hasattr(normEng, "buildEnergyNormForm"):
- if not hasattr(normEng, "energyNormMatrix"):
- normEng.buildEnergyNormForm()
- estimatorEnergyMatrix = normEng.energyNormMatrix
+ if hasattr(normEngn, "buildEnergyNormForm"):
+ if not hasattr(normEngn, "energyNormMatrix"):
+ normEngn.buildEnergyNormForm()
+ estimatorEnergyMatrix = normEngn.energyNormMatrix
else:
- estimatorEnergyMatrix = normEng
+ estimatorEnergyMatrix = normEngn
self.estimatorNormEngine = normEngine(estimatorEnergyMatrix)
def errorEstimator(self, mus:List[np.complex]) -> List[np.complex]:
"""
Standard residual-based error estimator with explicit residual
computation.
"""
self.setupApprox()
nmus = len(mus)
err = np.empty(nmus)
if self.HFEngine.nbs == 1:
RHS = self.getRHS(mus[0], homogeneized = self.homogeneized)
RHSNorm = self.estimatorNormEngine.norm(RHS)
for j in range(nmus):
res = self.getRes(mus[j], homogeneized = self.homogeneized)
err[j] = self.estimatorNormEngine.norm(res) / RHSNorm
else:
for j in range(nmus):
res = self.getRes(mus[j], homogeneized = self.homogeneized)
RHS = self.getRHS(mus[j], homogeneized = self.homogeneized)
err[j] = (self.estimatorNormEngine.norm(res)
/ self.estimatorNormEngine.norm(RHS))
return np.abs(err)
def getMaxErrorEstimator(self, mus:List[np.complex],
plot : bool = False) -> Tuple[Np1D, int, float]:
"""
Compute maximum of (and index of maximum of) error estimator over given
parameters.
"""
errorEstTest = self.errorEstimator(mus)
idxMaxEst = np.argmax(errorEstTest)
maxEst = errorEstTest[idxMaxEst]
if plot and not np.all(np.isinf(errorEstTest)):
from matplotlib import pyplot as plt
plt.figure()
plt.semilogy(np.real(mus), errorEstTest, 'k')
plt.semilogy(np.real(mus[[0, -1]]), [self.greedyTol] * 2, 'r--')
plt.semilogy(np.real(self.mus),
2. * self.greedyTol * np.ones(len(self.mus)), '*m')
plt.semilogy(np.real(mus[idxMaxEst]), maxEst, 'xr')
plt.grid()
plt.show()
plt.close()
return errorEstTest, idxMaxEst, maxEst
def greedyNextSample(self, muidx:int, plotEst : bool = False)\
-> Tuple[Np1D, int, float, complex]:
"""Compute next greedy snapshot of solution map."""
modeAssert(self._mode, message = "Cannot add greedy sample.")
mu = self.muTest[muidx]
self.muTest = np.delete(self.muTest, muidx)
if self.verbosity >= 2:
verbosityDepth("MAIN", ("Adding {}-th sample point at {} to "
"training set.").format(
self.samplingEngine.nsamples + 1, mu),
timestamp = self.timestamp)
self.mus = np.append(self.mus, mu)
self.samplingEngine.nextSample(mu, homogeneized = self.homogeneized)
errorEstTest, muidx, maxErrorEst = self.getMaxErrorEstimator(
self.muTest, plotEst)
return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
modeAssert(self._mode, message = "Cannot start greedy algorithm.")
self.computeScaleFactor()
if self.samplingEngine.samples is not None:
return
if self.verbosity >= 2:
verbosityDepth("INIT", "Starting computation of snapshots.",
timestamp = self.timestamp)
self.resetSamples()
self.mus, _ = self.trainSetGenerator.generatePoints(self.S)
muTestBase, _ = self.sampler.generatePoints(self.nTestPoints)
for j in range(len(self.mus) - 1):
if self.verbosity >= 2:
verbosityDepth("MAIN",
("Adding {}-th sample point at {} to training "
"set.").format(self.samplingEngine.nsamples+ 1,
self.mus[j]),
timestamp = self.timestamp)
self.samplingEngine.nextSample(self.mus[j],
homogeneized = self.homogeneized)
muTestBase = np.sort(pruneSamples(muTestBase, self.mus[: -1],
1e-10 * self.scaleFactor))
self.muTest = np.empty(len(muTestBase) + 1, dtype = np.complex)
self.muTest[: -1] = muTestBase
self.muTest[-1] = self.mus[-1]
self.mus = self.mus[: -1]
def _enrichTestSet(self, nTest:int):
"""Add extra elements to test set."""
muTestExtra, _ = self.sampler.generatePoints(2 * nTest)
muTestExtra = pruneSamples(muTestExtra,
np.append(self.mus, self.muTest),
1e-10 * self.scaleFactor)
muTestNew = np.empty(len(self.muTest) + len(muTestExtra),
dtype = np.complex)
muTestNew[: len(self.muTest)] = self.muTest
muTestNew[len(self.muTest) :] = muTestExtra
self.muTest = np.sort(muTestNew)
if self.verbosity >= 5:
verbosityDepth("MAIN", "Enriching test set by {} elements.".format(
len(muTestExtra)),
timestamp = self.timestamp)
def greedy(self, plotEst : bool = False):
"""Compute greedy snapshots of solution map."""
modeAssert(self._mode, message = "Cannot start greedy algorithm.")
if self.samplingEngine.samples is not None:
return
self._preliminaryTraining()
nTest = self.nTestPoints
errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(-1,
plotEst)
if self.verbosity >= 2:
verbosityDepth("MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(maxErrorEst),
timestamp = self.timestamp)
trainedModelOld = copy(self.trainedModel)
while (self.samplingEngine.nsamples < self.maxIter
and maxErrorEst > self.greedyTol):
if (1. - self.refinementRatio) * nTest > len(self.muTest):
self._enrichTestSet(nTest)
nTest = len(self.muTest)
muTestOld, maxErrorEstOld = self.muTest, maxErrorEst
errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(
muidx, plotEst)
if self.verbosity >= 2:
verbosityDepth("MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(maxErrorEst),
timestamp = self.timestamp)
if (np.isnan(maxErrorEst) or np.isinf(maxErrorEst)
or maxErrorEstOld < maxErrorEst * self.TOL_INSTABILITY):
RROMPyWarning(("Instability in a posteriori estimator. "
"Starting preemptive greedy loop termination."))
maxErrorEst = maxErrorEstOld
self.muTest = muTestOld
self.mus = self.mus[:-1]
self.samplingEngine.popSample()
self.trainedModel.data = copy(trainedModelOld.data)
break
trainedModelOld.data = copy(self.trainedModel.data)
if (self.interactive and maxErrorEst <= self.greedyTol):
verbosityDepth("MAIN", ("Required tolerance {} achieved. Want "
"to decrease greedyTol and continue? "
"Y/N").format(self.greedyTol),
timestamp = self.timestamp, end = "")
increasemaxIter = input()
if increasemaxIter.upper() == "Y":
verbosityDepth("MAIN", "Reducing value of greedyTol...",
timestamp = self.timestamp)
while maxErrorEst <= self._greedyTol:
self._greedyTol *= .5
if (self.interactive
and self.samplingEngine.nsamples >= self.maxIter):
verbosityDepth("MAIN", ("Maximum number of iterations {} "
"reached. Want to increase maxIter "
"and continue? Y/N").format(
self.maxIter),
timestamp = self.timestamp, end = "")
increasemaxIter = input()
if increasemaxIter.upper() == "Y":
verbosityDepth("MAIN", "Doubling value of maxIter...",
timestamp = self.timestamp)
self._maxIter *= 2
if self.verbosity >= 2:
verbosityDepth("DEL", ("Done computing snapshots (final snapshot "
"count: {}).").format(
self.samplingEngine.nsamples),
timestamp = self.timestamp)
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return (super().checkComputedApprox()
and len(self.mus) == self.trainedModel.data.projMat.shape[1])
def assembleReducedResidualGramian(self, pMat:Np2D):
"""
Build residual gramian of reduced linear system through projections.
"""
self.initEstimatorNormEngine()
if (not hasattr(self.trainedModel.data, "gramian")
or self.trainedModel.data.gramian is None):
gramian = self.estimatorNormEngine.innerProduct(pMat, pMat)
else:
Sold = self.trainedModel.data.gramian.shape[0]
S = len(self.mus)
if Sold > S:
gramian = self.trainedModel.data.gramian[: S, : S]
else:
gramian = np.empty((S, S), dtype = np.complex)
gramian[: Sold, : Sold] = self.trainedModel.data.gramian
gramian[: Sold, Sold :] = (
self.estimatorNormEngine.innerProduct(pMat[:, Sold :],
pMat[:, : Sold]))
gramian[Sold :, : Sold] = gramian[: Sold, Sold :].T.conj()
gramian[Sold :, Sold :] = (
self.estimatorNormEngine.innerProduct(pMat[:, Sold :],
pMat[:, Sold :]))
self.trainedModel.data.gramian = gramian
def assembleReducedResidualBlocksbb(self, bs:List[Np1D], pMat:Np2D,
scaling : float = 1.):
"""
Build blocks (of type bb) of reduced linear system through projections.
"""
self.initEstimatorNormEngine()
nbs = len(bs)
if (not hasattr(self.trainedModel.data, "resbb")
or self.trainedModel.data.resbb is None):
resbb = np.empty((nbs, nbs), dtype = np.complex)
for i in range(nbs):
Mbi = scaling ** i * bs[i]
resbb[i, i] = self.estimatorNormEngine.innerProduct(Mbi, Mbi)
for j in range(i):
Mbj = scaling ** j * bs[j]
resbb[i, j] = self.estimatorNormEngine.innerProduct(Mbj,
Mbi)
for i in range(nbs):
for j in range(i + 1, nbs):
resbb[i, j] = resbb[j, i].conj()
self.trainedModel.data.resbb = resbb
def assembleReducedResidualBlocksAb(self, As:List[Np2D], bs:List[Np1D],
pMat:Np2D, scaling : float = 1.):
"""
Build blocks (of type Ab) of reduced linear system through projections.
"""
self.initEstimatorNormEngine()
nAs = len(As)
nbs = len(bs)
S = len(self.mus)
if (not hasattr(self.trainedModel.data, "resAb")
or self.trainedModel.data.resAb is None):
resAb = np.empty((nbs, S, nAs), dtype = np.complex)
for j in range(nAs):
MAj = scaling ** (j + 1) * As[j].dot(pMat)
for i in range(nbs):
Mbi = scaling ** (i + 1) * bs[i]
resAb[i, :, j] = self.estimatorNormEngine.innerProduct(MAj,
Mbi)
else:
Sold = self.trainedModel.data.resAb.shape[1]
if Sold == S: return
if Sold > S:
resAb = self.trainedModel.data.resAb[:, : S, :]
else:
resAb = np.empty((nbs, S, nAs), dtype = np.complex)
resAb[:, : Sold, :] = self.trainedModel.data.resAb
for j in range(nAs):
MAj = scaling ** (j + 1) * As[j].dot(pMat[:, Sold :])
for i in range(nbs):
Mbi = scaling ** (i + 1) * bs[i]
resAb[i, Sold :, j] = (
self.estimatorNormEngine.innerProduct(MAj, Mbi))
self.trainedModel.data.resAb = resAb
def assembleReducedResidualBlocksAA(self, As:List[Np2D], pMat:Np2D,
scaling : float = 1.,
basic : bool = False):
"""
Build blocks (of type AA) of reduced linear system through projections.
"""
self.initEstimatorNormEngine()
nAs = len(As)
S = len(self.mus)
if (not hasattr(self.trainedModel.data, "resAA")
or self.trainedModel.data.resAA is None):
if basic:
MAEnd = scaling ** nAs * As[-1].dot(pMat)
resAA = self.estimatorNormEngine.innerProduct(MAEnd, MAEnd)
else:
resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
for i in range(nAs):
MAi = scaling ** (i + 1) * As[i].dot(pMat)
resAA[:, i, :, i] = (
self.estimatorNormEngine.innerProduct(MAi, MAi))
for j in range(i):
MAj = scaling ** (j + 1) * As[j].dot(pMat)
resAA[:, i, :, j] = (
self.estimatorNormEngine.innerProduct(MAj, MAi))
for i in range(nAs):
for j in range(i + 1, nAs):
resAA[:, i, :, j] = resAA[:, j, :, i].T.conj()
else:
Sold = self.trainedModel.data.resAA.shape[0]
if Sold == S: return
if Sold > S:
if basic:
resAA = self.trainedModel.data.resAA[: S, : S]
else:
resAA = self.trainedModel.data.resAA[: S, :, : S, :]
else:
if basic:
resAA = np.empty((S, S), dtype = np.complex)
resAA[: Sold, : Sold] = self.trainedModel.data.resAA
MAi = scaling ** nAs * As[-1].dot(pMat)
resAA[: Sold, Sold :] = (
self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
MAi[:, : Sold]))
resAA[Sold :, : Sold] = resAA[: Sold, Sold :].T.conj()
resAA[Sold :, Sold :] = (
self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
MAi[:, Sold :]))
else:
resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
resAA[: Sold, :, : Sold, :] = self.trainedModel.data.resAA
for i in range(nAs):
MAi = scaling ** (i + 1) * As[i].dot(pMat)
resAA[: Sold, i, Sold :, i] = (
self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
MAi[:, : Sold]))
resAA[Sold :, i, : Sold, i] = resAA[: Sold, i,
Sold :, i].T.conj()
resAA[Sold :, i, Sold :, i] = (
self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
MAi[:, Sold :]))
for j in range(i):
MAj = scaling ** (j + 1) * As[j].dot(pMat)
resAA[: Sold, i, Sold :, j] = (
self.estimatorNormEngine.innerProduct(MAj[:, Sold :],
MAi[:, : Sold]))
resAA[Sold :, i, : Sold, j] = (
self.estimatorNormEngine.innerProduct(MAj[:, : Sold],
MAi[:, Sold :]))
resAA[Sold :, i, Sold :, j] = (
self.estimatorNormEngine.innerProduct(MAj[:, Sold :],
MAi[:, Sold :]))
for i in range(nAs):
for j in range(i + 1, nAs):
resAA[: Sold, i, Sold :, j] = (
resAA[Sold :, j, : Sold, i].T.conj())
resAA[Sold :, i, : Sold, j] = (
resAA[: Sold, j, Sold :, i].T.conj())
resAA[Sold :, i, Sold :, j] = (
resAA[Sold :, j, Sold :, i].T.conj())
self.trainedModel.data.resAA = resAA
diff --git a/rrompy/solver/__init__.py b/rrompy/solver/__init__.py
index ff83f0c..391cfa5 100644
--- a/rrompy/solver/__init__.py
+++ b/rrompy/solver/__init__.py
@@ -1,28 +1,34 @@
# 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 .linear_solver import RROMPyLinearSolvers, setupSolver
+from .norm_utilities import (Np2DLike, Np2DLikeInv, Np2DLikeInvLowRank,
+ normEngine)
__all__ = [
'RROMPyLinearSolvers',
- 'setupSolver'
+ 'setupSolver',
+ 'Np2DLike',
+ 'Np2DLikeInv',
+ 'Np2DLikeInvLowRank',
+ 'normEngine'
]
diff --git a/rrompy/utilities/fenics/__init__.py b/rrompy/solver/fenics/__init__.py
similarity index 72%
rename from rrompy/utilities/fenics/__init__.py
rename to rrompy/solver/fenics/__init__.py
index 6b2e26a..1e1a506 100644
--- a/rrompy/utilities/fenics/__init__.py
+++ b/rrompy/solver/fenics/__init__.py
@@ -1,34 +1,38 @@
# 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 .fenics_constants import fenZERO, fenZEROS, fenONE, fenONES, bdrTrue, bdrFalse
-from .fenics_norms import L2NormMatrix, H1NormMatrix, elasticNormMatrix
+from .fenics_constants import (fenZERO, fenZEROS, fenONE, fenONES,
+ bdrTrue, bdrFalse)
+from .fenics_norms import (L2NormMatrix, H1NormMatrix, Hminus1NormMatrix,
+ elasticNormMatrix, elasticDualNormMatrix)
__all__ = [
'fenZERO',
'fenZEROS',
'fenONE',
'fenONES',
'bdrTrue',
'bdrFalse',
'L2NormMatrix',
'H1NormMatrix',
- 'elasticNormMatrix'
+ 'Hminus1NormMatrix',
+ 'elasticNormMatrix',
+ 'elasticDualNormMatrix'
]
diff --git a/rrompy/utilities/fenics/fenics_constants.py b/rrompy/solver/fenics/fenics_constants.py
similarity index 100%
rename from rrompy/utilities/fenics/fenics_constants.py
rename to rrompy/solver/fenics/fenics_constants.py
diff --git a/rrompy/solver/fenics/fenics_norms.py b/rrompy/solver/fenics/fenics_norms.py
new file mode 100644
index 0000000..a4e36c8
--- /dev/null
+++ b/rrompy/solver/fenics/fenics_norms.py
@@ -0,0 +1,76 @@
+# 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 scipy.sparse import csr_matrix
+import fenics as fen
+from rrompy.utilities.base.types import Np2D, FenFunc, DictAny, FenFuncSpace
+from rrompy.solver import Np2DLikeInv, Np2DLikeInvLowRank
+
+__all__ = ['L2NormMatrix', 'H1NormMatrix', 'Hminus1NormMatrix',
+ 'elasticNormMatrix', 'elasticDualNormMatrix']
+
+def _fen2sp(expr):
+ matFen = fen.as_backend_type(fen.assemble(expr)).mat()
+ return csr_matrix(matFen.getValuesCSR()[::-1], shape = matFen.size)
+
+def L2NormMatrix(V:FenFuncSpace, r_ : FenFunc = 1.) -> Np2D:
+ u = fen.TrialFunction(V)
+ v = fen.TestFunction(V)
+ return _fen2sp(r_ * fen.dot(u, v) * fen.dx)
+
+def H1NormMatrix(V:FenFuncSpace, w : float = 0., r_ : FenFunc = 1.,
+ a_ : FenFunc = 1.) -> Np2D:
+ u = fen.TrialFunction(V)
+ v = fen.TestFunction(V)
+ return _fen2sp((w * r_ * fen.dot(u, v)
+ + fen.dot(a_ * fen.grad(u), fen.grad(v))) * fen.dx)
+
+def Hminus1NormMatrix(V:FenFuncSpace, w : float = 0., r_ : FenFunc = 1.,
+ a_ : FenFunc = 1., solverType : str = "SPSOLVE",
+ solverArgs : DictAny = {}, compressRank : int = None,
+ compressOversampling : int = 10,
+ compressSeed : int = 420) -> Np2D:
+ if compressRank is None:
+ return Np2DLikeInv(H1NormMatrix(V, w, r_, a_), L2NormMatrix(V, r_),
+ solverType, solverArgs)
+ return Np2DLikeInvLowRank(H1NormMatrix(V, w, r_, a_), L2NormMatrix(V, r_),
+ solverType, solverArgs, compressRank,
+ compressOversampling, compressSeed)
+
+def elasticNormMatrix(V:FenFuncSpace, l_:FenFunc, m_:FenFunc,
+ w : float = 0., r_ : FenFunc = 1.) -> Np2D:
+ u = fen.TrialFunction(V)
+ v = fen.TestFunction(V)
+ epsilon = lambda f: 0.5 * (fen.grad(f) + fen.nabla_grad(f))
+ sigma = (l_ * fen.div(u) * fen.Identity(u.geometric_dimension())
+ + 2. * m_ * epsilon(u))
+ return _fen2sp((w * r_ * fen.dot(u, v)
+ + fen.inner(sigma, epsilon(v))) * fen.dx)
+
+def elasticDualNormMatrix(V:FenFuncSpace, l_:FenFunc, m_:FenFunc,
+ w : float = 0., solverType : str = "SPSOLVE",
+ solverArgs : DictAny = {}, r_ : FenFunc = 1.,
+ compressRank : int = None,
+ compressOversampling : int = 10,
+ compressSeed : int = 420) -> Np2D:
+ if compressRank is None:
+ return Np2DLikeInv(elasticNormMatrix(V, l_, m_, w, r_),
+ L2NormMatrix(V, r_), solverType, solverArgs)
+ return Np2DLikeInvLowRank(elasticNormMatrix(V, l_, m_, w, r_),
+ L2NormMatrix(V, r_), solverType, solverArgs,
+ compressRank, compressOversampling, compressSeed)
diff --git a/rrompy/solver/norm_utilities.py b/rrompy/solver/norm_utilities.py
new file mode 100644
index 0000000..a86b3ba
--- /dev/null
+++ b/rrompy/solver/norm_utilities.py
@@ -0,0 +1,71 @@
+# 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 deepcopy as copy
+from rrompy.utilities.base.types import Np1D, Np2D, DictAny
+from rrompy.solver.linear_solver import setupSolver
+from rrompy.utilities.exception_manager import RROMPyException
+
+__all__ = ['Np2DLike', 'Np2DLikeInv', 'Np2DLikeInvLowRank', 'normEngine']
+
+@abstractmethod
+class Np2DLike:
+ def dot(self, u:Np2D) -> Np2D:
+ pass
+
+class Np2DLikeInv(Np2DLike):
+ def __init__(self, K:Np2D, M:Np2D, solverType:str, solverArgs:DictAny):
+ self.K, self.M, self.MH = K, M, M.T.conj()
+ self.solver, self.solverArgs = setupSolver(solverType, solverArgs)
+ def dot(self, u:Np2D) -> Np2D:
+ return self.MH.dot(self.solver(self.K, self.M.dot(u), self.solverArgs))
+
+class Np2DLikeInvLowRank(Np2DLike):
+ def __init__(self, K:Np2D, M:Np2D, solverType:str, solverArgs:DictAny,
+ rank:int, oversampling : int = 10, seed : int = 420):
+ if rank > M.shape[1]:
+ RROMPyException(("Cannot select compressed rank larger than "
+ "original size."))
+ if oversampling < 0:
+ RROMPyException("Oversampling parameter must be positive.")
+ HF = Np2DLikeInv(K, M, solverType, solverArgs)
+ np.random.seed(seed)
+ xs = np.random.randn(M.shape[1], rank + oversampling)
+ samples = HF.dot(xs)
+ Q, _ = np.linalg.qr(samples, mode = "reduced")
+ R = HF.dot(Q).T.conj() # assuming HF (i.e. K) hermitian...
+ U, s, Vh = np.linalg.svd(R)
+ self.L = Q.dot(U[:, : rank]) * s[: rank]
+ self.R = Vh[: rank, :]
+ def dot(self, u:Np2D) -> Np2D:
+ return self.L.dot(self.R.dot(u))
+
+class normEngine:
+ def __init__(self, energyNormMatrix:Np2D):
+ self.energyNormMatrix = copy(energyNormMatrix)
+
+ def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False) -> Np2D:
+ if onlyDiag:
+ return np.sum(self.energyNormMatrix.dot(u) * v.conj(), axis = 0)
+ return v.T.conj().dot(self.energyNormMatrix.dot(u))
+
+ def norm(self, u:Np2D) -> Np1D:
+ return np.abs(self.innerProduct(u, u, onlyDiag = True)) ** .5
+
diff --git a/rrompy/utilities/base/__init__.py b/rrompy/utilities/base/__init__.py
index c236f15..4e3ae8f 100644
--- a/rrompy/utilities/base/__init__.py
+++ b/rrompy/utilities/base/__init__.py
@@ -1,46 +1,48 @@
# 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 .find_dict_str_key import findDictStrKey
from .get_new_filename import getNewFilename
from .pickle_utilities import pickleDump, pickleLoad
from .purge_dict import purgeDict
from .purge_list import purgeList
from .number_theory import (squareResonances, primeFactorize,
getLowestPrimeFactor)
from .sobol import sobolGenerate
from .low_discrepancy import vanderCorput, lowDiscrepancy
from . import types as Types
from .verbosity_depth import verbosityDepth
__all__ = [
'findDictStrKey',
'getNewFilename',
'pickleDump',
'pickleLoad',
'purgeDict',
'purgeList',
'squareResonances',
'primeFactorize',
'getLowestPrimeFactor',
'sobolGenerate',
+ 'vanderCorput',
+ 'lowDiscrepancy',
'Types',
'verbosityDepth'
]
diff --git a/rrompy/utilities/base/types.py b/rrompy/utilities/base/types.py
index 4c8b273..15fdaef 100644
--- a/rrompy/utilities/base/types.py
+++ b/rrompy/utilities/base/types.py
@@ -1,52 +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','ScOp','Np1D','Np2D','Np1DLst',
'N2FSExpr','FenExpr','FenFunc','FenFuncSpace','HFEng','ROMEng',
- 'sampleEng','GenExpr','strLst','BfSExpr']
+ 'sampleEng','normEng','GenExpr','strLst','BfSExpr']
# ANY
TupleAny = Tuple[Any]
ListAny = List[Any]
DictAny = Dict[Any, Any]
# SCIPY
ScOp = TypeVar("Scipy sparse matrix for space operator")
# NUMPY
Np1D = TypeVar("NumPy 1D array")
-Np2D = TypeVar("NumPy 2D array")
+Np2D = TypeVar("NumPy 2D array-like")
Np1DLst = TypeVar("NumPy 1D array or list of NumPy 1D array")
N2FSExpr = TypeVar("NumPy 2D array, float or str")
# FENICS
FenExpr = TypeVar("FEniCS expression")
FenFunc = TypeVar("FEniCS function")
FenFuncSpace = TypeVar("FEniCS function space")
# ENGINES
HFEng = TypeVar("High fidelity engine")
ROMEng = TypeVar("ROM engine")
sampleEng = TypeVar("Sampling engine")
+normEng = TypeVar("Norm engine")
# OTHERS
GenExpr = TypeVar("Generic expression")
strLst = TypeVar("str or list of str")
BfSExpr = TypeVar("Boolean function or string")
diff --git a/rrompy/utilities/fenics/fenics_norms.py b/rrompy/utilities/fenics/fenics_norms.py
deleted file mode 100644
index f2ec21f..0000000
--- a/rrompy/utilities/fenics/fenics_norms.py
+++ /dev/null
@@ -1,48 +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 fenics as fen
-from rrompy.utilities.base.types import Np2D, FenFunc, FenFuncSpace
-from scipy.sparse import csr_matrix
-
-__all__ = ['L2NormMatrix', 'H1NormMatrix', 'elasticNormMatrix']
-
-def _fen2sp(expr):
- matFen = fen.as_backend_type(fen.assemble(expr)).mat()
- return csr_matrix(matFen.getValuesCSR()[::-1], shape = matFen.size)
-
-def L2NormMatrix(V:FenFuncSpace) -> Np2D:
- u = fen.TrialFunction(V)
- v = fen.TestFunction(V)
- return _fen2sp(fen.dot(u, v) * fen.dx)
-
-def H1NormMatrix(V:FenFuncSpace, w : FenFunc = 0.) -> Np2D:
- u = fen.TrialFunction(V)
- v = fen.TestFunction(V)
- return _fen2sp((w * fen.dot(u, v)
- + fen.dot(fen.grad(u), fen.grad(v))) * fen.dx)
-
-def elasticNormMatrix(V:FenFuncSpace, l_:FenFunc, m_:FenFunc,
- w : FenFunc = 0.) -> Np2D:
- u = fen.TrialFunction(V)
- v = fen.TestFunction(V)
- epsilon = lambda f: 0.5 * (fen.grad(f) + fen.nabla_grad(f))
- sigma = (l_ * fen.div(u) * fen.Identity(u.geometric_dimension())
- + 2. * m_ * epsilon(u))
- return _fen2sp((w * fen.dot(u, v) + fen.inner(sigma, epsilon(v))) * fen.dx)
-
diff --git a/tests/test_1_utilities/fenics_const.py b/tests/test_1_utilities/fenics_const.py
index 54e834b..96d611e 100644
--- a/tests/test_1_utilities/fenics_const.py
+++ b/tests/test_1_utilities/fenics_const.py
@@ -1,23 +1,20 @@
# 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 .
#
def test_fenics_loads():
- from rrompy.utilities.fenics import fenZERO, fenZEROS, fenONE, fenONES
- from rrompy.utilities.fenics import (L2NormMatrix, H1NormMatrix,
- elasticNormMatrix)
-
+ from rrompy.solver.fenics import fenZERO, fenZEROS, fenONE, fenONES
diff --git a/tests/test_1_utilities/fenics_norms.py b/tests/test_1_utilities/fenics_norms.py
index bc1b4f0..8340afb 100644
--- a/tests/test_1_utilities/fenics_norms.py
+++ b/tests/test_1_utilities/fenics_norms.py
@@ -1,66 +1,95 @@
# 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.utilities.fenics import (L2NormMatrix, H1NormMatrix,
- elasticNormMatrix)
+from rrompy.solver.fenics import (L2NormMatrix, H1NormMatrix,
+ Hminus1NormMatrix, elasticNormMatrix,
+ elasticDualNormMatrix)
def test_fenics_L2():
V = fen.FunctionSpace(fen.UnitSquareMesh(3, 3), "P", 3)
u = fen.interpolate(fen.Constant(3.), V)
v = fen.interpolate(fen.Expression("x[0]*x[0]+x[1]", degree = 2), V)
uv = u.vector()[:]
vv = v.vector()[:]
mass = L2NormMatrix(V)
inner = fen.assemble(fen.dot(u, v) * fen.dx)
assert np.isclose(uv.T.dot(mass.dot(vv)), 2.5, rtol = 1e-8)
assert np.isclose(inner, 2.5, rtol = 1e-8)
def test_fenics_H1():
V = fen.FunctionSpace(fen.UnitSquareMesh(3, 3), "P", 2)
u = fen.interpolate(fen.Expression("x[0]+exp(x[1])", degree = 4), V)
v = fen.interpolate(fen.Expression("x[0]*x[0]+x[1]", degree = 2), V)
uv = u.vector()[:]
vv = v.vector()[:]
stiffness = H1NormMatrix(V)
helmholtz = H1NormMatrix(V, 12)
inners = fen.assemble(fen.dot(fen.grad(u), fen.grad(v)) * fen.dx)
innerh = 12 * fen.assemble(fen.dot(u, v) * fen.dx)
assert np.isclose(uv.T.dot(stiffness.dot(vv)), np.exp(1), rtol = 1e-6)
assert np.isclose(inners, np.exp(1), rtol = 1e-6)
assert np.isclose(uv.T.dot(helmholtz.dot(vv)), 5 * np.exp(1) + 14,
rtol = 1e-3)
assert np.isclose(inners + innerh, 5 * np.exp(1) + 14, rtol = 1e-3)
+def test_fenics_Hminus1():
+ V = fen.FunctionSpace(fen.UnitSquareMesh(3, 3), "P", 2)
+ u = fen.interpolate(fen.Expression("x[0]+exp(x[1])", degree = 4), V)
+ v = fen.interpolate(fen.Expression("x[0]*x[0]+x[1]", degree = 2), V)
+ uv = u.vector()[:]
+ vv = v.vector()[:]
+ energyFull = Hminus1NormMatrix(V, 12)
+ energyLR = Hminus1NormMatrix(V, 12, compressRank = 20)
+ assert np.isclose(uv.T.dot(energyFull.dot(vv)),
+ uv.T.dot(energyLR.dot(vv)), rtol = 1e-2)
+ assert np.isclose(uv.T.dot(energyFull.dot(vv)), .1641618355, rtol = 1e-6)
+
def test_fenics_elastic():
V = fen.VectorFunctionSpace(fen.UnitCubeMesh(5, 5, 5), "P", 1)
l_ = 1.
m_ = fen.Expression(".5*x[0]+1.", degree = 1)
u = fen.interpolate(fen.Expression(("exp(x[1])", "x[0]-x[2]", "3."),
degree = 4), V)
v = fen.interpolate(fen.Expression(("x[0]*x[0]+x[2]", "1.", "-1. * x[1]"),
degree = 2), V)
uv = u.vector()[:]
vv = v.vector()[:]
energyMat = elasticNormMatrix(V, l_, m_, 10)
epsilon = lambda f: 0.5 * (fen.grad(f) + fen.nabla_grad(f))
sigma = (l_ * fen.div(u) * fen.Identity(3) + 2. * m_ * epsilon(u))
energy = fen.assemble((10 * fen.dot(u, v)
+ fen.inner(sigma, epsilon(v))) * fen.dx)
assert np.isclose(uv.T.dot(energyMat.dot(vv)), energy, rtol = 1e-8)
+
+def test_fenics_elastic_dual():
+ V = fen.VectorFunctionSpace(fen.UnitCubeMesh(5, 5, 5), "P", 1)
+ l_ = 1.
+ m_ = fen.Expression(".5*x[0]+1.", degree = 1)
+ u = fen.interpolate(fen.Expression(("exp(x[1])", "x[0]-x[2]", "3."),
+ degree = 4), V)
+ v = fen.interpolate(fen.Expression(("x[0]*x[0]+x[2]", "1.", "-1. * x[1]"),
+ degree = 2), V)
+ uv = u.vector()[:]
+ vv = v.vector()[:]
+ energyFull = elasticDualNormMatrix(V, l_, m_, 10)
+ energyLR = elasticDualNormMatrix(V, l_, m_, 10, compressRank = 50)
+ assert np.isclose(uv.T.dot(energyFull.dot(vv)),
+ uv.T.dot(energyLR.dot(vv)), rtol = 1e-1)
+ assert np.isclose(uv.T.dot(energyFull.dot(vv)), -.00804628936, rtol = 1e-6)
diff --git a/tests/test_2_hfengines/helmholtz_external.py b/tests/test_2_hfengines/helmholtz_external.py
index e297ebc..4ff624a 100644
--- a/tests/test_2_hfengines/helmholtz_external.py
+++ b/tests/test_2_hfengines/helmholtz_external.py
@@ -1,45 +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.hfengines.linear_problem import (
HelmholtzCavityScatteringProblemEngine, HelmholtzBoxScatteringProblemEngine)
def test_helmholtz_square_scattering():
solver = HelmholtzCavityScatteringProblemEngine(kappa = 4, gamma = 2.,
n = 20, verbosity = 0)
mu = 5
uh = solver.solve(mu)
assert np.isclose(solver.norm(uh), 20.719752682674923, rtol = 1e-5)
assert np.isclose(solver.norm(solver.residual(uh, mu)), 4.25056407e-13,
rtol = 1e-1)
+def test_helmholtz_scattering_copy(capsys):
+ solver1 = HelmholtzCavityScatteringProblemEngine(kappa = 4, gamma = 2.,
+ n = 20, verbosity = 0)
+ mu = 5
+ uh1 = solver1.solve(mu)
+ solver2 = HelmholtzCavityScatteringProblemEngine(kappa = 4, gamma = 2.,
+ n = 20, verbosity = 100)
+ assert solver1.As[0] is not None and solver1.bs[0] is not None
+ assert solver2.As[0] is None and solver2.bs[0] is None
+ solver2.setAs(solver1.As)
+ solver2.setbs(solver1.bs)
+ uh2 = solver2.solve(mu)
+ assert np.allclose(uh1, uh2, rtol = 1e-8)
+
+ out, err = capsys.readouterr()
+ assert ("Assembling operator term" not in out
+ and "Assembling forcing term" not in out)
+ assert len(err) == 0
+
def test_helmholtz_box_scattering():
solver = HelmholtzBoxScatteringProblemEngine(R = 2, kappa = 10.,
theta = np.pi * 30 / 180, n = 20, verbosity = 0)
mu = 15
uh = solver.solve(mu)
solver.plotmesh(show = False, figsize = (7, 7))
assert np.isclose(solver.norm(uh), 63.98946657389119, rtol = 1e-5)
assert np.isclose(solver.norm(solver.residual(uh, mu)), 9.62989935e-13,
rtol = 1e-1)
from matplotlib import pyplot as plt
plt.close('all')
diff --git a/tests/test_3_reduction_methods/rational_interpolant_greedy.py b/tests/test_3_reduction_methods/rational_interpolant_greedy.py
index d4ddc6b..cea03ee 100644
--- a/tests/test_3_reduction_methods/rational_interpolant_greedy.py
+++ b/tests/test_3_reduction_methods/rational_interpolant_greedy.py
@@ -1,67 +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 matrix_fft import matrixFFT
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as RIG
def test_lax_tolerance(capsys):
solver = matrixFFT()
params = {"POD": True, "muBounds": [1.5, 6.5], "S": 4,
"polybasis": "CHEBYSHEV", "greedyTol": 1e-2,
"errorEstimatorKind": "bare"}
approx = RIG(solver, 4., params, verbosity = 10)
approx.greedy()
out, err = capsys.readouterr()
assert "Done computing snapshots (final snapshot count: 10)." in out
assert len(err) == 0
assert np.isclose(approx.normErr(0), .0077041389, rtol = 1e-3)
def test_samples_at_poles():
solver = matrixFFT()
params = {"POD": True, "muBounds": [1.5, 6.5], "S": 4, "nTestPoints": 100,
"polybasis": "CHEBYSHEV", "greedyTol": 1e-5,
"errorEstimatorKind": "exact"}
approx = RIG(solver, 4., params, verbosity = 0)
approx.greedy()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu) / (1e-15 + approx.normHF(mu)),
0., atol = 1e-4)
poles = approx.getPoles()
for lambda_ in range(2, 7):
assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-3)
assert np.isclose(np.min(np.abs(approx.mus - lambda_)), 0.,
atol = 1e-1)
def test_maxIter():
solver = matrixFFT()
params = {"POD": True, "muBounds": [1.5, 6.5], "S": 5, "nTestPoints": 500,
"polybasis": "CHEBYSHEV", "greedyTol": 1e-6, "maxIter": 10,
"errorEstimatorKind": "basic"}
approx = RIG(solver, 4., params, verbosity = 0)
approx.input = lambda: "N"
approx.greedy()
assert len(approx.mus) == 10
_, _, maxEst = approx.getMaxErrorEstimator(approx.muTest)
assert maxEst > 1e-6
+def test_load_copy(capsys):
+ mu = 3.
+ solver = matrixFFT()
+ params = {"POD": True, "muBounds": [1.5, 6.5], "S": 4, "nTestPoints": 100,
+ "polybasis": "CHEBYSHEV", "greedyTol": 1e-5,
+ "errorEstimatorKind": "exact"}
+ approx1 = RIG(solver, 4., params, verbosity = 100)
+ approx1.greedy()
+ err1 = approx1.normErr(mu)
+ out, err = capsys.readouterr()
+ assert "Solving HF model for mu =" in out
+ assert len(err) == 0
+ approx2 = RIG(solver, 4., params, verbosity = 100)
+ approx2.setApprox(approx1)
+ approx2.setHF(mu, approx1.uHF)
+ err2 = approx2.normErr(mu)
+ out, err = capsys.readouterr()
+ assert "Solving HF model for mu =" not in out
+ assert len(err) == 0
+ assert np.isclose(err1, err2, rtol = 1e-10)
diff --git a/tests/test_3_reduction_methods/rb_centered.py b/tests/test_3_reduction_methods/rb_centered.py
index dcc0abe..44a04b8 100644
--- a/tests/test_3_reduction_methods/rb_centered.py
+++ b/tests/test_3_reduction_methods/rb_centered.py
@@ -1,50 +1,69 @@
# 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 matrix_fft import matrixFFT
from rrompy.reduction_methods.centered import RBCentered as RBC
def test_R():
mu = 1.5
mu0 = 2. + 1.j
solver = matrixFFT()
uh = solver.solve(mu)
params = {"POD": True, "R": 5, "E": 10, "sampleType": "Krylov"}
approx = RBC(solver, mu0, params, verbosity = 0)
approx.setupApprox()
uhP = approx.getApprox(mu)
errP = approx.getErr(mu)
errNP = approx.normErr(mu)
assert np.allclose(np.abs(errP - (uhP - uh)), 0., rtol = 1e-3)
assert np.isclose(errNP, 0.023691832, rtol = 1e-1)
poles = approx.getPoles()
assert np.isclose(np.min(np.abs(poles - 2.)), 0., atol = 1e-4)
assert np.isclose(np.min(np.abs(poles - 1.)), 0., atol = 1e-2)
assert np.isclose(np.min(np.abs(poles - 3.)), 0., atol = 1e-2)
-
+
def test_moments():
mu0 = 2. + 1.j
solver = matrixFFT()
params = {"POD": True, "E": 10, "sampleType": "Krylov"}
approx = RBC(solver, mu0, params, verbosity = 0)
approx.setupApprox()
assert np.isclose(approx.normErr(mu0), 0., atol = 1e-10)
+def test_load_copy(capsys):
+ mu = 1.5
+ mu0 = 2. + 1.j
+ solver = matrixFFT()
+ params = {"POD": True, "E": 10, "sampleType": "Arnoldi"}
+ approx1 = RBC(solver, mu0, params, verbosity = 100)
+ approx1.setupApprox()
+ err1 = approx1.normErr(mu)
+ out, err = capsys.readouterr()
+ assert "Solving HF model for mu =" in out
+ assert len(err) == 0
+ approx2 = RBC(solver, mu0, params, verbosity = 100)
+ approx2.setApprox(approx1.trainedModel)
+ approx2.setHF(mu, approx1.uHF)
+ err2 = approx2.normErr(mu)
+ out, err = capsys.readouterr()
+ assert "Solving HF model for mu =" not in out
+ assert len(err) == 0
+ assert np.isclose(err1, err2, rtol = 1e-10)